Comparing the Sensitivity of Bank Retreat to Changes in Biophysical Conditions between Two Contrasting River Reaches Using a Coupled Morphodynamic Model

Morphodynamic models of river meandering patterns and dynamics are based on the premise that the integration of biophysical processes matching those operating in natural rivers should result in a better fit with observations. Only a few morphodynamic models have been applied to natural rivers, typically along short reaches, and the relative importance of biophysical parameters remains largely unknown in these cases. Here, a series of numerical simulations were run using the hydrodynamic solver TELEMAC-2D, coupled to an advanced physics-based geotechnical module, to verify if sensitivity to key biophysical conditions differs substantially between two natural meandering reaches of different scale and geomorphological context. The model was calibrated against observed measurements of bank retreat for a 1.5 km semi-alluvial meandering reach incised into glacial till (Medway Creek, Ontario, Canada) and an 8.6 km long sinuous alluvial reach of the St. François River (Quebec, Canada). The two river reaches have contrasting bed and bank composition, and they differ in width by one order of magnitude. Calibration was performed to quantify and contrast the contribution of key geotechnical parameters, such as bank cohesion, to bank retreat. Results indicate that the sensitivity to key geotechnical parameters is dependent on the biophysical context and highly variable at the sub-reach scale. The homogeneous sand-bed St. François River is less sensitive to cohesion and friction angle than the more complex Medway Creek, flowing through glacial-till deposits. The latter highlights the limits of physics-based models for practical purposes, as the amount and spatial resolution of biophysical parameters required to improve the agreement between simulation results and observations may justify the use of a reduced complexity modelling approach.


Introduction
Despite increasing reliance on numerical modelling to simulate flow hydraulics, sediment transport, and bank erosion in rivers [1][2][3][4], several challenges remain when investigating natural channels.These are attributed to the complexity and spatial heterogeneity of processes related to soil properties, bank morphology, hydrology, and riparian vegetation [5,6].
Knowledge gains on lateral adjustments in natural channels have often emerged from studies undertaken at the scale of a single bank or river reach (e.g., [7,8]).The manner in which findings are presented in these studies may give the impression that alluvial river channels are affected in a similar way by external forces, independent of their scale and biophysical context.For instance, generalizations have been made such as soil cohesion increasing river bank stability [2,9] or vegetation stabilizing banks due to mechanical reinforcement [10][11][12].On the contrary, the surcharge imposed by mature trees on a river bank can have a destabilizing effect [13], in particular during the falling limb of a hydrograph [14].However, because each finding likely applies to limited river contexts (i.e., similar to the river from which they were drawn), the relative importance of biophysical variables may in fact differ considerably between river channels [15].Given the diversity of soil characteristics and heterogeneity of the floodplain with respect to biophysical conditions [16], there is an urgent need to develop tools that can be used with ease to evaluate the evolution of a diversity of alluvial and semi-alluvial river reaches [17].
In the last 10-15 years, several laboratory studies have examined the role of vegetation on increased cohesion in meander formation (e.g., [18][19][20]).Previous laboratory studies also used cohesive substrate, such as clay (e.g., [21][22][23]).However, most modelling (e.g., [4,24]) and flume-based studies (e.g., [25,26]) on meandering rivers involve sand beds and banks, whereas the meandering process typically occurs in cohesive floodplains [27].In particular, the mechanisms responsible for the development of meandering channels in a laboratory channel (e.g., [2]) may differ from those observed in nature, even if these channels share similar physiological traits.In addition, model validation is often achieved broadly with respect to visual elements captured using airborne imagery, such as channel planform and dynamics (e.g., [28,29]), with only a few modelling parameters representing the broad characteristics of soil and vegetation cover.For instance, bank material can be attributed to an erodibility coefficient [28], while plants are represented by a density value.The problem of how transferrable findings obtained at a given scale for a particular river type are to other channels remains scarcely documented.For example, are the key physical parameters contributing equally to bank stability in all river contexts?Very few studies have attempted to calibrate a model against data from a natural river channel/floodplain to help answer this question ( [30][31][32] being exceptions), and none, to the best of our knowledge, have compared parameter values between river channels.This could be attributed to weaknesses in the physics behind meander dynamics models, which do not take into account channel morphology (including bars) and vegetation [33].It may also be related to the substantial computational power required by these models [34] and to the scarcity or incompleteness of field datasets, owing to limitations in financial resources, time, available technologies, and data accessibility [35,36].This is particularly the case for bank retreat models [31].
River bank retreat involves a sequence of processes operating simultaneously at different timescales [37,38].For example, bend migration process is generally associated with periods that are much larger than that of the flow [39].Indeed, river banks are subject to slow transformations involving tension cracking [40], basal erosion by the flow [41], and riparian vegetation cover [42,43] and assemblage [13].Bank failures occur as soon as bank strength drops below a critical value.From a modelling point of view, this duality is challenging, and it has been addressed by at least two contrasting approaches.A linear framework-relying on near-bank excess velocity and an empirical representation of bank resistance (known as the HIPS formulation, from [44,45])-makes it possible to examine channel dynamics at long spatiotemporal scales (e.g., [46]).However, the consideration of basic flow and soil properties (velocity, cohesion), combined with the lack of a groundwater hydrology model and the impossibility to take into account complex channel and floodplain bathymetry [33,47], are serious limitations when analysing lateral retreat along natural river channels [16].A second popular approach is the use of enhanced computational fluid dynamics (CFD) models (e.g., [4,5,31]).The primary limitations of this approach are that it is computationally demanding and applicable to short spatiotemporal scales.They are also usually designed for curvilinear meshes (e.g., [48,49]), which often prevents their use with multithreaded channels [50].
This paper compares the sensitivity of river bank retreat to geotechnical parameter values between two natural river channels of different scale and geomorphological contexts in Quebec and Ontario (Canada).Both modelling investigations were undertaken with the hydrodynamic solver TELEMAC-2D [51] and sediment transport solver SISYPHE [52] from the suite open TELEMAC-MASCARET [53] v7.0.They were coupled to a physics-based geotechnical module that considers a broader set of soil properties than is commonly included in most bank erosion models [54].This research was based on three methodological novelties: (i) the use of a coupled CFD-geotechnical numerical model to examine the morphodynamics of a multithreaded river reach at the kilometre scale, (ii) the use of sedimentological and bathymetric data to calibrate models of river bank retreat, and (iii) the inclusion of groundwater hydrology into a coupled CFD-geotechnical numerical model and study.These lad to two main novel applications: (i) the identification of a set of biophysical conditions that fit observations of bank retreat for two different natural river channels, and (ii) a comparison of simulated bank retreat evolution between two natural river channels of different scales and geomorphological contexts.

