## 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

*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.

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.

*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

*et al.*(2016), the friction stress exerted by the solid phase on the bed may be expressed as: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 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 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.

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.

*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.

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.

*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

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}

_{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).

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}

_{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).

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

*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.

. | h
. _{s} | . | . | 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 |

. | h
. _{s} | . | . | 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 |

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.

#### Dam break on a wide erodible floodplain made of sand

*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.

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

*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.

*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.

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).

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.

^{®}, Version 16. SAS Institute Inc., Cary, NC, 1989–2023