A comprehensive method to estimate flood levels of rivers subject to ice jams: A case study of the Chaudière River, Québec, Canada

Flood delineation for cities and urban development is a common task to assess risk management. Unfortunately, flood hazard analysis is mostly done for open-water cases, and very little is done for ice-affected flooding. As highlighted by Kovachis et al. (2017), to limit mapping to only open-water events is to ignore a type of riverine flooding that has severe consequences in Canada. Only a few examples take into consideration both open-water and ice-induced annual exceedance probability (AEP) levels. Such an example is given by Tuthill et al. (1996), where the combination of peak stages under ice jam and no ice jam is analysed to produce a combined probability distribution on the Missisquoi River in Vermont. However, there are certain deficiencies due to the lack of historical data that can be collected regarding both ice jam and open-water events. Some towns and municipalities have complete records of each known event, while others have no such information. This logging can include ice jam length, location, gauge station, flow discharge, ice thickness (IT), water extent, etc. Therefore, specific measurements for the delineation of ice-related floods seem lacking (Lindenschmidt et al. 2018) in certain towns. Hence, these methodologies need to be extended to include a more encompassing flood analysis incorporating both ice jam and open-water AEP to have a broader view of what certain rivers can be entitled to in flood management (Tuthill et al. 1996; Tucotte et al. 2017). Examples of case studies which utilize the 1:100 ice-induced AEP have been worked on in past years, such as Tuthill et al. (1996), and more specifically, Lindenschmidt et al. (2016, 2021), Lindenschmidt (2023) and Rokaya et al. (2019) have used a stochastic modelling framework using RIVICE to produce flood maps, ice jam water level forecasting and jam mitigation analysis using the Monte–Carlo framework (MOCA). Ice jam analysis using HEC-RAS has been undertaken by Ahopelto et al. (2015) and Altonen & Huokuna (2017) using Visual Basic to control HEC-RAS to produce a MOCA framework of ice jam scenarios based on the normal distribution of ice jam parameters to map ice-induced floods. Other examples that incorporate both open-water and ice-induced flood risk analysis are described by Zufelt (2005) and Zufelt & Lefers (2019) and which utilize the FEMA (2018) probability equation for rivers located in the states of New York and North Dakota.

As part of a governmental flood hazard programme, assessing ice jam and open-water flood delineation of the Chaudière River, located south of Quebec City, Canada, has been done with the simulation model HEC-RAS with an external script using a MOCA framework written in Python to generate a distribution of ice jam flood events along the study reach while open-water flood delineation has been done by unsteady flow simulation over a 90 km reach of the river.

The objectives of this paper are to show the applicability of historical data in determining the 1:100 AEP ice-induced water levels with a MOCA framework; to assess the sensitivity of mechanical parameters and historical data, i.e., ice jam toe location, ice jam length, flow and IT, on simulated water levels and to determine which geomorphological features cause higher ice-flood levels than those for the open-water 1:100 AEP in this study reach.

Among these watercourses, two are major tributaries and 16 are minor tributaries that are selected to carry out the hydraulic modelling to include their lateral water supply which can be significant, given the surface area of their watershed (Table 1). Table 1 presents all the tributaries (ID) and the respective watershed area, as well as the watershed area of the Chaudière River at the confluence of each tributary.

The two major tributaries that contribute to the total volume of ice in the Chaudière River during a breakup are the Famine River and the Bras-Victor River. The Famine River (10) joins the Chaudière from the east at downtown Saint-Georges-de-Beauce (SGB) with a drainage area of 714 km². At the confluence, the drainage area of the Chaudière increases to 3,820 km². Approximately 7 km upstream from its confluence, a gauging station was built in 1964 and transmits data every 15 min. The Famine River can emit its ice into the Chaudière River since there is no obstruction and can therefore contribute to the volume of ice on the Chaudière River. Further downstream, the Bras-Victor River (8) reaches the Chaudière River between the municipalities of Beauceville (BCE) and Saint-Joseph-de-Beauce (SJB). Its watershed is situated west of the Chaudière River and covers an area of 735 km². The ice volume along the Bras-Victor River can also contribute significantly to the volume of ice on the Chaudière River.

The sixteen minor rivers and creeks selected for this study are located between the towns of SGB and Scott-Junction (SJ). The area of these tributaries varies between 11.86 and 216.4 km².