Study Sites
Two contrasting fluvial environments with field observations on bank retreat were examined in this study (Figure 1a).The river reaches differ in terms of geological and geomorphological history, type of channel, spatial scale, and bank composition.The two sites were selected based on the availability of data against which to validate modelling outcomes.The first site is the downstream 8.6 km reach of the alluvial St. François River (SF hereafter), St. François-du-lac (Quebec, Canada), near its confluence with the St. Lawrence River (Figure 1b).A 3.4 km secondary branch forms an island within this reach.The primary branch has a sinuosity of 1.22.This channel has a bankfull (2-year return period) discharge of 1227 m 3 /s and a bankfull width of ~151 m when combining both branches.The longitudinal bed slope is around 0.04% in the main branch, and therefore markedly smaller than that of MC.The nature and lateral extent of riparian vegetation cover varied according to land use and to the occurrence of recent bank erosion episodes.A few banks were stabilized with riprap to protect roads near the island on both sides of the river.A narrow riparian zone-consisting of an assemblage of herbaceous plants, shrubs, and immature trees-is present next to agriculture lands.The trees were able to reach late succession stages in areas with less human impact near the deltaic area, near the downstream end of the island, and along the upstream left bank.A thorough analysis of the geotechnical properties of bank material revealed a sandy soil with several sandy silt, silty sand, and silt layers [55]; the silt material has a mean grain size of 0.075 mm.Mean grain size of the bed is 0.3 mm [56].The average recorded retreat rates in the two bends downstream of the confluence were 1.53 and 1.13 m/year between September 2008 and 2010; these values are lower than the average of 3.53 and 2.07 m/year between 1964 and 2008 [55].Channel bathymetry was surveyed from a boat using sonar and DGPS (differential global positional system) technologies [57].A total of 22,934 topographic points were acquired along 91 cross-sections, with an average distance of ~140 m between two consecutive cross-sections.Bed elevation values were interpolated linearly in the downstream direction.Floodplain topography was defined based on 194,631 LiDAR points acquired in 2001 and combined to the bathymetry dataset to create a digital elevation model (DEM).The SF river banks were also surveyed in the spring and fall of 2009 and 2010 to identify the nature of bank retreat.
The second site is a 1.5 km long reach of the semi-alluvial river Medway Creek (MC hereafter), London (Ontario, Canada) (Figure 1c).The study reach has a sinuosity of 2.31 and is located in a postglacial valley covered by different assemblages of deciduous and coniferous trees, shrubs, and herbaceous species.The lower bank soil layer consists of glacial till, which is buried under a thick sand layer covered with a thin organic layer.Bed substrate comprises gravel, sand, and till patches.The mean grain size of the unconsolidated material is 103.7 mm, although substantial variations existed between geomorphic units.Bankfull discharge, based on data collected at a gauging station just downstream of the study reach, is 43 m 3 /s, which corresponds to a channel width of ~20 m.Longitudinal bed slope is 0.27%.Further details on this site are available in [54].Over 5000 topographic points of the channel bed and banks were acquired in the study reach using a high-resolution DGPS instrument in November 2012.This manual technique is appropriate for the acquisition of morphological data in a channel with vegetated banks, as plants and large woody debris create visual obstructions.These points were combined with 1 m resolution LiDAR data [58] to create a DEM of the channel and floodplain.The field site was visited before and after each flood (>15 m 3 /s) between 2012 and 2015 to determine changes in morphology and vegetation cover along four sub-reaches of interest, herein referred to as A, B, C, and D (Figure 1c).The geometric properties of these sub-reaches are presented in (Table 1).Sites SF and MC, and associated sub-reaches, are shown in Figure 1.

Model Description
Physics-based geotechnical algorithms were coupled to the hydrodynamic module TELEMAC-2D and to the sediment transport module SISYPHE to include river bank retreat due to mass wasting.Only the basic features of the model are provided in this section.Further details are available in [59].
The hydrodynamic model TELEMAC-2D, which is described in [59], was set up with Smagorinsky's [60] formula to simulate turbulence while minimizing runtime.The default advection scheme (method of characteristics) was used to determine flow velocities and depth.The bed load was enabled using Meyer-Peter and Müller's [61] equation, included in the sediment transport module SISYPHE (described in [62]).The effects of local topography on transport magnitude and direction were taken into account using Koch and Flokstra's [63] equation.
The developed geotechnical module is powered by a genetic algorithm similar to the one proposed by [64] and by the equilibrium of forces described by [65] slope stability equation: where F s = safety factor; c = soil cohesion; L i = length of slice base i out of n; W i = weight of soil material; F cp,i = confining pressure exerted by the flow; δ i = angle between the result of hydrostatic confining force and normal to failure plane; U i = hydrostatic uplift force due to pore water pressure at slice base, basal angle β i , φ = friction angle of the soil material; and m i = m-term (Figure 2).A river bank hydrology model calculates water table elevation (z wt ) at time t using where t 0 = previous time, z fs = flow surface elevation at time t 0 , and k = rate of convergence of the water table elevation toward z fs .The constant k is adjusted according to the hydraulic conductivity of the bank material, and thus represents the rapidity by which water table adapts to a change in flow stage.Two k values are required: one for the rising flow stage, and one for the falling stage.

