Physics-based numerical models often depend on several parameters to close. Some of them can be expressed using established theoretical or empirical closure formulations. However, some others aggregate complex physical processes and are hence left as tuneable parameters, and can only be calibrated by trial and error. Yet, calibration data are not always available to do so, which prevents these models from being applied to wide ranges of laboratory or river flows. We hence propose a machine learning-based methodology to close any group of unclosed and correlated parameters, applied here to a two-phase/two-layer (2P2L) morphodynamical model. The methodology combines a numerical experiment with a known theoretical solution and machine learning. It is applied to the considered model to close two friction parameters for which generalizable and vastly acknowledged closure formulations lack in the literature. The resulting hybrid model, combining the original 2P2L model and the closure models, is tested against two laboratory dam break test cases. Despite excessive smoothness and underestimation of the concentration in sediment, the hybrid model performed similarly to other models from the literature requiring trial and error calibration and showed high stability and accuracy regarding the estimation of the water-sediment mixture's inertia.

  • A methodology to close tuneable parameters applicable to many morphodynamical models.

  • Methodology based on machine learning and self-generated data from a numerical experiment with theoretical solution.

  • Presentation of a two-phase/two-layer model for unsteady flows highly laden with sediment, used as a proof of concept of the methodology.

  • New model's performance is similar to evenly complex models requiring trial and error to calibrate tuneable parameters.

ANN

Artificial neural network

CM

Closure model

CV

Cross-validation

HM

Hybrid model

ML

Machine learning

MPM

Meyer-Peter & Müller

PBNM

Physics-based numerical model

RMSE

Root mean square error

2P2L

Two-phase/two-layer

The shallow-water equations constitute a physics-based numerical model (PBNM) widely used to numerically model river flows, either for reproduction or prediction purposes. For simple channel configurations with limited sediment transport and vegetation development, they can provide very accurate results (Soares-Frazão & Zech 2002). However, when these conditions are not met, they can be associated with other equations, e.g. the Exner equation (Soares-Frazão & Zech 2011; Meurice & Soares-Frazão 2020) or extended to account for specific processes, e.g. multi-phase flows or debris flows (Di Cristo et al. 2016; Martínez-Aranda et al. 2020). In particular, these extended equations often account for several parameters that need to be closed.

Many theoretical or empirical closure formulations have been developed to include physical processes in PBNMs, e.g. bedload transport capacity (Meyer-Peter & Müller 1948, MPM). Yet, some processes may be too complex to be described physically and implemented numerically, a fortiori at a reasonable computational cost. Ad-hoc parameters encompassing the intractable physics may constitute a satisfying solution, provided that sufficiently rich data can be used to calibrate them. Unfortunately, poorly investigated rivers hence fall out of the scope of these models.

Machine learning (ML) has recently become an increasingly powerful fluid dynamics tool (Vinuesa & Brunton 2022). It is mainly used for three purposes: accelerating direct numerical simulations, developing reduced-order models, and improving closure models (CMs), e.g. turbulence models. Zhu et al. (2022) classify its use into three categories: data-driven models fed with large amounts of data, knowledge-informed models that are physically constrained, and hybrid models (HMs) combining classical PBNMs with ML either to predict some result used by the PBNM, develop a new CM, or tune unclosed parameters (Goldstein et al. 2019). In morphodynamics specifically, ML has mainly been used to predict sediment transport, to build morphological models (e.g. shoreline profile), and as part of HMs (Goldstein et al. 2019).

In this paper, we propose a novel methodology which generates data from PBNMs and then uses them with ML to develop CMs for unclosed and correlated parameters of any morphodynamical PBNM in a hybrid-modelling approach. This methodology aims to get rid of trial and error for tuning parameters of PBNMs and to bridge the gap between complex physical processes and numerical implementation. To exemplify this procedure, we apply it to a two-phase/two-layer (2P2L) PBNM (Meurice & Soares-Frazão 2022) and find closure for static and dynamic friction coefficients for which CMs lack in the literature. We train and validate these CMs on data generated by a numerical experiment with known expected solutions and covering a large range of flows. To preserve the interpretability of the CMs, the ML method employed is a supervised stepwise weighted linear regression method. Furthermore, this method is combined with a k-fold cross-validation (CV), which suits the limited nature of the generated dataset, contrary to more advanced deep learning methods, less interpretable. Such explicit CMs can hence be easily implemented in the HM. The resulting HM not only provides results with a similar degree of accuracy as the original PBNM with optimized parameters does, but also expands the PBNM's applicability to a wider range of flow scenarios, including those where calibration data are scarce or unavailable. This study thus contributes to the ongoing efforts in leveraging ML for the advancement of hydromorphodynamical modelling.

First, the 2P2L model that is used to present the proposed methodology is described. Then, the numerical experiment and the methodology used to find a proper CM are explained. The CMs developed for two parameters are then presented. Finally, the resulting HM is applied to experimental test cases and conclusions are drawn.

Two-phase/two-layer model

The PBNM discussed in this paper and illustrated in Figure 1 is a 2D finite-volume morphodynamical model based on the work of Di Cristo et al. (2016). It was introduced by Meurice & Soares-Frazão (2022) and is described in more detail below. Not only can the model simulate water and bed level evolution but also that of the sediment concentration in the water column, split into two layers, referenced to as upper and lower. It is suited for transient flows with high suspension and bedload transport.
Figure 1