The analysis of the 1:100 AEP open-water levels was done by Bessar (2021). This analysis was first done by the calibration of a 1D model in unsteady flow with HEC-RAS using the April 2019 event, which covered flow rates from 500 to 2,100 m³/s. By achieving the water mass budget between the gauging stations of Sartigan Dam and Saint-Lambert, the flow distribution along the model was set up to have 60% of the total flow starting at Sartigan dam while having a uniformly distributed lateral flow equal to the remaining 40%. A more detailed description in Bessar et al. (2020, 2021a, 2021b) is available. Upon calibrating the 1D unsteady model, the open-water 1:100 AEP water levels were determined using a regional flood frequency analysis based on instantaneous peak data as described by Ricard et al. (2021) for the Chaudière River. The 1:100 AEP open-water levels determined along this reach are those that will be compared to the 1:100 AEP ice-induced water levels.

The HEC-RAS 1D ice model of the study reach took into account the properties of the bridges, the roughness coefficients for the channel, the banks and the floodplains. By integrating feasible roughness coefficients, the model can reproduce the water levels observed along the river for a given flow range.

To calibrate the model in open-water to have a full-bank flow, such as it occurs during an ice jam, the flow rate would need to range from 700 to 1,000 m³/s. The calibration of the hydraulic model was carried out in a steady state because HEC-RAS can only simulate ice jams under steady flow conditions.

In addition to the choice of the flow range, the time series for flow must precede that of the breakup period to take into account the Mannings' roughness of the flood plains which may be in the presence of snow and/or ice.

Table 2 shows the root mean squared error (RMSE) between the simulated and observed water levels along the reach. A value of 0.086 m was obtained for the RMSE. The roughness values are based on those provided in the HEC-RAS manual and those proposed by Chow (1959). The roughness of the riverbed is related to the size of the granular materials present, i.e., the larger the granular materials, the greater the resistance. In addition, the sinuosity of the river brings additional resistance, the more sinuous, the greater will be the resistance to the flow. Also, sections with a greater slope have a greater roughness because the shear induced by the flow velocities is the greatest and the sediments that resist the flow are the largest. Areas where there were discrepancies between simulated water surface (SWS) and the observed water surface (OWS) were calibrated until desired agreement. Roughness values ranged between 0.025 and 0.035 for simple geometry xsections and reached up to 0.068 for reaches that had rapids or boulders.

The model was calibrated for ice cover conditions using data available from the watershed management committee of the Chaudière River (COBARIC) and gauging stations from the Quebec Hydrology Expertise Centre (CEHQ) stations and data obtained during the 2019 and 2020 sampling campaigns.

A drilling campaign on the river totalling 62 sections was carried out in 2019 and 114 sections in 2020. IT and water levels were measured at each sampling site, including six along each of the 62 cross sections geolocated every 1,000 m between the municipalities of SJ and SGB and 10 for each of the 114 sections geolocated every 500 m between SJ and BCE. The ice thicknesses and the water levels were recorded concurrently with the bathymetric survey campaign in the winter of 2020 to complete the profile of the ice cover along the river in the study area between SJ and BCE.

Moreover, considering that the thickness of the ice cover varies considerably at the same section, for example from 0.18 to 0.99 m downstream of the Devils rapids and from 0.60 to 1.59 m at Île Ronde (BCE), the average obtained from each section was assigned to the corresponding section of the model. In the situation where there is no reading near a section of the model, a linear interpolation of the IT was carried out from the averages of each section.

The roughness of the ice was determined from the observations made according to Nezhikovskiy (1964). The three types of ice cover, namely sheets of ice, frozen frazil and unconsolidated frazil are the ones that are referenced. Frozen frazil is the type of ice that best represents the ice cover present on the Chaudière River at the start of freeze-up or following a winter breakup. A relationship between the ice cover thickness and the roughness of the ice cover is established with these observations.

From these results, we also note that the flow changes according to the watershed area can vary about the observed water levels, particularly at 62.5 and 61.7 km located near the Nadeau, Lessard and Belair Rivers and at 75.1 km located near the confluence of the Bras-Victor River. Because there are no flow gauging stations at these locations, flow estimation is difficult, especially at low flow. Rather than increasing the roughness to correct the simulated water levels, adjustments were made to the lateral flow by up to 20% to obtain water levels close to those observed at 84.7 and 75.1 km.

A second simulation was done for the 2019 ice cover conditions to confirm the plausible causes for differences in water levels. The sections chosen to compare the results are those which coincided with the sections of the bathymetric survey carried out in 2020. As for the flow rate, the mean value for the days at which the IT measurements were taken was 19.33 m³/s. Table 4 presents the 2019 ice cover results.