Simulated Flows
The largest flow discharge to be recorded at each site during the observation period was selected for subsequent numerical simulations.A verification was made to ensure that the selected flow discharges were indeed associated with bank retreat.Although flood hydrographs vary in magnitude, duration, and shape in nature, simple single-peak hydrological events were used in numerical simulations to limit the level of complexity, as morphodynamic modelling results are known to be affected by a large number of variables [66].The duration of simulations was substantially reduced, relative to the duration of natural events being simulated, to limit simulation times.
The numerical simulations analysed in this study did not use the flood hydrographs recorded at the examined field sites due to the extended period of time over which observations took place (2 and 3 years, respectively, for sites SF and MC).A flood with a peak discharge of 1928 m 3 /s, recorded on the St. François River at Hemming-Falls gauging station (ID 02OF002; Water Office, Environment Canada), occurred between 24 September and 15 October 2010.It corresponds to a return period of ~22 years.The single-peak hydrograph was compressed in time to become an 8 h simulation in which peak discharge is reached between t = 1 h and t = 4 h (Figure 3).The 3 h peak discharge is essential owing to the substantial amount of time required for a new hydrodynamic equilibrium to be reached following a change in imposed flow discharge (travel time from inlet to outlet is large due to the long channel combined with low average velocity) (Table 2).For Medway Creek, hydrograph shape is approximated using a γ function.Peak discharge is 60 m 3 /s-that is, similar to the maximum value recorded during the observation period (66.2 m 3 /s on 12 March 2013)-and event duration is 2.75 h (Figure 3).The selected discharge corresponds to a return period of ~2.5 years.Legend: t = simulation time (h), B = channel width (m), H = flow depth (m), V = flow velocity (m/s), Q = flow discharge (m 3 /s), Fr = V/(g H) 0.5 , Froude number, S 0 = longitudinal slope of free surface (%), and τ 0 = ρ g R S = bed shear stress (N/m 2 ) (ρ is mass density of water, g is gravitational acceleration and R is hydraulic radius).For site SF, H, and V are the average flow depth and velocity values for all points with H > 0.01 m, and B = Q/(H V).

Mesh Generation
The geotechnical module requires a large number of mesh nodes near river banks to allow for sediments eroded from the bank to deposit following a mass wasting event and correctly distinguish pre-and post-failure bank morphologies [54].Mesh resolution is slightly reduced in the areas that are unlikely to be affected by mass wasting (i.e., along channel centreline or away from river banks on the floodplain) (Table 3).Varying node density allows to minimize the total number of mesh nodes, and thus simulation time; this is essential to perform calibration on a complex model as it often requires running large numbers of numerical simulations.Therefore, a triangular mesh structure was built using the software BlueKenue [67].The software creates nodes using a dynamic moving front algorithm, connects them using unconstrained Delaunay triangulation, and smoothens the mesh with a Laplacian algorithm.For MC, elevation values from the DEM were assigned to a mesh structure that includes very small triangular elements (0.22 m 2 ) along unstable steep banks, small elements (0.87 m 2 ) along channel centreline, medium-size elements (2.16 m 2 ) in the riparian zone, and large ones (up to 10.83 m 2 ) away from the channel on the floodplain.The same strategy was adopted for SF, with the exception that a greater variety of element sizes was used to ensure smooth transitions between small elements along river banks (1 m 2 ) and larger ones near channel centre (25 m 2 ) and on floodplain (100 m 2 ).

Calibration and Boundary Conditions Flow
Calibrating the model required free surface elevation data at the inlet and outlet of the simulation domain.The method employed to measure free surface elevation values varied between field sites, as explained below.Following this, numerous hydraulic-only simulations were run with varying bed roughness values to adjust the energy slope so that it matched field measurements [68].
For SF, the calibration procedure was based on the analysis of high-resolution aerial photographs and hydrometric data.Flow width at the inlet was measured using an aerial photograph and was associated with the flow discharge recorded by the upstream Hemming Falls gauging station (ID 02OF002) on the day the photograph was taken (Q = 791.6 m 3 /s on 18 March 2016).A theoretical Manning's n bed roughness value of 0.0352 was calculated from flow conditions and cross-section geometry, which allowed for building a stage-discharge rating curve and imposing an unsteady flow discharge and level at the inlet.Near the outlet of the simulation domain, free surface elevation is markedly affected by the level of the St. Lawrence River.Historical hydrometric data (Q = 174.6 m 3 /s) from the Sorel gauging station (ID 02OJ022), located ~20 km upstream along the St. Lawrence River, combined with a second aerial photograph taken on 14 August 2009, was used to estimate free surface elevation (based on flow width) at the outlet of the simulation domain during the event that occurred on 24 September 2009.
For MC, flow measurements (depth and velocity) were taken along the inlet and outlet cross-sections of the study reach at low flow discharge (1.15 m 3 /s).A Manning bed roughness coefficient of n = 0.0153 produced an energy slope matching field conditions.This roughness value, combined with known cross-section morphology, were used to estimate water surface elevation at the outlet for larger discharge values.The anticipated difference in free surface elevation between low discharge (1.15 m 3 /s) and simulated peak flow discharge (60 m 3 /s) was 1.58 m.It compares well with the difference of 1.48 m recorded at a gauging station located 1.1 km downstream (ID 02GD008), considering the large longitudinal bed slope in this reach (4.17 m over 1.5 km).
The selected bed roughness values contrast with those obtained from a qualitative approach (e.g., [69]).For instance, Manning's n bed roughness coefficient at site SF is expected to be ~0.032.The model may require a slightly higher value due to the flow resistance created by the St. Lawrence River, located at the downstream end of the SF river.In the case of MC, the roughness value estimated using the methodology described in [69] is ~0.043.Determining bed roughness during calibration was based on the best adjustment with measured water levels.Differences between this numerical approach and qualitative estimations were therefore expected.