Schematic representation of the 2P2L model and its associated quantities.

Figure 1

Schematic representation of the 2P2L model and its associated quantities.

Close modal

The solid concentrations in both layers of the water column and are variable in time and space. Sediment and water mass exchanges are possible between the lower layer and the bed through the rate, ensuring a constant concentration in the bed. However, only sediment mass exchanges can occur between the two water column layers, through the rate.

The velocity of the solid phase in the lower layer may be different from that of the water while that of the solid phase of the upper layer is supposed to be equal to the velocity of the liquid phase . The latter is averaged over the whole water column to save computational time.

The water level , that of the interface between the two layers , and the bed level may also vary in time and space while the bedrock level is non-erodible and can only vary in space.

The 2P2L model counts four continuity equations and four momentum equations for there are two equations for each velocity and (Equations 1(d)–(1e)) since it is a 2D model:
formula
(1a)
formula
(1b)
formula
(1c)
formula
(1d)
formula
(1e)
formula
(1f)
where , , and are, respectively, the volumic fractions of the solid phases of the upper and lower layers and that of the liquid phase while , , and are the depths relative to the upper and lower layers and the whole water column, respectively. The porosity of the bed is referred to as and is such that . The constant densities of the liquid and wet solid phases are expressed as and , respectively, while the effective density of the solid phase is represented by ; g is the gravitational acceleration constant; and and are the friction stresses exerted, respectively, by the liquid and solid phases on the bed. The drag between the two phases is encapsulated in .

Equation (1) is similar to that obtained by Di Cristo et al. (2016), except for Equation (1e) which adds a pressure term due to the influence of sediment in suspension. Nevertheless, the influence of this term on the results is usually negligible. Some CMs used by Di Cristo et al. (2016) are the same as those used here, but some differ and are detailed in the Supplementary Material, Appendix A. Nevertheless, since this paper focuses on closing friction coefficients, these will be properly introduced hereafter.

Closure models of friction coefficients

According to Di Cristo et al. (2016), the friction stress exerted by the solid phase on the bed may be expressed as:
formula
(2)
with being the dynamic friction coefficient of the bed material and the collisional coefficient (see Di Cristo et al. (2016) for closure of ).

This formulation combines a Manning-like term accounting for the collisions between particles and a Mohr–Coulomb-like term depending on the weight of the particles acting on the bed and considered as a continuous medium rather than a granular one. Contrary to Di Cristo et al. (2016), we added to in the Mohr–Coulomb-like term to be consistent with our equations.

The resistance that the bed opposes to its shear can be modelled in a Mohr–Coulomb way as:
formula
(3)
with being the static friction coefficient of the bed material.

The value of the dynamic friction coefficient of the bed material is difficult to measure experimentally and its definition may vary (Al-Hashemi & Al-Amoudi 2018). Nevertheless, it should always be inferior to , which is often associated with the angle of repose of the granular material so that , even though Al-Hashemi & Al-Amoudi (2018) highlighted that and the static friction angle were not exactly equal in all conditions. Determining the relationship between the static and dynamic friction angles and is difficult and depends on many factors such as the water content, the grain size distribution, the shape of the grains, etc. Consequently, no empirical CM stands out as vastly recognized by the scientific community despite its continuous efforts (Buffington et al. 1992; Blau 2001; Combarros et al. 2014). Due to the complex granular constitution of bed materials, Al-Hashemi & Al-Amoudi (2018) acknowledge the role of these friction coefficients as candidates for model calibration, as Bassi et al. (2023) did regarding friction of movable beds.

As outlined in the Introduction, the objective of this paper is to propose a methodology that morphodynamical modellers can seize to get free from a time-consuming trial and error calibration of any unclosed parameter. This methodology will be described and applied to the 2P2L model, with the objective of proposing closure equations for the friction coefficients.

Numerical experiment

Many ML strategies rely on the use of large data sets from the literature. However, rather than using experimental data from the literature, limiting the training window of the CM to the conditions of these data, we propose to use the 2P2L model to directly generate data whose flow conditions fall within the intended scope of the model. Contrary to measured river flows for which governing quantities such as water depth and bed elevation are difficult to obtain or predict, these can be known exactly for the simple 1D numerical experiment we propose here.

The considered flow corresponds to a mixture of constant water and sediment discharges fed to a straight erodible slope with a knickpoint downstream, ensuring the complete transmission of the mixture beyond it and the rotation of the bed around it (point A, see Figure 2).
Figure 2

Sketch of the numerical experiment's configuration. Point A corresponds to the pivot point of the erodible material. and () correspond to the input (output) water and sediment discharges, respectively.

Figure 2

Sketch of the numerical experiment's configuration. Point A corresponds to the pivot point of the erodible material. and () correspond to the input (output) water and sediment discharges, respectively.

Close modal
The value of the sediment discharge corresponds to the transport capacity computed with whatever transport law from the literature (e.g. Meyer-Peter & Müller 1948) and considering the expected bed slope equal to the energy slope used to compute the value of the total shear stress involved in most transport laws:
formula
(4)
with being the hydraulic radius.