The RMSE obtained for the 2019 ice cover simulation was 0.37 m, which is above the tolerance of ±0.2 m as mentioned by Beltaos et al. (2012). At 87.3 km, the simulated water level is 0.95 m lower than the observed data and is most likely related to the in situ ice cover variation (0.34–0.81 m) which cannot be implemented into HEC-RAS, while combining the effect of an island downstream. For sections 61.7 and 48.3, the presence of frazil which was observed but difficult to measure was not considered in the simulated ice cover may have influenced the results.

The calibration parameters made it possible to highlight that the estimation of the flow in the winter period proves to be a crucial element for obtaining accurate values and that the lateral flows must be taken into consideration. Also, IT variation as seen in the simulation of the 2019 ice cover influenced simulated water levels when considered uniform. Therefore, these differences in elevation are caused by the limitations of the software, which include simulating the obstruction of the channel according to the real variations in IT at a section, the presence of frazil, the consolidation of ice sheets, the effect of the weight of the snow on the ice cover and the effect of cantilevered ice cover between rocks, particularly at low levels.

The HEC-RAS ice module is largely based on the principles of ice mechanics and ice jam equilibrium by Pariset et al. (1966), Uzuner & Kennedy (1976) and Flato & Gerard (1986). All of the equations used in connection with ice jams are presented in the HEC-RAS reference manual. The modelling of an ice jam using HEC-RAS is based on the following specific mechanical and physical parameters, which include IT, Mannings' n ice roughness, friction angle of ice rubble, ice porosity of rubble, k1 stress ratio, maximum water velocity for ice jam thickness iteration and ice cohesion.

For each ice jam simulation, the parameters were assigned to a minimum of two consecutive sections, considering that it is the minimal number of sections needed to simulate an ice jam and the total length of the jam depends on the number of cross sections to which the ice parameters are assigned.

In all cases, the role and influence of the range of values for each parameter chosen will influence the force balance equation and simulated water levels. Therefore, the global sensitivity analysis (GSA) will help discern which parameters have the greatest influence on the simulated back-water levels.

The location of ice jams in the study area was carried out using only the data listed since the commissioning of the ice control structure in 1967 to consider the new dynamics of the ice. Indeed, since the function of the structure is to retain ice upstream, it can no longer transit downstream of the ice control structure (ICS) as before and contribute to the total volume of ice available to generate ice jams. From 1967 to 2020, a total of 129 ice jams are listed. Some of these events are incomplete as to their known location, length or date. Among the 129 events, 87 have a known location which can be associated with a specific section of the hydraulic model.

The range of values for each mechanical and physical parameter required to simulate floods in the presence of ice using HEC-RAS was assigned based on historical data and the literature.

The IT is calculated based on cumulative degree days of freezing (CDDF) for the ice in the channel. CDDF are randomly determined based on the distribution of historical data. CDDF is related to IT using Stefan's equation. The Manning roughness ‘n’ of the ice cover is calculated using the equation adapted to the results of Nezhikovskiy (1964) to consider the initial IT. The ranges of values for the angle of friction (φ) and porosity (p) follow the range of recommended values (White 1999) which are, respectively, 40–60° and 0.35–0.8.

As for the lateral stress coefficient (K1), the HEC-RAS user manual does not offer any guidelines for determining the value of K1 (Beltaos 2018). Several authors, Beltaos (2018, 2019), Tuthill & Mamone (1998), Beltaos et al. (2012), and Rokaya & Lindenschmidt (2020) use 0.33 as the default value of K1 with HEC-RAS, as suggested by Flato & Gerard (1986). In the present study, a range of random values of 0.3–0.36 was chosen to consider the variability of the lateral stress coefficient.

Since the internal resistance coefficient (μ) is related to the angle of friction (φ), the porosity (p) and the lateral stress coefficient (K1), the combination of the three parameters respects the range of values recommended by White (1999), i.e., 1.0 ≤ μ ≤ 2.0. By iterative calculation, the values of the three parameters are chosen within their respective range to then calculate the internal resistance coefficient. In case the value of μ exceeds the limits of the range, the algorithm iterates new values for φ, p and K1 and again proceeds to calculate μ.

