https://orcid.org/0000-0002-3962-9602
ABSTRACT
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.
HIGHLIGHTS
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.
ABBREVIATIONS
INTRODUCTION
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.
METHODS
Two-phase/two-layer model
Schematic representation of the 2P2L model and its associated quantities.
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.
and 
(Equations 1(d)–(1e)) since it is a 2D model:
, 
, 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
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.
can be modelled in a Mohr–Coulomb way as:
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.
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.
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.
equal to the energy slope 
used to compute the value of the total shear stress 
involved in most transport laws:
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.
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.
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. 
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 
.
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 .
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.
:
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
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. 
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.
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.
, 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:
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.
RESULTS AND DISCUSSION
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).
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. 
Comparison between (a) 
and 
and (b)
and 
. The wider the circle, the heavier the situation. The smallest circle corresponds to 
and the largest to 
.
Comparison between (a) and and (b) and . The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to .
on its associated validation fold. The following CM was found:
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).
and 
, see Figure 6). The proposed CM is as follows:
Comparison between 
and 
The wider the circle, the heavier the situation. The smallest circle corresponds to 
and the largest to 
Comparison between and The wider the circle, the heavier the situation. The smallest circle corresponds to and the largest to
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
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). 
Initial configuration of the 1D dam break test case as described by Spinewine & Zech (2007).
Initial configuration of the 1D dam break test case as described by Spinewine & Zech (2007).
, 
, and 
. Considering that 
is representative of the whole grain size distribution and following Strickler (1923), the Manning coefficient is computed as:
in this case.
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 
. RMSE (mm) around the different depths highlighted by the measurements when the 2P2L model is applied with different constant 
couples
 | hs |  |  | Total |
|---|---|---|---|---|
 |  7 |  |  |  |
 |  |  |  |  |
 |  |  |  |  |
 |  |  |  |  |
 |  |  |  |  |
 |  |  |  |  |
| hs | | | Total |
|---|---|---|---|---|
| | 7 | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
| | | | | |
Application of the 2P2L model with constant 
couples to the 1D dam break. Results obtained at 
.
Application of the 2P2L model with constant couples to the 1D dam break. Results obtained at .
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. RMSE between the numerical and experimental results for the 1D dam break
 | hs |  |  | Total |
|---|---|---|---|---|
 | 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 |
| hs | | | Total |
|---|---|---|---|---|
| | 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 |
Comparison between the experimental and numerical estimations of the different interfaces of the 2P2L for the 1D dam break.
Comparison between the experimental and numerical estimations of the different interfaces of the 2P2L for the 1D dam break.
Spatio-temporal evolution of the friction coefficients and their difference along the centreline of the flume.
Spatio-temporal evolution of the friction coefficients and their difference along the centreline of the flume.
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.
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. 
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.
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.
Dam break on a wide erodible floodplain made of sand
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. 
Description of the experimental flume of the 2D dam break, with the courtesy of Soares-Frazão et al. (2012).
Description of the experimental flume of the 2D dam break, with the courtesy of Soares-Frazão et al. (2012).
The sand grain parameters are as follows: 
, 
, and 
. Equation (10) is used again to compute 
.
Plan view of the equilibrium topography of the 2D dam break as measured by Soares-Frazão et al. (2012).
Plan view of the equilibrium topography of the 2D dam break as measured by Soares-Frazão et al. (2012).
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. 
Plan views of the numerical equilibrium topography of the 2D dam break with (a) 
and 
and (b) the HM.
Plan views of the numerical equilibrium topography of the 2D dam break with (a) and and (b) the HM.
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).
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).
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).
Comparison of the RMSE around 
between several numerical models for three longitudinal profiles for the 2D dam break
| RMSE (mm) |  | y = 0.7 m | y = 1.45 m |
|---|---|---|---|
2P2L- |  |  |  |
| HM |  |  |  |
| Canelas et al. (2013) |  |  |  |
| Di Cristo et al. (2016) |  |  |  |
| Swartenbroekx et al. (2013) |  |  |  |
| RMSE (mm) | | y = 0.7 m | y = 1.45 m |
|---|---|---|---|
| 2P2L- | | | |
| HM | | | |
| Canelas et al. (2013) | | | |
| Di Cristo et al. (2016) | | | |
| Swartenbroekx et al. (2013) | | | |
CONCLUSIONS
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.
ACKNOWLEDGEMENTS
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.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=7142)