Geotechnical Analysis
The geotechnical calibration procedure comprised two steps and was applied to each one of the fluvial environments examined.The first step consisted of running a series of bank stability evaluations, assuming a straight bank profile (Figure 2c), to explore the model behaviour for a range of river bank dimensions (height and angle) and hydrological conditions (free surface and water table elevations) that can be encountered at the two field sites (row 'General' under 'Tree classification' in Table 4).The same procedure was followed with mass density (ρ) and friction angle (φ) parameter values set to known values for these sites (rows 'SF' and 'MC' under 'Tree classification' in Table 4).This calibration, which uses a simple bank geometry, served in the identification of key geotechnical parameters in the estimation of threshold parameter values leading to a safety factor near unity, and in creating a statistical model representing geotechnical processes for the two environments examined.The procedure was performed using machine learning algorithms-namely, random forest [70] and tree classification [71] (see Section 2.5).Here, random forest was employed to quantify variable importance.The binary tree-like structure resulting from tree classification is organized such that a splitting criterion (at any given node, based on the values associated with a variable with respect to the optimal threshold value) separates the data from its leaves and splits it into two sets, maximizing the increase in homogeneity from the node to its children.Here, the homogeneity of bank profiles with respect to safety factors was quantified using the Gini index.In addition to revealing divides and threshold values, a classification tree provides the frequency of occurrence of each leaf.In a second step, a three-dimensional calibration was performed using combinations of the two most influential parameters identified in the first step (i.e., soil cohesion, c, and friction angle, φ) (Table 5).Additional simulations were run while varying soil density ρ to verify if this could result in improved fit.Legend: c = cohesion (kPa), φ = friction angle ( • ) of bank material, and ρ = soil mass density (kg/m 3 ).Subscript correspond to identifiers of parameter sets for sites SF and MC.

Bank Retreat and Fitness
Transects placed at a regular interval, longitudinally along each reach (at each 100.0 m and 3.2 m, respectively, for sites SF and MC), were employed to classify bank locations in terms of stability; the categories employed were: stable, eroded by fluvial processes, eroded by geotechnical processes, and eroded by a combination of fluvial and geotechnical processes.For SF, the anthropogenically protected river banks that did not undergo any change were excluded from the analysis.
This study does not seek to replicate the bank retreat rates observed along sites MC and SF, but to identify the location of retreating river banks.The model is said to be successful, or to agree with an observation of river bank failure, if it predicts a retreat distance larger than a threshold value at a given transect.The selected threshold distances were set to 0.1 m for MC (0.5% of channel width) and to 5.0 m for SF (3.3% of channel width).River banks that have undergone a shorter retreat were considered stable.Note that the threshold value for SF corresponds to the actual average long-term retreat rate (see Section 2.1).For MC, in the absence of long-term historical data, the threshold value was selected because it maximizes the overall fit between observed and simulated retreat distances.
The fitness of a parameter set (all locations combined along the studied reach) with respect to bank retreat was quantified using confusion matrices and an alteration of Youden's J index [72]: where SN = sensitivity, SP = specificity, TP = number of true positive, TN = number of true negatives, FP = number of false positives, and FN = number of false negatives.Positive refers to the occurrence of a bank failure, whereas a negative refers to the lack of a failure.This rating method provided a constant range of indices, between 1.0 (correct prediction for every transect) and −1.0 (wrong prediction for every transect), which facilitated the comparison of fitness amongst sites and parameter sets.For the purpose of facilitating the interpretation of results, TP and TN were combined to indicate an agreement (or fit); the same was done with FP and FN to indicate a disagreement.
In the context of the present study, the overall similarity between observed and simulated retreat rates, all transects combined, is referred to as model fitness or agreement.The term accuracy could have been used if independent calibration and validation sets had been used as part of the modelling procedure.This is not the case here: the modelling experiment consisted of calibrating two models, one at each site, to determine the circumstances under which banks are unstable and describe differences in the influence of biophysical parameters between the models.The lack of a validation dataset implies that the models were not used to perform predictions.Therefore, the terminology employed reflects this situation.However, the metric employed here, Youden's J index, is compatible with both fitness or accuracy assessment.

Sensitivity Analysis Using Tree Classification
Tree classification was performed (within the software R [73] v4.1 and with the package rPart [71]) to build a visual representation of the combinations of parameter values (e.g., cohesion, compaction, heights of bank, friction angles) leading to similar safety factors, to identify key geotechnical variables, and to quantify the importance of each parameter.The set of rules behind a tree constitutes a statistical model.A simple straight slope with uniform bank material and spatial scales similar to the ones considered in this study was considered for this purpose (Figure 2c).A first tree was built by analysing river bank stability for the geotechnical and geometric properties encountered at the two field sites combined.Two additional trees were built (i.e., one for SF and one for MC) while imposing measured site-specific biophysical conditions.Here, each parameter set is narrowed down to the range of morphological and geotechnical conditions present at the associated field site.Therefore, the values of variables φ and ρ were set to those measured on site, and thus were kept constant.
The calibration strategy based on tree classification was selected to restrain the number of simulations to run with the coupled model.It is acknowledged that testing a larger number of combinations of parameters values would likely lead to a better fit with experimental data, but it would also be very time consuming in a context where each simulation takes several days to run (e.g., over 5 days for most the simulations presented here).In addition, this paper does not seek to accurately replicate the location and extent of bank failures, but rather to evaluate whether two contrasting fluvial environments are affected similarly by key biophysical conditions.This is why the model was calibrated without being validated.The calibration process, therefore, had a single purpose: demonstrate that parameters can be adjusted to fit observations.

Results
The presentation of results is done in three steps.In Section 3.1, the decision trees representing bank stability at both field sites are described, which shows the emergence of key parameters.
Note that these trees were built (see Section 2.5) by considering a large number of combinations of bank geometries and biophysical conditions (Figure 2c) and running the geotechnical stability model (without bank retreat) to evaluate the stability of each combination.In Section 3.2, decision trees are used to evaluate bank stability along transects, each of which is represented by a simplified geometry (Figure 2c).In Section 3.3, numerical simulations are run within the coupled hydraulic-sediment-geotechnical model (with bank retreat) for 12 combinations of parameter set values.