In the literature, Tuthill et al. (1998) mention a range of values for V_max between 1.21 and 2.43 m/s. In HEC-RAS, the default value of V_max is 1.524 m/s. In the context of this study, the range of values obtained following the reconstruction of the four historical events from 1.5 to 2.0 m/s. Beltaos & Tang (2013) recommend the use of a large value of V_max, which makes it possible to generate an acceptable ice jam profile since HEC-RAS does not simulate the flow through the interstices of the ice jam. More precisely, Beltaos & Burrell (2015) mention that the use of V_max = 10 m/s makes it possible to solve the local instabilities produced by the simulations in HEC-RAS with the default values. Studies carried out by Beltaos (1993, 1999) have shown that the flow under the jam occurs mainly through the interstices of the jam and that the flow velocity under the foot of the jam is relatively low. Therefore, the application of a large value of V_max under the foot of the ice jam compensates for the flow through the interstices. Therefore, the choice of the range retained for the values of V_max is between 1.524 and 10 m/s. Regarding cohesion, the default value used in HEC-RAS is 0 since in situ measurements cannot be performed. Moreover, since the pieces of ice are not consolidated, they are considered aggregates (White 1999).

The hydrotechnical method combined with the MOCA approach is based on the simulation of a multitude of scenarios to obtain different water level profiles. The generated scenarios were made using a script written in Python.

The simulation of the different scenarios begins with the recovery and integration of the input data which is divided into two groups, namely the processed historical data and the basic files of the calibrated hydraulic model which include geometry, flow, plan and project. The following process is an algorithm developed that allows historical data to be managed according to the goodness-of-fit, which is based on the Akaike information Criterion (AIC) (Akaike 1973). The processed data are presented in the form of histograms, including the mean and standard deviation.

In the fifth process, the HEC-RAS Controller in Python was called upon for each simulation by iteration and exports the raw values related to maximum water level, ice volume and flow at each river section of the model. The maximum number of simulations that can be made per batch file is N = 99. When the ‘N’ number of simulations is completed, it calls on an algorithm which retrieves the simulated water levels, ice volume and flow for each section. The ice volume and flow values are used to identify if HEC-RAS managed to produce an ice jam. In the initial steps of writing the algorithm, some scenarios did not simulate an ice jam or did not have an increase in flow along its reach. This leads to the sixth process. In the event no ice jam or flow change was simulated, the scenario was discarded from the results, as well as scenarios that produced an ice jam with a flow that was below or above a specific threshold. Scenarios that were kept could now be analysed for every section along the study reach.

The determination of the flood levels generated by ice jams over the study reach was carried out using an algorithm developed in Python making it possible to analyse all the simulations retained by the MOCA approach and generate the water level profile according to the 1:100 AEP.

Prior to the analysis of the scenarios, the algorithm first filtered the scenarios for which the simulated flow was lower than the minimum flow to initiate an ice jam, secondly, it filtered the scenarios for which the flow was greater than at which an ice jam can remain in place or collapse. The minimum and maximum flow values were determined based on available historical data and on in situ observations made during the winters of 2018, 2019 and 2020. The minimum flow for an ice jam to be initiated has been established to be 40 m³/s. Moreover, considering the variations in maximum flow values from one section of the river to another, a maximum flow threshold was used: 780 m³/s for the sections located between STL and SJ; 840 m³/s for the sections upstream of SJ to SJB; 1,200 m³/s from SJB to BCE; 850 m³/s for the towns of NDP and SGB. A total of 4,752 scenarios were simulated, i.e., 48 batches of 99 simulations. Of this number, nine of them were rejected due to model crashes. In addition, a total of five scenarios were removed since no ice jam was produced. The total number of simulations retained to determine the 1:100 AEP ice-induced water levels once filtered according to a minimum and maximum flow was 3,369 scenarios or 72.2% of the total scenarios produced.

In Equation (2), r is the rank of the water level in descending order and n is the number of values of the series. To consider the annual probability, the scaling factor of 53/129 was added to the denominator since there have been 129 ice jams recorded during the last 53 years, which gives an average of 2.6 events per year.