The goal of the experiment is to reach the same bed slope with the 2P2L model as that imposed to compute the solid discharge. The value of a given parameter to close is increased gradually to cover an appropriate range of values and the value that gives the closest slope to the expected one is retrieved as the optimal value for that situation. We refer to situations for every combination of hydraulic and morphological factor values.

Generation of situations

The situations used to develop the CM correspond to a set of 100 combinations of random values between two arbitrary limits for five independent factors: the median grain diameter – the 2P2L model using a single class of sediment – the relative density of the grains , the porosity of the bed , the bed slope , and the unit water discharge . Except for that corresponds to the objective result of the experiment, all other factors are directly provided to the 2P2L model as boundary conditions (i.e. ) or model input parameters. The chosen limits correspond to laboratory flows such as those discussed in this paper. These limits are listed in Table 1. The complete set of situations is provided in the Supplementary Material. This set was limited to 100 elements to limit the computational cost. Moreover, the situations with valuable results could be fewer than expected because no value led to a stable situation. Although many ML methods are often based on richer data sets, the ML method considered here is adapted to this limitation, as explained below. Other morphodynamical PBNMs may consider other factors that better suit their nature and other design plans according to their computational resources.

Table 1

Factor values used in the different situations of the numerical experiment

Factor (units)Limits
  
  
  
  
  
Factor (units)Limits
  
  
  
  
  

Also, this experiment is based on steady flows while the model aims to simulate transient flows. Neglected complex time-related processes might be a reason for a CM's performance lower than expected. Finally, the CM depends on a transport capacity formula, which was itself developed for a limited range of flows. If the training conditions of the CM fall out of this range, poor performance can be expected.

Data collection and treatment

The development of a CM for starts with the collection of the PBNM's simulation results. For every situation, a stable state (i.e. an equilibrium slope ) should ideally be met for at least one value of after some appropriate duration.

If several parameter values lead to stable slopes, the parameter value that gives the closest slope to the expected one is retrieved, as well as the difference between and (Figure 3). This difference can indeed be used to compute a weight w. A heavier w is attributed to the equilibrium slope that is the closest to , and that will lead to this situation becoming more influential on the regression model.
Figure 3

Example of two different situations S1 (blue) and S2 (red) with the same expected slope (black) leading to two optimal parameter values with their associated equilibrium slope . Situation S2 will have more weight on the CM than situation S1 because its equilibrium slope is closer to .

Figure 3

Example of two different situations S1 (blue) and S2 (red) with the same expected slope (black) leading to two optimal parameter values with their associated equilibrium slope . Situation S2 will have more weight on the CM than situation S1 because its equilibrium slope is closer to .

Close modal
Following the approach of a weighted linear regression, the difference between calculated and expected slopes is used to compute these weights:
formula
(5)
where the difference between slopes is divided by to treat results from situations with different equally.

Also, with the objective to later build a CM with a stepwise weighted linear regression method, the normality of the response variable must be ensured. A p-value in a Shapiro–Wilk test (Shapiro & Wilk 1965) can be considered satisfactory in this perspective. If the values cannot be considered as normally distributed, a transformation of the response variable, from to , can be performed or outliers removed, at the cost of information loss. The distribution of should not significantly deviate from normality.

Development of the closure model

Once the optimal values of were collected and treated for every situation, the situation set is divided into a training and a test set. The partition is made randomly, based on the values, so that the training and test sets have similar mean and standard deviations.

Then, any appropriate ML method can be applied to find a CM. Here, to preserve interpretability and ease of numerical implementation, we propose to use a stepwise weighted linear regression method. A stepwise method iteratively progresses and finds the CM that best satisfies any criterion provided. For instance, it can select the CM that gives the highest weighted :
formula
(6)
with corresponding to the number of elements in that set, to the observed optimal value of , to the computed value of according to the CM, and to the weighted average of the observed values y.

Also, because of the limited nature of the training set, a -fold CV is used. This technique consists in dividing the training set into k folds. The test set should never be used to train the model. Alternatively, each of the folds will be used as a validation set while the others will be used to train k models with any method, e.g. the stepwise method. The criterion used to select a model for each iteration (i.e. each fold being used as a validation set), is the maximization of the when applied to the validation set, , and considering the weights determined previously. Eventually, the model that gave the highest throughout all iterations is retrieved as the final CM. The number of folds k should split the data evenly and several repetitions of the CV can be performed to obtain an averaged model that can eventually become more generalizable.

Once a CM is determined, it can be applied to the test set, with consideration of the weights. Both and should be high enough to avoid underfitting. Moreover, they should ideally be rather similar to make sure that the training set was not overfitted. Playing with the number of folds k and the number of repetitions could lead to different CMs that can be compared with each other. This whole step can be achieved easily in JMP® Software (2023), a predictive analytics software with handy ML features.

Finally, the resulting HM can be applied to any other test case than the numerical experiment upon which it was built.

Resulting methodology applied to the 2P2L model

The 2P2L model described in (1) depends on two parameters and to calibrate. The methodology followed to develop CMs for these parameters is represented in Figure 4. Throughout the procedure, some situations might be excluded because the PBNM did not provide any stable solution.
Figure 4