Identification of Key Geotechnical Parameters
The trees built using the rPart package algorithm were simplified to keep only one instance per sequence, along with its frequency.Looking at the general tree (SF and MC sites combined), the most frequent sequence is c-φ-α-h, which includes 33% of stability assessments in this sample (Figure 4).The other frequent sequences are c-φ-ρ (23%), c-φ-α (17%), and c-φ-α-ρ (16%).Note that c-φ is present in 94% of sequences, whereas h FS only influences the safety factor for 3% of the sample.Similarly, soil compaction and h WT are not present in any sequence.Overall, the most important variables are h, α, c, and φ (column 'General' in Table 6).Therefore, given that h and α are imposed by bank geometry at the field sites, and that h FS and h WT depend on flow conditions, c and φ are the geotechnical parameters that the model is most sensitive to.Site-specific classification trees (Figure 5) are subsets of the 'General' tree that is summarized in Figure 4.In addition, all variable sequences, with the exception of c-α-h-h FS , are present in both classification trees, and are arranged identically.Similarly, the relative importance of the variables is similar between sites (columns SF and MC in Table 6).However, the frequency of sequences differs slightly.For instance, there are fewer instances of the c-α-h sequence with site SF, which is compensated by a greater number of c-α-h FS and c-h-α.Therefore, h FS appears to have a greater control on bank stability at the SF site, and h is less influential, relative to the other parameters.
This calibration, which is based on parameter values related to geotechnical processes, suggests that bank retreat patterns along two contrasted river channels are expected to be similar.However, it is not clear whether this similarity will remain when taking into account sediment transport and complex channel bathymetries within morphodynamic simulations.This is tested in Section 3.3.But first, the achievable fit, based on tree classification, is examined with respect to each study site.

Evaluation of River Bank Stability Based on Tree Classification
The rules behind the condensed decision trees presented in Figure 4 were used to evaluate the stability of river banks along St. François River and Medway Creek within the studied reaches, based on a simplified representation of bank morphology (i.e., according to bank height (h) and angle (α) only) (Figure 2c).The calculations were performed on the observed geometries and bathymetries.This required to a priori extract h and α (Table 1) at each transect (Figure 6).The free surface elevation (h FS ) values employed were those calculated by a hydraulic model (sediment transport and geotechnical processes disabled), whereas water table elevation (h WT ) was assumed to be equal to h FS .The best fit (with a Youden's J = 23.3%) is obtained with c = 1.0 kPa and φ = 10 • at the SF site (Table 7).The fit drops to J = 19.0%with c = 0.5 kPa.The φ-value is compatible with the soil samples analysis performed by [55] (φ = 9.2 • ).For the MC site, the best agreement (J = 26.6%)occurs with c = 0.1 kPa and φ = 30 • , when considering the whole reach.However, other combinations of parameters c and φ also result in a good agreement between simulated and observed failures, for example, c = 1 kPa and φ = 20 • (J = 21.9%), or c = 2.5 kPa and φ = 10 • (J = 15.3%).The largest φ-value is compatible with the steep banks present at the field site and with the properties of highly cohesive soils.However, there are important differences in the level of agreement between sub-reaches A, B, C, and D (Table 7).The best agreement (J = 47.0%) is obtained in sub-reach B (with c = 0.1 kPa and φ = 30 • ).The worst agreement is found in reach A (Youden's J = 0% with any c-φ combination) and is barely acceptable in sub-reach C (J = 8.3%) with c = 0.5 kPa and φ = 10 • .In sub-reach D, the agreement is also poor, with the highest Youden's J (0.0%) having a high friction angle (φ = 40 • and c = 0.1 kPa) or high cohesion (c = 2.5 kPa and φ = 20 • ) (Figure 7).Note that the parameter ρ did not influence fitness, although it affected the safety factor values along a few analysed transects.
Although the overall fitness is fairly similar between SF (23.3%) and MC (26.6%), the range of combinations resulting in the best agreement with field observations is narrower for SF (0.5 ≤ c ≤ 1 kPa and φ = 10 • ) than it is for MC (a few combinations within 0.1 ≤ c ≤ 2.5 kPa and 10 ≤ φ ≤ 30 • kPa).A simple geometry (such as shown in Figure 2) is assumed with the dimensions provided in Table 1.The location of each sub-reach is shown in Figure 6.Legend: c = soil cohesion (kPa) and φ = friction angle (     5 and 9), respectively for sub-reaches A (transects 746-774), B (transects 775-827), C (transects 352-409), and D (transects 929-954).Numbers refer to transect identifiers.Numbers refer to transect identifiers.For SF, sub-reaches A and B include transects 1-88 and 89-123, respectively.

St. François River
The strongest fit with observations of lateral erosion along the banks of the SF channel (J = 50.2%) is obtained with parameter set SF 03 (Table 8), with a moderate friction angle (φ = 20 • ) and low cohesion (c = 0.125 kPa).This is similar to the soil texture at the field site (sand with φ = 9.2 • ).Most sets agree about the instability of the left meander bend downstream of the confluence (transects ~60-70 in box A L of Figure 8).Good agreement is observed, to a limited extent, with respect to the instability of the second meander bend (transects ~75-80 in box A R of Figure 8), despite a few false negatives at transects ~70-74.The agreement seems stronger along the internal and external river banks in acute bends (transects 29-42, 60-69, 75-83, and 107-123 in Figure 7a).However, the longitudinal extent of the unstable banks at these locations is greater than observed at the field site (e.g., see transects 57-58, 71-74, and 106).Finally, the most important discrepancy is found along the right bank near the bifurcation.5).White rows correspond to unmonitored transects.The subscript corresponds to the river side (L = left; R = right).