The 1:100 AEP ice-induced water levels for each town are presented in Table 5. Upon comparing the results with the historical analysis of water levels at each town, the correction (ΔH) applied to the values calculated by MOCA was done until an agreement was obtained with the historical analysis. A total of three out of the eight towns did not need a correction to be applied to the MOCA values. Certain external factors not considered in the Monte–Carlo approach may explain the correction (ΔH) applied. The MOCA model does not generate double ice jams. Yet, during winter observations in 2018, 2019, and 2020, winter ice breakup generated ice jams at various locations along the river that remained in place all winter. The spring breakup can then generate a double ice jam effect as observed in SJB (2020), BCE (2019) and SGB (2020). Moreover, the script of the MOCA model assumes a permanent flow (steady flow) along the reach and does not consider the dynamic aspect of the release of the ice, as well as the incidence of the effect of a jave on the flow. This is one of the limitations of the ice module in HEC-RAS and may cause increased uncertainty in the results. Using a dynamic model may have helped simulate unsteady flow and jave events, but the computational time and external control of the software would have been difficult to manage.

GSA is used to assess the influence of the different parameters of a simulation model. As described by Wagener et al. (2004) and Pianosi et al. (2016), GSA is used to determine the boundary conditions of a model and to assess to what extent each parameter can affect the output variables (X_GSA). This information can then be used to identify important parameters that significantly influence model response (Wagener et al. 2004; Jahandideh-Tehrani et al. 2020). GSA has been used in the delineation of flood zones and the resolution of topographic Lidar on flood mapping by ice jam, as well as that of the geomorphology (Lindenschmidt & Chun 2013; Das & Lindenschmidt 2021).

In the context of this research, the output variable of the simulation model is the water level reached following the formation of an ice jam simulated by the nine randomly generated mechanical and physical parameters. For each of the sections, the water level generated is classified by descending order for each ice jam. The top 10% of the SWS having the highest elevations, called the ‘behavioural’ series, is used as a reference to compare with the remaining 90% of the results, called ‘non-behavioural’ series, which represents, as adequately as possible, the effects of an ice jam at a specific location along the river, to assess the influence of each parameter on the output variable which is the SWS.

The GSA consists of performing a comparative analysis of the distribution obtained for the ‘behavioural’ series with the ‘non-behavioural’ series (Saltelli et al. 2008). The GSA, therefore, serves as a reference for comparing the results obtained from the model by varying the different parameters by a classification of sensitivity for all the parameters (Wagener et al. 2004; Lindenschmidt & Chun 2013).

The magnitude of the difference between the two distributions is estimated using a two-sample Kolmogorov–Smirnov (KS) statistical test. The greater the difference (α_GSA), the greater the influence of this parameter on the water level.

In the context of this study, a total of nine parameters were chosen, namely the location of the foot of the jam (XS), length of the jam (L_j), IT, Manning's ‘n’ of the ice (n), angle of friction (AF), porosity (p), K1, V_max and the flow (Q). The GSA was carried out on two sections of the Chaudière River, namely the section of SGB which has narrow flood plains with a steeper slope and the presence of an important tributary (Famine River) and the section of SJB which has wide floodplains with a low slope gradient. The behavioural series at both locations are the simulations which produced the top 10% maximum water levels, and the remaining simulations are classified as non-behavioural to analyse which parameters have the greatest influence on maximum water levels for each location.