Procedure followed to close models' parameters. The blue boxes correspond to steps performed with the help of JMP®. These operations could yet be achieved with any other method. White boxes correspond to manual operations while green boxes refer to 2P2L simulations. Red arrows indicate that some situations were excluded from the analysis for they did not deliver relevant results.

Figure 4

Procedure followed to close models' parameters. The blue boxes correspond to steps performed with the help of JMP®. These operations could yet be achieved with any other method. White boxes correspond to manual operations while green boxes refer to 2P2L simulations. Red arrows indicate that some situations were excluded from the analysis for they did not deliver relevant results.

Close modal
Figure 4 describes the application of the method to two dependent parameters. This is why there are prediction (Sim 1) and correction (Sim 3) steps. Closing a single parameter allows skipping the correction step (cf. bypass). This can be achieved here by assuming an existing CM of , for example based on the empirical relations found by Cheng & Zhao (2017) who studied the difference between static and dynamic repose angles of uniform sediment grains. Although repose angles are different from friction angles (Al-Hashemi & Al-Amoudi 2018), we assume a negligible difference since we are dealing with a prediction step. We then further associate the upper limit of the repose angle according to Cheng & Zhao (2017) with the static friction angle and the dynamic repose angle with the dynamic friction angle so that we have:
formula
(7)

Thanks to this preliminary CM, it is possible to provide the 2P2L model with couples that can be run through gradually. Since this formulation keeps simple enough, finds numerical counterparts in the 2P2L model for its parameters, and correlates the two friction coefficients, it stands out as a good candidate for initial estimations.

Two blocks of simulations (Sim 1 and 2) are hence dedicated to assessing the relevance of the CM of Cheng & Zhao (2017) for , and to improve it if not satisfying. In Sim 1, the values of are increased gradually from to , corresponding to repose angles from to . Coefficient is then computed with (7). Once the optimal couples were found, Sim 2 starts with set to its optimal value for that situation while varies from to , which corresponds to a repose angle of The optimal values coming from that second iteration are then compared with , i.e. those computed with (7). If and show excessive under- or overfitting, a new CM is developed for , based on . For generalizability, the explanatory variables of the regression model should be non-dimensional. Based on the aforementioned five independent factors and considering a steady flow, we use the following eight non-dimensional factors as possible explanatory variables: s, , , the Froude number , the ratio , the non-dimensional grain diameter where is the kinematic viscosity of water supposed to be constant, the Reynolds number where is the shear velocity, and the non-dimensional shear stress .

Sim 3 can then perform the same job as Sim 1, but is now computed with the CM found previously rather than with (7). The optimal provided by Sim 3 (or Sim 1 if and were not significantly different) can then be used to develop a new CM for . Eventually, the HM including these new CMs can be applied to real test cases.

First, the developed CMs for and are presented, considering the sediment transport formula of MPM for comparison purposes with the literature. The values of the parameters provided by the formulae of the CMs below should be bounded by the minimum and maximum values of and . Also, they should not be applied directly to other PBNMs for they were developed with data generated by this particular PBNM. Secondly, the resulting HM is applied to two laboratory test cases.

Closure model for μs

From the 100 situations, Sim 1 and 2 could not provide valuable results for only one situation. The eventual distribution of the resulting optimal has a p-value in the Shapiro–Wilk test. The normality of the distribution is hence not rejected. The 99 remaining situations were split randomly, based on , between a training set (77 situations) and a test set (22 situations).

Figure 5(a) shows that the empirical expression (7) is not generalizable to the whole dataset. While can be considered acceptable, shows a complete underfitting for the test set. Although the situations presenting the most significant weight, that are mostly included in the training set (large blue circles in Figure 5(a)), align well with (7), it looks to underestimate in many situations, especially when is low. Developing a new CM for hence appears necessary.
Figure 5

Comparison between (a) and and (b) and . The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to .

Figure 5

Comparison between (a) and and (b) and . The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to .

Close modal
A 7-fold CV with a single repetition was applied to the stepwise method to find the CM with the best on its associated validation fold. The following CM was found:
formula
(8)

Equation (8) led to limited over- and underfitting: and (Figure 5(b)). These results can be considered acceptable and the developed CM was used for Sim 3 to develop a CM for . Besides and that were already considered by (7), (8) also receives a significant contribution from h through the ratio . For river flows however, this factor can be expected to have less influence, for h is usually much greater than . Eventually, Figure 5 illustrates how the methodology can yield meaningful results, despite the simplicity of the initial CM. Moreover, an inappropriate choice concerning the initial CM can be overcome through the correction step of the proposed methodology.

Closure model for μd

From the remaining 99 situations that led to (8), 63 situations provided valuable results. The p-value in the Shapiro–Wilk test is satisfying () and no transformation of the response variable is hence necessary. Based on , this set was randomly split between training (45 situations) and testing (18 situations).

The same stepwise method was applied and used a 5-fold CV. It led to the development of a CM explaining the vast majority of the response's variability ( and , see Figure 6). The proposed CM is as follows:
formula
(9)
Figure 6

Comparison between and The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to

Figure 6

Comparison between and The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to

Close modal