Medway Creek
The best overall fit, when combining sub-reaches A, B, C, and D, are associated with lower cohesion (c = 0.125 kPa) and high friction angle (35 9).This is compatible with the steep-angled banks of the semi-alluvial MC that are weakened by cracks and bioturbation.Two combinations of the parameters c, φ, and ρ also result in a good agreement between simulations and field observations (i.e., parameter sets MC 02 , MC 06 , and MC 11 ) (Table 5).The enhanced bank stability, owing to a high cohesion (from MC 06 to MC 11 ), is compensated by a decrease in φ to maintain the level of agreement.A similar observation can be made for set MC 09 (Table 9); the level of agreement is maintained by decreasing ρ (from 2100 to 1900 kg/m 3 ) to compensate for a reduction in c (from 0.500 to 0.375 kPa).The model is very sensitive to the value of all three variables and selecting the parameter value c is particularly critical; with φ = 40 • , any c-value ≥ 1 results in no mass failure.In addition, the agreement is stronger when φ = 35 • is selected in combination with 0.125 ≤ c ≤ 0.5 kPa (MC 11 ); selecting a lower or higher φ-value adversely affects the fit.Finally, conversely to the observations made for site SF, the model is sensitive to variations in ρ; overall, Youden's index drops when increasing ρ.
The model's performance with respect to any given parameter set varies substantially with the sub-reach considered (Table 9).Looking at the overall fit hides disagreements within sub-reaches A, C, and D. Indeed, the agreement is only good in sub-reach B (for all scenarios) and A (for MC 07 and MC 10 ).For instance, parameter set MC 05 (c = 0.25 kPa, φ = 35 • ) results in the best fit in sub-reach B, but is one of the worst parameter sets in sub-reach D. In the latter reach, the strongest agreements are obtained with c ≥ 1 kPa and φ = 40 • .Similarly, c ≥ 0.50, φ ≥ 40 • , and ρ ≥ 2000 kg/m 3 are required to maximize agreement in sub-reach C.This seems to indicate a spatial variation in geotechnical properties along the channel.This could be explained, at least partially, by the lack of consideration of the impact of vegetation and sedimentological layers in the model (see Discussion section).
Sub-reach A differs from the other sub-reaches due to much higher banks (20 m, relative to 2-4 m elsewhere).The fit varies between J = −62.5% and +25.0% in this sub-reach (Table 8).The best agreements are obtained with (0.125 ≤ c ≤ 0.375 kPa, φ = 35 • ) and (c = 0.50 kPa, φ = 40 • , ρ = 2000 kg/m 3 ).A slight increase in parameter value can substantially alter fit-for example, increase in φ by only 5 • between parameter sets MC 02 and MC 03 -which indicates that the model is very sensitive to geotechnical parameter values in this sub-reach.This variation could also be explained by the fact that the sum of TP and TN is very low (box A R in Figure 9).In contrast, the range of fitness values is lower in sub-reaches B and C. Legend: c = cohesion (kPa), φ = friction angle ( • ) of bank material, ρ = soil mass density (2100 kg/m 3 ), and J = Youden's index.The location of each sub-reach is shown in Figure 6.Maximum values are identified in bold.

Sensitivity
The SF model is more sensitive to changes in φ than in c (Table 8).The standard deviation (σ) of mean (µ) Youden's J is 0.9% for c values, compared to 16.1% for φ values.The first value is obtained by calculating µ at each row and calculating the standard deviation of the mean values.The second value is obtained by doing the same with columns.For site MC, σ is equal to 3.6% for c values, compared to 5.4% for φ values, when considering all sub-reaches of site MC.The same trend is observed for individual sub-reaches B, C, and D. However, sensitivity to φ is more important, with σ being equal to 11.2% for c values, compared to 23.7% for φ values.

Geotechnical Stability Model (without Bank Retreat) versus Coupled Model (with Bank Retreat)
A better agreement is obtained for SF with the coupled model (J = 50.2%;Table 8) than with the geotechnical stability model (J = 23.3%;Table 7).The optimal parameter values are also slightly different with φ = 10 • and c = 1.0 kPa for the geotechnical stability model (Table 7), relative to φ = 20 • and c = 0.125 kPa for the coupled model (Table 9).For MC, the best fits are similar with respect to Youden's index and optimal parameter values between the two model types.The indices are J = 26.6%without bank retreat (Table 7) and J = 22.5% with bank retreat (Table 8).The optimal parameters are φ = 30 • and c = 0.1 kPa for the model when strictly evaluating bank stability, compared to φ = 35 • and c = 0.25 kPa for the bank retreat model.However, the fit is relatively poor in sub-reach A with the geotechnical stability model for all values of φ and c (J = 0.0%; Table 7), which is less the case with the coupled model (maximum J = 25.0%;Table 8), although most indices are negative.However, the fit is stronger with the geotechnical stability model in sub-reach C (maximum J = 8.3% compared to 0.0%).5).The subscript corresponds to the river side (L = left; R = right).