As for the V_max parameter, its value in the KS test is α_GSA = 0.13 at SJB. The largest difference is found when V_max is close to 5.5 m/s. The influence behind this parameter is related to the hydraulic thickening of the ice cover of the jam. It could be the geometry of the river at this location that would be in question, such as the presence of a deep pool followed by a rapid rise of the riverbed. Concerning the five remaining parameters (IT, Mannings' ice roughness (n), ice porosity (IP), angle of friction (AF), K1), their influence on the sensitivity of the model is practically equal (α_GSA ≈ 0.08–1.2) so their variability does not play as important as a role as the parameters discussed previously.

From the results obtained for the 1:100 AEP ice-induced water levels, when compared to the open-water, certain locations across the study area are more prone to seeing higher water levels induced by ice than open-water. This is true for towns located at the upstream end of the reach. As for the towns located where the flood plains are at their largest (Figure 15: 41.8–85 km), the open-water flood level is higher than that of the ice-induced flood. Figure 15 presents the 1:100 AEP water level profiles for ice-induced and open-water levels. Along the superimposed profiles, a total of three subsectors have an ice-induced water level profile that is higher than open-water. The chainage of the subsectors is between 17.5 and 23 km (STL), 85 and 93 km (BCE and NDP) and 103.3 and 106.6 km (SGB). Table 6 presents the values for the 1:100 AEP for ice-induced and open-water levels. Values marked with an (*) indicate water levels induced by ice.

From a hydrological point of view, the lateral supply of the tributaries differs greatly upstream and downstream of VJ. The tributaries in the section upstream of VJ include several major tributaries that contribute significantly to the flow of the Chaudière River, including the Famine River (22.5%) and the Bras-Victor River (24%) including Plante River and Calway River. For the section downstream of VJ, the contribution of the eight tributaries varies from 3 to 6% each.

Moreover, in regards to the geomorphological aspect, when the slope gradient is significant between the tributaries and the Chaudière River, as is the case for the major tributaries of the study reach upstream of VJ, both water and ice supplied by these tributaries can significantly influence the ice dynamics in the main channel. On the one hand, the supply of frazil contributes to the formation of the ice cover and, on the other hand, the ice from the tributaries can abut against a still competent ice cover, which favours the initiation or consolidation of ice jams. It should be noted that such a phenomenon was observed during the three winters (2018, 2019 and 2020) without, however, evaluating the extent of the effects. Ettema et al. (1999) also report that confluences play an important role in the creation of ice jams and bring additional complexity to the flow at these locations. This is also pointed out for both locations where the GSA was done, downstream of a major tributary.

Regarding flood plains, the GSA for the study reach showed that the flow (Q) is the primary parameter to influence water levels when an ice jam was simulated at SJB. Since the maximum flow at which an ice jam can stay in place cannot exceed a certain threshold (section 4.1), the back-water effect will be routed by the flood plains. This phenomenon will therefore affect the influence of the mechanical parameters and low sensitivity since they are related to the hydraulic thickening of the jam (V_max). Therefore, the severity of the ice jam is limited by the imposing width of the flood plains, 300 m, and the flow will control the water level with minimal effect from the ice jam location (XS).

Continuing with the geomorphological aspects is the presence of islands or narrowing points. All subsectors that higher ice-induced flood water levels are in the presence of islands of varying importance, which range from a single island positioned in the middle of the channel to a braided channel, between the towns of NDP and STG, as well as for the towns of STL and BCE. In the case of BCE, an island is present in the middle of the channel and situated in the centre of the town, thus dividing flow, and reducing the hydraulic conveyance of the ice. In the case of narrowing points, the total width of the channel sees itself being reduced by up to 10% while seeing the flood plains metamorphosize to an elevated terrain with steep riverbanks will inevitably influence the hydraulic conveyance capabilities of the river when being obstructed by ice. As described in the GSA at Saint-Georges, the lack of flood plains combined with the funnel effect located 200 m downstream on the confluence of the Famine River comes into play for the hydraulic thickening of the ice jam imposed by V_max. The capacity for water to transit through laterally, via the plains, is close to none and the simulated water level will increase until the ice balance equation reaches equilibrium. Inarguably, ice jam length is also accountable for the back-water effect when positioned in such an area. The same explanation can be applied to BCE and STL, which both have elevated riverbanks.

A MOCA approach to determine the 1:100 AEP ice-induced water levels along an 80 km stretch of the Chaudière River was attained by a Python script controlling HEC-RAS. The MOCA framework was fed by historical ice jam parameters which included ice jam length, ice jam toe location, breakup flow and CDDF. By comparing historical water levels with the 1:100 AEP ice-induced water levels, a correction (ΔH) was applied to the values calculated by MOCA until an agreement was obtained with the historical analysis. This correction makes it possible to modulate the calculated values of the flood levels according to the external factors that were not considered by the MOCA analysis, which include double ice jams, unsteady flow and jave effects. Also, software limitations, such as unsteady flow ice jam simulation and jave effects, which are not handled by HEC-RAS, play a key role in the analysis of ice jam simulation.

Based on two separate locations along the Chaudière River, more precisely at SGB and SJB, a GSA was made to determine which ice jam parameters influenced the simulated water level in addition to discerning local geomorphological elements that affect local hydraulics. The most important hydrological parameters that made the 1:100 AEP ice-induced flood levels higher than 1:100 AEP open-water in this study were flow, ice jam length and location. Based on the topography of each site, it has been shown for this study reach, ice jam location and length are important factors when channel width is reduced while having elevated riverbanks. This will in turn influence the hydraulic conveyance of ice in the force balance equation of the jam, thus, affecting the sensitivity of ice jam parameters in such locations. When in the presence of wide flood plains, as in SJB, the flow will play the predominant role in maximum water level while seeing the ice jam length and location be of secondary effect. Further research is to be done in studying the effects of double ice jams in the AEP analysis for areas where mid-winter breakup occurs often and remains in place until spring breakup, such as the Chaudière River.

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

The authors declare there is no conflict.