Equation (9) shows a dependency on a wider set of factors – including their interactions – than (8) does. This illustrates well the difficulty of theoretically describing the physics behind friction with particles in movement.

Application of the hybrid model to laboratory test cases

Equations (8) and (9) were then implemented into the 2P2L model. The resulting HM was first applied to a 1D dam break on a bed of polyvinyl chloride (PVC) pellets. Spinewine & Zech (2007) provided the limits of the different layers identified in the flow while Spinewine & Capart (2013) delivered concentration results. This test case is hence used to evaluate the capacity of the HM to correctly predict the evolutions of the bed and water surfaces, the depth of the lower layer, and its concentration. Secondly, the performance of the HM is compared to other models from the literature on a benchmark test: a dam break occurring on a wide erodible floodplain made of sand (Soares-Frazão et al. 2012).

One-dimensional dam break on a bed of PVC pellets

The experiment referred to here was described by Spinewine & Zech (2007). They studied different dam break flows for six different configurations in a 6 m long and 0.25 m wide flume. These configurations always consisted of a 3 m long reservoir filled up to while the erodible bed material up- and downstream from the dam as well as the water level downstream of it could vary. Here, only configuration (a) is considered, i.e. a 0.1 m deep PVC bed spreading over the whole flume and no water downstream (Figure 7).
Figure 7

Initial configuration of the 1D dam break test case as described by Spinewine & Zech (2007).

Figure 7

Initial configuration of the 1D dam break test case as described by Spinewine & Zech (2007).

Close modal
The grain properties of the PVC pellets are as follows: , , and . Considering that is representative of the whole grain size distribution and following Strickler (1923), the Manning coefficient is computed as:
formula
(10)
which gives in this case.
Figure 8 reports different results obtained at by the 2P2L model with constant couples. Following Fowler & Wyatt (1960) who reported a difference ranging from to , six possible couples are displayed. These show the extent of the variability of the results and the large dependence of the PBNM on the friction coefficients. It appears that higher values of the friction coefficients tend to underestimate the erosion and the depth of the lower layer, and vice versa (Table 2). Larger differences also look to mishandle the inertia of the two phases for the wave front arrives too quickly when .
Table 2

RMSE (mm) around the different depths highlighted by the measurements when the 2P2L model is applied with different constant couples

hsTotal
    
     
     
     
     
     
hsTotal
    
     
     
     
     
     
Figure 8

Application of the 2P2L model with constant couples to the 1D dam break. Results obtained at .

Figure 8

Application of the 2P2L model with constant couples to the 1D dam break. Results obtained at .

Close modal
The results of the 2P2L model used with constant couples can be compared with those of the HM (Figure 9). The latter performs particularly well to estimate and h as well as the velocity of the wave front but underperforms when evaluating (Table 3). It actually provides results rather comparable to the couple , which makes sense as the wave front is dominated by values around and as can be seen from Figure 10. At the back of the wave front, and rise while shrinks. Interestingly, where velocity is non-zero, the range of is close to – as estimated by Fowler & Wyatt (1960), which advocates for the physical relevance of the developed CMs (8) and (9). Yet, the morphodynamic changes look mostly influenced by the wave front with as no striking difference can be noticed with the couple. Globally, the HM delivers an average performance when compared to the results obtained with constant couples. It hence provides valuable results that do not require the time-consuming trial and error step.
Table 3

RMSE between the numerical and experimental results for the 1D dam break

hsTotal
 7.92 11.68 14.14 33.75 
 10.32 16.37 17.46 44.14 
 5.54 16.78 11.14 33.45 
 6.58 17.91 10.93 35.42 
hsTotal
 7.92 11.68 14.14 33.75 
 10.32 16.37 17.46 44.14 
 5.54 16.78 11.14 33.45 
 6.58 17.91 10.93 35.42 
Figure 9

Comparison between the experimental and numerical estimations of the different interfaces of the 2P2L for the 1D dam break.

Figure 9

Comparison between the experimental and numerical estimations of the different interfaces of the 2P2L for the 1D dam break.

Close modal
Figure 10

Spatio-temporal evolution of the friction coefficients and their difference along the centreline of the flume.

Figure 10

Spatio-temporal evolution of the friction coefficients and their difference along the centreline of the flume.

Close modal

Figure 9 further illustrates the evolution of the flow, experimental results being available every quarter of a second. At , the hydrostatic and depth-averaging assumptions of the model fail. The analysis of the results is hence made from to . It appears that the front progression is correctly simulated by the HM, which indicates that the inertia of the two phases is sufficiently well captured. The position of the bed surface is rather accurately evaluated throughout the simulation while the estimation of h improves with time, as vertical and turbulent processes vanish. Yet, the error keeps superior to that around because of a local depression close to the gate in the numerical results that is not observed in reality. Finally, there is no surprise in observing a rather large error around since the 2P2L model is not meant to reproduce 3D processes like bedforms which lead to the oscillations of . Moreover, the confidence around its exact real position is certainly the lowest as it depends on chaotic processes. Nevertheless, it turns out the HM generally overestimates the level of , especially close to the wave front. Furthermore, it appears that this error rises with time. Yet, this error seems acceptable as the order of magnitude of the lower layer is correctly estimated by the HM.