Discussion
The two conceptual novelties of this study are the identification of a set of biophysical conditions that fit observations of bank retreat for two different natural river channels, and a comparison of simulated bank retreat between two natural river channels of different scales and geomorphological contexts.
The fact that the biophysical parameter values producing the best fit between observed and simulated bank retreat at site MC varies between sub-reaches suggests a variation in bank material composition over a relatively short distance.Indeed, each sub-reach at MC exhibits distinctive texture and layering (Figure 10).In sub-reach B, a sandy layer overlays glacial till, which also forms the channel bed.In sub-reach D, sediments are more thoroughly mixed vertically.Conversely, the analysis of two soil columns by [55] revealed a vertical profile consisting primarily of sand, but separated by a few silty sand, sandy silt, and silt layers.The greater uniformity of bank material in the latter case may explain that a stronger fit was obtained for SF (J = 50.2 in Table 9) than for MC (J = 22.5 in Table 8).Note that the presence of sedimentary strata was not considered in the model.The cohesion values that maximize model fit with river bank evolution at study sites are at least one order of magnitude lower than those commonly encountered in nature despite the fact that the model's ability to quantify stability was a priori tested against known problems of translational and rotational failures.For instance, [55] measured a cohesion of 13 kPa at the toe of a bank at SF (for a silty sample), which is ~100 times larger than the value associated with the parameter set that produced the best fit (Table 9).This seems to indicate that the model, even if it is physics-based, does not take into account a number of aspects of natural river banks that contribute to bank failures.The analysis of relatively homogeneous soil samples free of tension cracks in a laboratory may overestimate the overall strength of bank material found in nature.Finally, the fact that the observed lateral retreat corresponds to a timescale (2 and 3 years, respectively, for SF and MC) much larger than the simulation time (8 and 2.75 h, respectively, for SF and MC) may explain why a reduced cohesion value is required for the simulation outcome to fit observations.A marked difference in the cohesion value that maximizes fit was noted between MC sub-reaches (Tables 7 and 8).This may, at least partially, be attributed to the mechanical effects of the plant cover, which are lumped into the soil cohesion parameter.The model did not introduce spatial variations in soil characteristics; cohesion was identical at all nodes but varied between simulations.In natural rivers, the vertical variation in root density is such that the apparent cohesion is greater in the upper soil layer for several species [74,75], whereas the distribution of plant species and assemblage introduces a horizontal variation in soil strength, even assuming homogeneous soil material [13].The most unstable sub-reach at MC is B; it is densely vegetated but associated with the lowest cohesion value (Table 8).At this location, the removal of vegetation by floods, combined with the surcharge imposed by mature trees, seems to have triggered bank retreat.Conversely, sub-reach D is associated with the highest cohesion and friction angle (Table 8), perhaps owing to a vegetation stand consisting of herbaceous plants and young trees (Figure 10d), offering enhanced cohesion for little surcharge.However, an exposed bank section (upstream of the location shown in Figure 10d) was very stable during the observation period.Therefore, riparian vegetation seems to have a greater effect on bank stability in sub-reach B, even if it is less than 500 m away from sub-reach D. Similarly, spatial variations in root reinforcement may be responsible for the large number of false positives obtained in sub-reach A of site MC (transects 752-772).This steep bluff is topped by mature trees that may prevent its collapse.Failures have been observed in this area (Figure 10), but they are usually limited in extent to the area around a tree trunk falling off the bluff.For this sub-reach, the parameter set associated with the strongest agreement indicates that soil is moderately cohesive (c = 0.5) in this sub-reach, compared to sub-reaches B and D.
A better fit with observations is obtained by associating a set of parameter values to each river bank segment or floodplain patch with homogeneous biophysical conditions (soil characteristics, vegetation assemblage, etc.).For instance, the best overall fitness climbs from J = 22.5% to 41.8% for MC when considering the best parameter set at each sub-reach (A, B, C, and D).For SF, Youden's J only increases by 5.3% when varying parameter values spatially, which suggests more homogeneous soil characteristics at the SF field site, relative to MC.The suggestion to integrate floodplain heterogeneity into planform evolution models [26] would not only serve in improving the fit between observations and simulation results but would make the model more independent of input parameter values; the parameter set would remain valid throughout a simulation.However, care is required to avoid overfitting a model, which would otherwise compromise its capacity to be validated against a second dataset.Here, the possibility of overfitting comes from the fact that soil cohesion, despite the physics basis of the model, still represents the combined effect of at least two components, namely, soil and vegetation.This model characteristic would render the model less representative of the system with time if vegetation cover was to change.More work on complex river systems such as MC is needed to determine if, for practical purposes, a reduced complexity modelling approach would be more appropriate; it may not be realistic to attempt gathering biophysical data at the level of details required to obtain an acceptable fit.
Downstream conditions also differ between the sites.Both the inlet and outlet of the MC site are directly related to flow discharge.However, although the free surface at the inlet of SF is also related to the imposed flow discharge, the second half of the reach depends on the level of the St. Lawrence River in which the SF river drains.The selected free surface elevation at the outlet of this domain corresponds to the elevation recorded at a gauging station during the simulated event.During the observation period, a large number of combinations of inlet-outlet free surfaces have been encountered, but a single one of these combinations is examined in this study.We acknowledge that the differences in imposed hydrographs and boundary conditions between sites SF and MC, arising from the differences in scale and location of the sites compared, may affect results.However, free surface elevation (h FS ) and water table elevation (h WT ) were found to be less influential than geotechnical parameters (Table 6).It is thus unlikely that selecting a different event would have led to significantly different results.In addition, a good fit was obtained with respect to bank retreat for site SF.
One of the most important limitations of the coupled model is that it assumes that water table adjusts solely as a function of variations in the free surface of the flow.By doing so, it neglects the fact that banks may be fully saturated.This situation may explain the presence of false negatives between transects 380-395 at site MC.The area north of sub-reach C is partially submerged during the spring due to snowmelt, and a pond drains to the river by the remnants of a meander bend that was abandoned prior to 1942 (entering the floodplain at transects 377-384 and exiting at transects 403-410 and beyond) (Figures 6 and 10c).It is quite possible that the disagreement between simulated and observed failures in this area be attributed to the lack of a physics-based hydrological component that can be set up to consider realistic water table elevations across the floodplain.The presence of vegetation at the field site during failure events could have been detrimental when the soil was saturated, the mechanical effects of plants being outweighed by hydrological effects, as suggested by [16].Therefore, the fact that soil moisture content did not vary longitudinally along the bank according to preferential groundwater flows in the coupled model may have contributed to overestimating bank stability.
Calibrating a morphodynamic model can be time consuming, and detailed field data on bank retreat are seldom available.The integration of a geotechnical module, combined with large uncertainty regarding the value of geotechnical parameters to be used for a field site, renders the process more tedious as a larger number of trial simulations needs to be run to adjust model outcomes to observations.Here, a statistics-assisted calibration based on tree classification was completed to determine the range of parameter values leading to a safety factor near unity.Similar agreements were obtained with the geotechnical stability (without bank retreat) model and coupled model (with bank retreat) (see Section 3.4), which seem to indicate that initial channel bathymetry and biophysical conditions can be sufficient to estimate the location of bank retreat, without the need to consider hydrodynamics and sediment transport.However, caution is required, as these two model types generated slightly different parameter values (see Section 3.4).This difference could be attributed to the lag effect between water table and free surface elevations and to the consideration of complex bank morphologies by the coupled model.In contrast, water table elevation is equal to the free surface elevation, and bank profile is always straight with the geotechnical stability model.The coupled model also allows for one or several subsequent failures to occur in different portions of the slope.A potential application of machine learning algorithms would be to allow the geotechnical module to recognize and use the rules emerging from a decision tree directly into a coupled model.This would substantially reduce computation time.The module would remain partially physics-based, as the rules would have been pre-established using the geotechnical module.

Conclusions
This study sought to identify the most sensitive parameters in a morphodynamic model capable of taking into account mass wasting in a physics-based manner, and to verify whether the sensitivity to key geotechnical parameters differs between two contrasting fluvial environments.Our results indicate that lateral erosion is very sensitive to soil cohesion and friction angle and, to a lesser extent, to mass density of the bank material.A few combinations of these parameters resulted in a good agreement between simulation results and field observations, particularly in the alluvial St. François River case, where biophysical parameters of river banks and floodplain are more homogeneous.However, the agreement with field observations, and thus sensitivity, varies substantially from a parameter set to another, between sub-reaches, and between the two study sites.Sensitivity was also greater in some sub-reaches of the more complex semi-alluvial channel of Medway Creek.
The secondary objective was to devise a calibration method adapted to morphodynamic modelling that requires running as few simulations as possible.A pre-calibration phase, which used tree classification, and based on the assumption that the combined selected parameters values must bring the safety factor near unity, was used to estimate parameter values that are likely to result in an agreement with field observations.An iterative process was then used to run morphodynamic simulations within the coupled model, each time slightly varying the value of key geotechnical parameters in order to explore model behaviour and to improve fit with observations.The primary implication of the substantial degree of sensitivity found at both field sites is that morphodynamic models must account for spatial variations in geotechnical properties along a channel and must be reductionist enough to describe the complexity of the fluvial environment that they represent.For a highly complex semi-alluvial channel incised in glacial till, such as Medway Creek, it may be unrealistic to achieve this level of reductionism, and a reduced-complexity modelling approach may be considered.For less complex alluvial systems, such as the St. François River, the model employed in this study possesses these desirable characteristics, in particular, its physical base and ability to manipulate geotechnical variables to better represent bank material.In all cases, allowing to vary geotechnical properties across the floodplain should be the next step to improve morphodynamic models.

Figure 1 .
Figure 1.(A) Location map showing the two study sites.Numerical domain of (B) the St. François River (herein referred as SF) (Quebec; 72.908762 • W, 46.103081 • N), near its confluence with the St. Lawrence River and (C) a sinuous reach along Medway Creek (herein referred as MC) (Ontario; 81.289621 • W, 43.008107 • N).The white rectangle corresponds to the monitored MC sub-reaches.Arrows indicate flow direction.Elevation values are relative to mean sea level.Subscripts identify sub-reaches.

Figure 2 .
Figure 2. (A) Fragmentation of slump block into slices during bank stability assessment with Bishop's simplified method of slices.The dotted arc represents a hypothetic slip surface.(B) Stability was calculated by analysing the forces acting on slice base, namely pore-water pressure, confining pressure exerted by the flow, and soil weight.Variables are included in Equations (1)−(3).(C) Bank profile following deposition of the failure block at the friction angle of the bank material.In this study, this simplified bank profile also served to represent bank geometry during the comparison of model behaviour between the SF and MC reaches.

Figure 3 .
Figure 3. Hydrographs imposed to the coupled model during numerical simulations for (A) site SF and (B) site MC.

Figure 4 .
Figure 4. Simplification of tree classifications showing all decision paths leading to a similar safety factor.Each rectangle is associated with a decision along a path.Each red rectangle is a terminal node along a decision path.Percentages indicate the occurrence of the path within the dataset.

Figure 5 .
Figure 5. Simplification of tree classifications showing decision paths using fixed values of density (ρ) and friction angle (φ) based on field data for (A) SF sites, where ρ = 1950 kg/m 3 and φ = 9.2 • and (B) MC, where ρ = 2100 kg/m 3 and φ = 40 • .Each rectangle is associated with a decision along a path, e.g., c ≤ 2.5 kPa for the MC site.Each red rectangle is a terminal node along a decision path.Percentages indicate the occurrence of the path within the dataset.For instance, 76% of decision paths related to the MC site have the sequence c-α-h.Soil compaction was assumed to be equal to 0.75.

Figure 8 .
Figure 8. Agreement with observations for SF based on the components of the confusion matrix, i.e., FN = false negative, FP = false positive, TN = true negative, and TP = true positive.Rows indicate channel transects (see Figure 6), columns indicate parameter sets (see Table5).White rows correspond to unmonitored transects.The subscript corresponds to the river side (L = left; R = right).

Figure 9 .
Figure 9. Agreement with observations for MC based on the components of the confusion matrix, i.e., FN = false negative, FP = false positive, TN = true negative, and TP = true positive.Rows indicate channel transects (see Figure 6), columns indicate parameter sets (see Table5).The subscript corresponds to the river side (L = left; R = right).

Table 1 .
Channel and bank morphology.

Table 2 .
Hydraulics conditions for each site during the simulated flood hydrographs.

Table 4 .
Parameter values employed during calibration.

Table 5 .
Parameter sets and values.

Table 6 .
Importance of parameters (biophysical conditions and bank dimensions), expressed as Gini indices (or mean decreases in impurity).Impurity measures how well the trees included in a forest split a dataset.A large index indicates great importance.Note that the indices (importance) can be compared between variables, but not between the different samples (i.e., General, SF, and MC).
Legend: α = bank angle ( • ); h = bank height (m), h FS and h WT are proportions of h, respectively, for free surface and water table elevations; c = cohesion (kPa); ρ = soil mass density (kg/m 3 ); φ = friction angle ( • ) of bank material; k ASC and k DESC = water table adjustment rates with respect to free surface elevation, respectively, for the rising and falling limbs of the hydrographs.

Table 8 .
Fitness (Youden (1950)'s J indices) of the calibrated coupled model (with bank retreat) for site SF (maximum value identified in bold).

Table 9 .
Fitness (Youden (1950)'s J indices) of the calibrated coupled model (with bank retreat) for site MC.