The concentration measurements were published by Spinewine & Capart (2013). They are represented alongside numerical results at in Figure 11. These experimental results correspond to the aggregation of vertical concentration profiles from several experimental runs. Most of the time, two profiles from two different runs at the same position are represented as a range to illustrate the variability of the experiment. Figure 11 highlights the almost systematic underestimation of the concentration in the lower layer by the HM, despite the inertia being well captured. This can be partly explained by the fact that the HM tends to overestimate . The same conclusions can be drawn throughout the whole experiment.
Figure 11

Comparison between the numerical estimation of and the experimental concentration values for the 1D dam break. The experimental measurements correspond to one or two runs, depending on the profiles' locations. When two runs are concerned, the boundary values of concentration are plotted in black while the range between these boundaries is plotted in light grey. A single run is represented in black. The concentration values C are added to the position x of the vertical profile they correspond to, which is indicated by the vertical dashes. The coefficient 0.1 is arbitrary and chosen for readability reasons.

Figure 11

Comparison between the numerical estimation of and the experimental concentration values for the 1D dam break. The experimental measurements correspond to one or two runs, depending on the profiles' locations. When two runs are concerned, the boundary values of concentration are plotted in black while the range between these boundaries is plotted in light grey. A single run is represented in black. The concentration values C are added to the position x of the vertical profile they correspond to, which is indicated by the vertical dashes. The coefficient 0.1 is arbitrary and chosen for readability reasons.

Close modal

Dam break on a wide erodible floodplain made of sand

The experimental data of this second test case were provided by Soares-Frazão et al. (2012). They reproduced a dam break on a bed made of sand with a constant dam breach, narrower than the flume, up- and downstream from the dam (Figure 12). Upstream of the dam, the water level reached while the bed kept dry downstream of it. The 8.5 cm sand bed extended from 1.5 m ahead of the dam's gate up to 9 m downstream of it.
Figure 12

Description of the experimental flume of the 2D dam break, with the courtesy of Soares-Frazão et al. (2012).

Figure 12

Description of the experimental flume of the 2D dam break, with the courtesy of Soares-Frazão et al. (2012).

Close modal

The sand grain parameters are as follows: , , and . Equation (10) is used again to compute .

Soares-Frazão et al. (2012) provided the equilibrium topography of the bed downstream of the gate, after drainage of the water (Figure 13). An erosion pit can be found directly downstream of the gate while deposits form along the edges of the flume and then collide 7 m downstream of the gate, forming an arrow.
Figure 13

Plan view of the equilibrium topography of the 2D dam break as measured by Soares-Frazão et al. (2012).

Figure 13

Plan view of the equilibrium topography of the 2D dam break as measured by Soares-Frazão et al. (2012).

Close modal
Figure 14 shows two equilibrium topographies computed with the 2P2L model. In Figure 14(a), only constant friction coefficients and were used, as proposed by Di Cristo et al. (2016) for that test case run with the transport law of MPM. Figure 14(b) shows the results obtained by the HM, i.e. using the new CMs (8) and (9). Considering the whole topography, the root mean square error (RMSE) is a bit lower with constant friction coefficients (2P2L-), even though both simulations lead to the same order of magnitude. However, although the erosion pit and the collision between the lateral deposits are observed in both cases, none of them can correctly reproduce the equilibrium distance of the deposit front, as highlighted by Figure 15 ( and ). Yet, the HM predicts an equilibrium deposit front further from the gate and smoother than does 2P2L-. However, the amplitude of the deposits along the edges looks highly underestimated by the HM. This is clearly visible in Figure 15 ( where the HM barely reproduces the peak of deposit. Yet, no model comes really close to the real deposit's amplitude and they both appear smoother than reality. Unlike the case of 2P2L-, the depth of the erosion pit looks underestimated by the HM (Figure 15, ). Nevertheless, no model can reproduce the local features of the pit for they depend on 3D processes and are hence out of the scope of the considered models.
Figure 14

Plan views of the numerical equilibrium topography of the 2D dam break with (a) and and (b) the HM.

Figure 14

Plan views of the numerical equilibrium topography of the 2D dam break with (a) and and (b) the HM.

Close modal
Figure 15

Comparison between the experimental results and several numerical models around and along different longitudinal profiles for the 2D dam break. External results extracted from Di Cristo et al. (2016).

Figure 15

Comparison between the experimental results and several numerical models around and along different longitudinal profiles for the 2D dam break. External results extracted from Di Cristo et al. (2016).

Close modal

These two variations of the 2P2L model were compared to models from the literature also applied to this test case (Figure 15 and Table 4): a simpler mixture model (Canelas et al. 2013), a two-layer model (Swartenbroekx et al. 2013), and the 2P2L model of Di Cristo et al. (2016). Generally, the mixture model behaves the worst, both in terms of stability and RMSE. It appears that 2P2L- and the 2P2L model of Di Cristo et al. (2016) perform similarly in terms of RMSE, even though the 2P2L- model is more stable. A similar performance was expected since the theoretical model is the same and the transport and friction parameters are set equal. The HM and the two-layer model of Swartenbroekx et al. (2013) globally behave the best, except for in the case of the HM which underestimates the deposit as explained before. Altogether, the RMSE's order of magnitude is similar for the two variations of the 2P2L model and those of Di Cristo et al. (2016) and Swartenbroekx et al. (2013).

Table 4

Comparison of the RMSE around between several numerical models for three longitudinal profiles for the 2D dam break

RMSE (mm)y = 0.7 my = 1.45 m
2P2L-    
HM    
Canelas et al. (2013)     
Di Cristo et al. (2016)     
Swartenbroekx et al. (2013)     
RMSE (mm)y = 0.7 my = 1.45 m
2P2L-    
HM    
Canelas et al. (2013)     
Di Cristo et al. (2016)     
Swartenbroekx et al. (2013)     

The objectives of this paper were twofold. First, it presented a methodology aiming to circumvent the necessary and time-consuming trial and error calibration on past-obtained data of complex PBNMs containing unclosed and correlated parameters. Thanks to this novel procedure, these PBNMs can be run on any flume or river for which no calibration data exist. Secondly, a morphodynamical PBNM geared towards the simulation of transient flows highly laden with sediment, i.e. the 2P2L model, was presented. It served as a proof of concept for the proposed methodology and was applied to two different laboratory test cases.

The discussed methodology consists in closing a model for a parameter and is based on a numerical experiment involving the aggradation/degradation of an inclined slope of sediment (Figure 2). Depending on the parameters of the flow, a theoretical solution can be expected. For each flow situation, the optimal value that leads to the equilibrium slope the closest to the expected one is retrieved. It can then be used to develop a CM with ML based on many different non-dimensional flow parameters, each describing part of the flow's physics.

The PBNM described is a 2P2L model inspired by Di Cristo et al. (2016). This model treats separately the water and sediment inertia, which are related through drag (Figure 1). It also considers two layers of fluid with different sediment concentrations. It is based on eight governing equations (Equation (1)) which themselves rely on several parameters. CMs exist for most of them (see Di Cristo et al. 2016 and Supplementary Material, Appendix A) but static and dynamic friction coefficients and , involved in the stresses at the bed interface, lack convincing CMs from the literature. The methodology described above is then applied to these two correlated parameters to find appropriate closure formulations.

Based on an initial set of 100 different flow situations (Table 1), the initially assumed CM of (Equation (7)) proposed by Cheng & Zhao (2017) was rejected for it lacked generalizability (Figure 5). An alternative CM (Equation (8)), developed thanks to ML – a stepwise weighted linear regression method combined with a -fold CV – was proposed and emphasized the role of h in the resistance opposed by the bed to the flow. This new CM of was used to determine the optimal values of , following the same procedure. A CM of was then developed (Figure 6) and highlighted the multi-factorial dependence of that friction coefficient on flow parameters. These two CMs depended on the transport law that was chosen, i.e. MPM. Other transport laws could lead to different results. Moreover, even though these CMs were built with a generalizability objective, they are not expected to be used by other PBNMs. Besides the physics they try to mimic, they also compensate for the flaws of their associated numerical resolution. Nevertheless, the proposed methodology surely constitutes an interesting opportunity to get insight into the quantities involved in the physical processes of the concerned flows.

Finally, these CMs were included in the 2P2L to form a HM, which was applied to a 1D dam break on a PVC bed and a 2D dam break on a sand bed. The 1D dam break test case showed that the resulting HM was able to correctly capture the inertia of the flow's mixture and precisely simulated the progress of the wave front (Figure 9). It also proved to be accurate regarding the estimation of the bed level and the depth of the lower layer, despite a mild overestimation of the lower layer's depth. Consequently, the concentration was also underestimated most of the time though (Figure 11). Overall, it delivered a similar performance as the 2P2L model used with relevant constant couples (Figure 8). The application to the 2D dam break test case showed that the performance of the HM was also of a similar level as that of the 2P2L model whose friction coefficients would have been set by trial and error (Figure 14), as done by Di Cristo et al. (2016). Moreover, the HM performed similarly to other evenly complex models from the literature. Despite an excessive smoothness, it particularly stands out as the most stable model among the considered contenders. The HM approach is thus attractive for it eliminates the need for trial and error while ensuring satisfactory results.

Although the numerical experiment is meant for steady flows, the methodology presented in this paper proved to be effective when the 2P2L model was subject to transient flows. This methodology can actually be applied to any morphodynamical PBNM with a similar scope as the 2P2L. There must either be a single parameter to close or a group of correlated parameters such as and . A priori assumptions about parameters are hence necessary. Nevertheless, if computational results are limited, we advise not to go beyond two or three parameters to close. The number of simulations indeed rises quickly. If all assumed CMs are rejected, the number of simulations is , e.g. 10 simulations for , to be multiplied by the number of situations and values of . Other ML methods such as genetic programming or artificial neural networks (ANNs) could certainly provide similar or even better results, provided that the computational resources are sufficient (Goldstein et al. 2019). Yet, deep learning methods such as ANNs would certainly lead to a decrease in interpretability and ease of implementation in the PBNM.

The authors acknowledge the Fonds de la Recherche Scientifique (F.R.S.-FNRS) for the financial support of Robin Meurice as a Research Fellow. They would like also to thank Hervé Capart for the transmission of his experimental data and his kind support.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Al-Hashemi
H. M. B.
&
Al-Amoudi
O. S. B.
2018
A review on the angle of repose of granular materials
.
Powder Technology
330
,
397
417
.
Blau
P. J.
2001
The significance and use of the friction coefficient
.
Tribology International
34
(
9
),
585
591
.
Buffington
J. M.
,
Dietrich
W. E.
&
Kirchner
J. W.
1992
Friction angle measurements on a naturally formed gravel streambed: implications for critical boundary shear stress
.
Water Resources Research
28
(
2
),
411
425
.
Canelas
R.
,
Murillo
J.
&
Ferreira
R. M.
2013
Two-dimensional depth-averaged modelling of dam-break flows over mobile beds
.
Journal of Hydraulic Research
51
(
4
),
392
407
.
Cheng
N. S.
&
Zhao
K.
2017
Difference between static and dynamic angle of repose of uniform sediment grains
.
International Journal of Sediment Research
32
(
2
),
149
154
.
Combarros
M.
,
Feise
H. J.
,
Zetzener
H.
&
Kwade
A.
2014
Segregation of particulate solids: Experiments and DEM simulations
.
Particuology
12
,
25
32
.
Di Cristo
C.
,
Greco
M.
,
Iervolino
M.
,
Leopardi
A.
&
Vacca
A.
2016
Two-dimensional two-phase depth-integrated model for transients over mobile bed
.
Journal of Hydraulic Engineering
142
(
2
),
04015043
.
Fowler
R. T.
&
Wyatt
F. A.
1960
The effect of moisture content on the angle of repose of granular solids
.
Australian Journal for Chemical Engineers
1
,
5
8
.
Goldstein
E. B.
,
Coco
G.
&
Plant
N. G.
2019
A review of machine learning applications to coastal sediment transport and morphodynamics
.
Earth-Science Reviews
194
,
97
108
.
JMP®, Version 16. SAS Institute Inc., Cary, NC, 1989–2023
.
Martínez-Aranda
S.
,
Murillo
J.
&
García-Navarro
P.
2020
A robust two-dimensional model for highly sediment-laden unsteady flows of variable density over movable beds
.
Journal of Hydroinformatics
22
(
5
),
1138
1160
.
Meurice
R.
&
Soares-Frazão
S.
2020
A 2D HLL-based weakly coupled model for transient flows on mobile beds
.
Journal of Hydroinformatics
22
(
5
),
1351
1369
.
Meurice
R.
&
Soares-Frazão
S.
2022
A two-phase/two-layer model for fast and heavily loaded transient flows
. In
39th IAHR World Congress
,
20–24 June
,
Granada
.
Meyer-Peter
E.
&
Müller
R.
1948
Formulas for bed-load transport
. In
IAHSR 2nd Meeting, Stockholm, Appendix 2
,
Stockholm
.
Shapiro
S. S.
&
Wilk
M. B.
1965
An analysis of variance test for normality (complete samples)
.
Biometrika
52
(
3/4
),
591
611
.
Soares-Frazão
S.
&
Zech
Y.
2002
Dam break in channels with 90 bend
.
Journal of Hydraulic Engineering
128
(
11
),
956
968
.
Soares-Frazão
S.
&
Zech
Y.
2011
HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed
.
International Journal for Numerical Methods in Fluids
66
(
8
),
1019
1036
.
Soares-Frazão
S.
,
Canelas
R.
,
Cao
Z.
,
Cea
L.
,
Chaudhry
H. M.
,
Die Moran
A.
&
Zech
Y.
2012
Dam-break flows over mobile beds: Experiments and benchmark tests for numerical models
.
Journal of Hydraulic Research
50
(
4
),
364
375
.
Spinewine
B.
&
Capart
H.
2013
Intense bed-load due to a sudden dam-break
.
Journal of Fluid Mechanics
731
,
579
614
.
Spinewine
B.
&
Zech
Y.
2007
Small-scale laboratory dam-break waves on movable beds
.
Journal of Hydraulic Research
45
(
sup1
),
73
86
.
Strickler
A
.
1923
Beitrage zur Frage der Geschwindigheits-Formel und der Rauhigkeitszahlen fur Strome Kanale und Geschlossene Leitungen
, In:
Mitteilungen des Eidgenossischer Amtes fur Wasserwirtschaft
.
Bern, Switzerland
.
Swartenbroekx
C.
,
Zech
Y.
&
Soares-Frazão
S.
2013
Two-dimensional two-layer shallow water model for dam break flows with significant bed load transport
.
International Journal for Numerical Methods in Fluids
73
(
5
),
477
508
.
Vinuesa
R.
&
Brunton
S. L.
2022
Enhancing computational fluid dynamics with machine learning
.
Nature Computational Science
2
(
6
),
358
366
.
Zhu
L. T.
,
Chen
X. Z.
,
Ouyang
B.
,
Yan
W. C.
,
Lei
H.
,
Chen
Z.
&
Luo
Z. H.
2022
Review of machine learning for hydrodynamics, transport, and reactions in multiphase flows and reactors
.
Industrial & Engineering Chemistry Research
61
(
28
),
9901
9949
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Supplementary data