RIVICE—A Non-Proprietary, Open-Source, One-Dimensional River-Ice Model

Abstract

Currently, no river ice models are available that are free and open source software (FOSS), which can be a hindrance to advancement in the field of modelling river ice processes. This paper introduces a non-proprietary (conditional), open-source option to the scientific and engineering community, the River Ice Model (RIVICE). RIVICE is a one-dimensional, fully-dynamic wave model that mimics key river ice processes such as ice generation, ice transport, ice cover progression (shoving, submergence and juxtapositioning) and ice jam formation, details of which are highlighted in the text. Three ice jam events at Fort McMurray, Alberta, along the Athabasca River, are used as case studies to illustrate the steps of model setup, model calibration and results interpretation. A local sensitivity analysis reveals the varying effects of parameter and boundary conditions on backwater flood levels as a function of the location of ice jam lodgment along the river reach and the location along the ice jam cover. Some limitations of the model and suggestions for future research and model development conclude the paper.

1. Introduction

River ice can be a major influence on the fluvial hydraulics and geomorphology of northern, high-latitude rivers. An important feature of river ice is that its presence on rivers raises water levels higher than those occurring during open-water conditions, given the same water discharge. The increased roughness and wetted perimeter length introduced by the ice contributes to the increased staging. Rubble or slush ice also has the potential to jam and dam up a river section, causing even more backwater staging and river banks and levees to overflow and flood the surrounding floodplain. From an ecological perspective, many floodplain and wetland ecosystems require such periodic flooding for their replenishment of moisture, sediment and nutrients, especially in perched ponds and lakes of inland deltas such as, for example in Canada, the Peace-Athabasca Delta in Alberta [1], the Slave River Delta in the Northwest Territories [2,3] and the Saskatchewan River Delta along the Saskatchewan/Manitoba boundary [4,5]. However, ice jams and the flooding they induce can pose threats to, and wreak damage on, many communities located along these rivers. To mention only a few cases in Canada, examples include the Town of Peace River along the Peace River in Alberta [6,7], Fort McMurray on the Athabasca River in Alberta (this study), Winnipeg and Selkirk on the Red River in Manitoba [8,9], St. Raymond on the St. Anne River in Quebec [10] and many communities along the Saint John River in New Brunswick [11]. Floods account for the greatest number of hydrological/meteorological natural disaster events in Canada [12] and, for most Canadian rivers, the annual peak water levels are due to ice jams [13]. Ice jam floods have also been reported on other continents, for instance along the Tornionjoki River in Finland and Sweden [14], the Oder River in Germany and Poland [15] and many rivers in Russia [16]. A review of successful river ice modeling cases is provided in [17]. Other important features of ice jam flooding are its unpredictability and the rapidity of inundation, often making it more dangerous than open-water flood events.

River ice modelling has been a useful tool in understanding key river ice processes that lead to ice jamming and subsequent flooding. The models help us to determine backwater levels, leading to a more comprehensive assessment of the flood hazard and risk posed to communities [6,7,18,19]. They assist in the design of flood protection measures such as crest elevations for river dikes [20], and in the mitigation or exacerbation of the effects of hydro-power dam construction and flow regulation on the severity and frequency of such events [21]. Also, much effort has been dedicated to attempting to forecast ice jams and ice jam flood events [22], with some success, and this remains a field with much research potential. Many of the processes leading to ice jam formation are probabilistic in nature (e.g., ice volume and location of an ice jam); hence, inserting such models in a stochastic framework, such as Monte-Carlo analysis, allows the random nature of these events to be captured and described in terms of frequency distributions. The model must be complex enough to capture all of the main processes leading to ice cover and jam formation, but not be so computationally demanding to discount the model’s applicability in a stochastic simulation environment. The River Ice Model (RIVICE) has proven to be very robust, capturing many of the essential river ice processes in a computationally efficient structure to yield high performing results. Examples include good simulation performances for flood water level profiles along the Peace River [7] and dike crest elevation design along the Dauphin River [20]. The model is a non-proprietary, open-source model currently housed at and distributed from the Global Institute for Water Security at the University of Saskatchewan, with permission from Environment Canada. User registration to access the software can be done at http://giws.usask.ca/rivice/. The author is not aware of any other river ice models that are free and open source software (FOSS), which may be a hindrance to the advancement of the field. Hence, an important aim of this paper is to promote a non-commercial option for scientists and engineers who wish to carry out river ice modelling.

In 1988, Environment Canada, along with several consulting firms, government agencies and a consortium of hydro-power companies, initiated the development of the River Ice Model with the aim of providing a non-proprietary numerical model to simulate river ice processes dynamically given a set of hydraulic, meteorological, fluvial geomorphological conditions and the thermal and ice regimes of the river water. These factors steer the simulations to capture key ice processes such as ice generation, ice transport, ice cover formation and progression, hanging dam formation and ice cover melting and ablation. The development of the model was completed by Maurice Sydor ([email protected]), then at Environment Canada, and Rick Carson from KGS Group in Winnipeg. Beta-testing was carried out by the author, then with Manitoba Water Stewardship for ice jam flooding on the Red River [8,23] and a flood risk mitigation project on the Dauphin River [20,24], both rivers in Manitoba, Canada. The first version of the model was released in January 2013. The model has since been applied successfully to winter flow testing of the Qu’Appelle River in Saskatchewan [25,26], ice-jam flood risk assessment and mapping of the Town of Peace River in Alberta [6,7], predicting the impact of climate change on ice-jam frequency at Fort McMurray on the Athabasca River in Alberta [27], determining areas of increased vulnerability to ice cover breakup from hydro-peaking along the South Saskatchewan River in Saskatchewan [28] and investigating ice jam characteristics in the Slave River Delta in the Northwest Territories [29].

2. River Hydraulic and Ice Processes Modelled in RIVICE

To convey an impression of an ice jam, a panoramic view of a jam is shown in Figure 1. The river and ice jam are small enough to allow the jam to be viewed in its entirety from one vantage point on the ground. The figure shows the toe of the jam lodged against a downstream intact ice cover. Bridge piers often slow down and arrest the flow of ice runs to create a rubble ice accumulation. The head of the jam shows the elevated water level surface due to the backwater effect imposed by the jam’s damming of water.

For an open water case, RIVICE solves the full dynamic wave of the Saint Venant equations in one dimension (hydraulic variation in the longitudinal direction; lateral and vertical variations are averaged), solved using an implicit finite difference scheme. To close the momentum and continuity equations of the Saint Venant formulations, Manning’s equation is applied:

$$v=\frac{1}{n_b}{r_H}^{2/3}{s}^{1/2}$$

(1)

which relates the mean flow velocity v at a cross-section to the river bed slope s, the cross-sectional hydraulic radius r_H (= cross-sectional flow area/wetted perimeter) and the Manning’s roughness coefficient n_b. RIVICE can also simulate the flow dynamics for an ice-covered river by increasing the wetted perimeter to include the width of the underside of the cover. For a solid ice cover, the roughness coefficient now becomes a composite value n_c of both the bed and ice roughness coefficients, n_b and n_i respectively:

$$n_c={\left(\frac{n_b^{3/2}+n_i^{3/2}}{2}\right)}^{2/3}$$

(2)

n_b is assumed constant for all in-stream hydraulic conditions of a river reach. For ice covers formed by ice rubble or slush the roughness n can be set as a function of ice thickness h through a relationship established by Nezhikhovskiy [30]:

$$n=0.0302*\ln (h)+0.0445$$

(3)

The equation is extended in RIVICE to allow the user to define the magnitude of the roughness coefficient at an ice thickness of 8 m (n_8m) [31]:

$$n=\frac{n_{8m}}{0.105}\left(0.0302*\ln (h)+0.0445\right)$$

(4)

The following summary of the ice processes modelled in RIVICE refers to Figure 2.

In the initial stages of river freeze-up, border ice often forms along the river banks, particularly in areas where flow conditions are more quiescent allowing a layer of skim ice and a subsequent thermally thickened bottom layer to form in the slower moving water. The border ice can extend towards the middle of the river from both banks, and can eventually meet mid-river, if flow conditions permit. In its simplest form, the user can indicate a maximum width and the simulation time for the border ice to attain that width in RIVICE. A dependency on cumulative degree days of freezing [32] or a more complex dependency on flow velocities and heat flux from the water to the atmosphere [33] are also options offered in RIVICE. The thickness of the border ice is defined by the user. If the water level drops or rises a certain input displacement, as defined by the user, the border ice breaks off and flows downstream to accumulate at the next ice cover front. A more detailed description of border ice formulation is provided in [34]. However, an increase in complexity does not necessarily yield more accurate results, due to the inherent uncertainty in the structure and parameters of border ice formulations. Nonetheless, the equations do provide an orientation to the magnitude of border ice volumes, which can be a large source of rubble ice for downstream ice cover formation. This is particularly important in wider rivers with low gradients of velocity currents along the banks and large fluctuations of the water level, where the border ice can potentially break off at a hinge crack along the bank and float as rubble ice downstream. Large fluctuations can occur through hydro-peaking of the flow by upstream dam operations, for example, or through the backwater staging from a downstream ice jam, particularly in low sloping river reaches.

The border ice covers can extend from both banks until they meet within the channel middle to form a bridging where the flow of ice formed upstream can lodge to initiate further ice cover formation. Lodgments may also occur at river constrictions, such as between bridge piers, at islands and at reaches that are narrower than the rest of the river. Other than through border ice bridging, the location of lodgments cannot be predicted by RIVICE but must be set by the user, indicated by the cross-section number. Efforts have been made, with some success through geospatial modelling, to correlate geomorphological features of a river system (e.g., sinuosity, channel width, slope) to lodgment locations [35,36] which may provide an initial orientation to where lodgments are most likely to occur.

Ice progression from a lodgment requires an upstream ice source to provide ice that can accumulate at the lodgment. One source, rubble ice broken off from border ice, was already discussed above. A volume of rubble ice per time step can also be explicitly provided as input, at the upstream boundary to represent any ice volume that has floated from upstream of the modelled river section. This may include fragments from upstream broken up ice covers or breached ice jams. Frazil in the form of slush pans is the third source of ice and is generated along the open water stretch upstream of the ice cover front, when the water temperature T_w is at 0 °C and the overlying air temperature T_a is below freezing. The heat transfer q simply becomes:

$$q=H\left(T_w-T_a\right)$$

(5)

with H, the heat transfer coefficient, typically ranging between 15 and 25 W/m²/°C, depending on how conducive the site conditions are to heat transfer. Many factors can affect the water-to-air transfer of heat, including wind speed, the degree of sheltering of the river from the wind (high sloping banks with trees provide more sheltering) and the amount of longwave radiation (cloudy conditions) that can be an important heat source that slows down the production of frazil ice. Only the open-water areas contribute to the heat loss; heat loss through moving ice is considered negligible.

The suspended frazil ice conglomerates to form frazil flocs, which become buoyant enough to float to the water surface. Here they further conglomerate to form slush pans that float downstream and add to the existing ice accumulation at the downstream lodgment. The main ice properties of the rubble ice and slush pans parameterised in RIVICE are their thickness ST and porosity PS.

Once the ice floes reach the lodgment or ice accumulation, the ice cover progression can occur in three ways—shoving, submergence and juxtaposition—depending on flow velocity and ice-cover forcing criteria:

Shoving: During the ice simulations, a balance of the external forces exerted on the ice cover is compared to the resistive forces within the cover. External forces include the thrust of the flowing water on the ice cover front F_T, the drag force exerted by the flow on the ice cover underside F_D, the vertical component of the weight of the ice cover along the direction of the river’s slope F_W, the friction of the ice along the river banks F_F, and the cohesive forces F_C along the banks when freezing of the ice to the banks occurs. Friction is parameterised by K1TAN which quantifies the amount the longitudinal compressive forces in the cover shed laterally to the river banks. The resistive force within the ice cover F_I is parameterised by the redistribution of the longitudinal forces vertically within the ice cover as it thickens K2. When the summation of the external forces exceed the resistive forces within the ice cover (F_T + F_D + F_W − F_F − F_C > F_I), the ice shoves downstream and thickens the cover. This “telescoping” of the ice cover continues until the ice is thick enough for the internal resistive force to exceed the summation of the external forces. Properties of the ice cover include its porosity PC, thickness of the ice front FT, and the thickness of the ice cover downstream of the lodgment h.

A quantitative description of the forces is provided elsewhere [24,31,37]. The cohesive and friction forces have been extended in the model to include further distribution of the force components to additional “banks” at cross-sections with islands or bridge piers:

$$F_C=2c·h·L·SHED$$

(6)

$$F_F=2f·h·L·K1TAN·SHED$$

(7)

where c is the area of cohesion at the ice-bank interface, f is the compressive stress within the ice cover, h is the average ice thickness, L is the distance between cross-sections and SHED is a shed factor indicating the number of bank pairs to which loads can be shed to increase the resistance to the external forces. SHED is set to 1 for no additional banks other than the left and right bank of the cross-section, 2 if an island/pier (four banks) is situated in the river, 3 for two islands/piers (six banks), and so on.

Submergence: If the densimetric Froude number F_r at the ice cover front is high enough, the rubble ice and slush pans submerge under the ice front and travel downstream in transit under the ice cover. F_r is a non-dimensional ratio of the inertial forces of the flow (represented by flow velocity v) to the gravitational forces (represented by the square-root of the gravitational acceleration g and the water depth at the ice cover front D). F_r can also be formulated as a function of the densities of water and ice, ρ and ρ_i respectively, the porosity of the ice cover PC and the ratio of the ice thickness h to the water depth D [38]:

$$F_r=\frac{v}{\sqrt{gD}}>\sqrt{2\frac{\left(\rho -{\rho }_i\right)}{\rho }\left(1-PC\right)\frac{h}{D}}\left(1-\frac{h}{D}\right)$$

(8)

The ice can then deposit on the ice underside if the mean flow velocity under the ice drops below a threshold value v_dep. Deposited ice can also erode and be transported downstream again if the flow velocity exceeds a higher threshold value v_erode.

Juxtaposition: If the Froude number is low enough and the ice cover’s internal resistance exceeds the summation of the external forces applied to it, then the rubble ice and slush pans will stack up against each other along the water surface and protract the ice accumulation in the upstream direction.

Subroutines are also embedded in the RIVICE source code to simulate water quality. Variables that can be simulated include salinity, water temperature, dissolved oxygen, biochemical oxygen demand, organic and inorganic components of nitrogen and phosphorus, phytoplankton, zooplankton, fecal coliforms and conservative and decaying lignins. Reaeration, substance decay rates, phytoplankton respiration and nutrient-limited growth, zooplankton grazing, excretion and respiration, are some of the parameters that steer the water-quality simulations. Further model development and testing is required for the coupling of the water quality components to the river ice processes

3. Study Site

A map of the site, Athabasca River at Fort McMurray in Alberta, Canada, is provided in Figure 3. Fort McMurray is particularly vulnerable to ice jam flooding since the Athabasca River’s geomorphological setting at that location makes this river section conducive to ice jamming, for the following reasons:

• There is a substantial change in river slope (from approximately 0.0010 to 0.0003 m/m) and widening of the river cross-section (from approximately 300 to 700 m) immediately at the town;
• There are many rapids along the upstream stretch to the Town of Athabasca (approximately 300 km) along which much ice rubble and water accumulates which, once released, can discharge a surge of water and ice into the Fort McMurray area;
• The confluence of Clearwater River lies at Fort McMurray and can be a source of additional ice and water, exacerbating an already hazardous ice flood situation along the Athabasca River (see Figure 4); water can also back up in Clearwater River from a jam on the Athabasca River to flood Fort McMurray’s downtown area;
• Many obstacles, both natural (e.g., islands and sand bars) and anthropogenic (e.g., bridges) in origin, can impede the flow of ice, increasing the area’s susceptibility to ice jam formation;
• The Athabasca River’s flow is in a northerly direction, an orientation favourable to ice jam formation, since the ice cover in the upstream, southerly areas generally breaks up sooner due to warmer weather conditions (increased ice cover ablation and runoff from melting snow), transporting ice to downstream, northerly areas where the temperatures are cooler and the ice cover remains intact and competent longer to resist the flow of ice.

Beltaos tested his model assuming a rectangular channel using approximations of the ice jam parameters typical for breakup jams at Fort McMurray [39]. Other modelling exercises have been carried out along the Athabasca River but with a different focus (ice jam releases) in a different reach along the Athabasca River (between the Town of Athabasca and Fort McMurray, a stretch that extends upstream from the one studied here) [40,41]. The upstream stretch is much steeper and generally geomorphologically different than the one presented in this paper.

4. Model Setup

Three ice jam events were investigated in each of the years 1977, 1978 and 1979, to lay out the steps for setting up the RIVICE model. The configurations of each jam, in regard to the river hydraulic and ice regimes and the extent of the subsequent flooding, are quite different, allowing the performance of the RIVICE model to be tested.

Bathymetry: Cross-sections of the river are required to capture the cross-sectional flow areas along the river, an important pre-requisite for accurate simulations of flows, current velocities and water level elevations. An example of a cross-section is shown in Figure 5, indicating the difference in width and flow area between the steeper stretch immediately upstream of Fort McMurray and the less steep stretch immediately downstream of the town. The cross-sections consist of surveyed elevation points along transects of a river bed.

Lateral inflows: Lateral inflows, which include tributary inflows and overland runoff, are simulated as a discharge per unit length along a river stretch. In our test case, the inflow from Clearwater River is incorporated as its discharge is divided up along the 100 m length, which corresponds to the approximate width of its mouth.

Boundary conditions: Typically, water flow rate serves as an upstream boundary condition whereas a water level elevation designates the downstream boundary condition, although rating curves can also define the boundaries. These values were derived from the records of the Water Survey of Canada gauges “Athabasca River below McMurray” (# 07DA001) and “Clearwater River at Draper” (# 07CD001).

Surveyed water level profiles: During the ice jam events, water levels were surveyed along the banks of the river to provide a longitudinal profile of the flood water level elevation [42]. In difficult to access areas near the jam, aerial photographs were taken of flood water along the banks with reference features that could then be surveyed after the event had passed. A similar method was carried out by the author along the Dauphin River in Manitoba using a road surface adjacent to the river as a reference elevation for the encroached ice jam flood water, as shown in Figure 6. A road surface profile was surveyed before the flood event to provide backwater level elevations during the flight.

Parameters: Important parameters used for model calibration can be grouped according to the characteristics or function of a component of the hydraulic/ice regime (

Table 1. Parameters and boundary conditions for model calibration.

Parameter	Description	Units	1977	1978	1979
rubble ice:
PS	porosity of rubble ice	–	0.3	0.3	0.3
ST	thickness of rubble pans	m	2.5	1.2	1.5
ice cover:
Vice	volume of incoming ice	m3/Δt	2000	2000	2000
PC	porosity of ice cover	–	0.5	0.5	0.5
FT	thickness of ice cover front	m	2.5	1.0	1.5
ice jam lodgment:
h	thickness of ice downstream of jam	m	1.0	0.75	1.0
x	cross-section number of lodgement 1	–	550	491	783
ice transport:
vdeposit	threshold of deposition velocity	m/s	1.2	1.2	1.2
verode	threshold of erosion velocity	m/s	1.8	1.8	1.8
hydraulic roughness:
n8m	ice roughness	s/m⅓	0.14	0.16	0.11
nbed	river bed roughness	s/m⅓	0.025	0.025	0.025
strength properties:
K1TAN	lateral: longitudinal stresses	–	0.18	0.18	0.18
K2	longitudinal: vertical stresses	–	8.4	8.4	8.4
boundary conditions:
Q	upstream discharge	m3/s	1000	1500	1366
W	downstream water level	m a.s.l.	235.0	235.0	235.5

Note: ¹ corresponding chainage = 27,500, 24,500 and 39,150 m for 1977, 1978 and 1979, respectively.

). Parameter values and ranges were extracted from the literature [42,43] and then calibrated for each event to fit the simulations to the water level observations.

5. Results and Discussion

5.1. Calibration and Validation

Figure 7 shows the longitudinal simulated profiles of the ice cover (black) during the ice jam events of each year studied; 1977, 1978 and 1979. Included, for reference, is the thalweg (orange) and the open-water profile (blue line) attained with the same upstream and tributary discharge and downstream water level boundary conditions. The open-water profile reveals how much staging occurs by the ice cover and jam. A Manning’s roughness coefficient n_b of 0.025 was calibrated to fit the simulated open-water level elevation at the gauge with the water level elevation recorded by the gauge (blue dot).

Surveyed water level elevations attained at the peak of each ice jam flood event (pink diamonds) were included to guide the calibration of the remaining river ice parameters and assess model performance. The optimum values are provided in Table 1. Generally, there is good agreement between the ice jam water profile and the surveyed water elevation points. The 1978 ice jam event was difficult to calibrate and yielded larger deviations between simulated and observed water level profiles. This is due to the fact that the discharge upstream of the jam (1850 m³/s, sensu [42]) was far greater than that recorded at the gauge downstream of the jam (544 m³/s). This is a large discrepancy in flows. To fulfil continuity, the extra upstream water may have been diverted around the jam, which has not been considered in the model. A similar diversion was reported during the 1977 event by Doyle [43], who stated that “flood wave [washed] over ice on right side of island downstream of bridges but ice sheet remains intact while sheet on left side of island is broken up and all flow is directed around left side of island”. Increasing the discharge to fit the simulated ice jam water profile to the surveyed water level points required the ice cover downstream of the jam toe to lift slightly in order to allow for additional flow under the jam toe and maintain computational stability. Using an upstream boundary discharge of 1500 m³/s, a “Pareto” balance was attained between slightly underestimated upstream water level profiles and a slightly overestimated ice cover elevation downstream of the ice jam. There was some discrepancy between the upstream discharge and downstream gauge flow reading for the other studied events, but not of the same magnitude as the 1978 ice jam. Andres and Doyle [42] estimated a discharge range of 1135–1600 m³/s upstream of the 1977 ice jam, with 934 m³/s recorded at the gauge. For 1979, the upstream flow was estimated to range between 1300 and 1850 m³/s when the gauge reading showed 1480 m³/s.

The 1978 ice jam event is a special case because its toe was situated at the bridge piers, which provided additional support and strength to the jam to maintain its high volume of rubble ice and high backwater staging levels. Some of the extra strength required to support such a large jam can be incorporated in the model using an increased SHED factor of 5, to account for the bridge piers in the water, which shed much of the longitudinal compressive stresses in the ice cover laterally to the banks and piers. However, not all of them may have been shed, and some of the longitudinal stresses of a jam of that size were potentially shed vertically onto the banks and, perhaps, on points along the jam where the ice grounded on the river bed or a sand bar. This shedding of forces has not been included in the model and is a topic to be tackled in future model development.

5.2. Local Sensitivity Analysis

Each jam lodgment was also situated at different locations along the river, each location having slight differences in their fluvial geomorphology. The 1977 jam occurred the furthest upstream in the steepest section of the studied reach and the 1979 ice jam was located furthest downstream in the least steep slope of the reach. These differences in the jams’ geomorphological settings will have different impacts on the river ice processes and parameters. Such effects can be quantified using a sensitivity analysis to show how simulated backwater levels react to slight changes in the parameter and boundary condition settings.

Figure 8 provides the local sensitivities of each parameter to the output variable backwater level elevation, both 5 km and 15 km upstream of the ice jam lodgment of each studied event. Positive sensitivity means that an increase/decrease in the parameter yields an increase/decrease in the output variable. When a parameter induces a negative sensitivity to the output variable, the variable increases or decreases opposite to the value change in the parameter.

Rubble ice: Ice jam formation is not sensitive to the parameters related to the ice rubble characteristics, PS and ST, regardless of the location along the jam. Any rubble or slush ice simply adds to the accumulation of the ice volume in the jam, regardless of the dimensions of the ice blocks or pans.

Ice jam cover: Backwater levels 15 km upstream of a lodgment are positively sensitive to the inflowing ice volume V_ice and negatively sensitive to the ice cover porosity PC. Hence, an increase in the compaction of the ice in the jam due to more ice of lesser porosity will lead to increased water level staging far upstream of the jam. However, in close proximity to the jam lodgment (5 km upstream from the jam toe), the sensitivities reverse due to the added weight of the upstream ice abutted to the lodgment. Staging all along the ice jam is positively sensitive to the ice front thickness FT since more thrust force against the cover causes more shoving against the jam ice cover and subsequent backwater staging.

Lodgment: An indication as to the thickness of the intact ice cover h before the ice jam events is provided in [42]. However, variations in the thickness of the intact ice cover downstream of the ice jam have little effect on the staging upstream of the jam. Geomorphology is also an important influence on floodwater levels, as reflected in the variation of sensitivities of the ice jam lodgment location x. This highlights the importance of accurately representing the bathymetry in the model setup.

Ice transport: The ice transport parameters, v_d and v_e, have little effect on the ice jam morphology or water levels since the main mechanism of ice cover thickening is the retraction of the ice cover by shoving.

Hydraulic roughness: Backwater level rises are very sensitive to the ice jam underside roughness parameter n₈, more so than river bed roughness n_b, particularly in close proximity to the ice jam lodgment. The point closer to the lodgment (5 km upstream of the ice jam toe) induces a negative sensitivity of the parameters K1TAN and K2 on backwater levels, indicating the importance of the frictional and resistive forces in reducing further thickening and stabilizing the ice jam. Increases in these parameters will increase friction of the ice rubble against the banks and amongst each other within the ice accumulation. This becomes less prevalent at the point further upstream of the lodgment (15 km upstream of the ice jam toe) where sensitivities lessen and, particularly for the 1977 event, even become slightly positive.

Boundary conditions: Flow Q is one of the most sensitive factors influencing water level elevation and becomes more sensitive as we move further upstream along the jam. More flow will induce greater thrust and drag forces on the ice cover. Further upstream, the river slope is steeper, which increases the weight of the ice on the cover to promote shoving. The downstream water level boundary W only has an influence on the 1977 water levels.

6. Discussion on the Broader Applicability of RIVICE

The applicability of RIVICE has been shown here using one case study, that of ice jam flooding along the Athabasca River at Fort McMurray at latitude 56.7° N. The model has also been successfully applied to a more northerly river, the Slave River near Fort Resolution in the Northwest Territories of Canada at latitude 61.3° N and more southerly regions such as the Red River at Winnipeg and Selkirk in Manitoba at the approximately latitude 50° N. The mean winter flows for the months between December and February of the Athabasca River at Fort McMurray is 174 m³/s but the model has also been successfully applied to rivers with mean January to February flows as high as 1123 m³/s, such as at the Town of Peace River on the Peace River and as low as 4 m³/s for the upper Qu’Appelle River in Saskatchewan. Although these examples all stem from Canadian rivers, more recently RIVICE has also been applied to the Oder River, the lower reach of which constitutes the border between Germany and Poland. Although these examples all stem from regions with continental climates, currently RIVICE is being implemented for the Exploits River in Newfoundland, a region which has a more maritime climate.

Due to its relatively short computational times and the availability of the source code, RIVICE can be embedded in a Monte-Carlo framework to allow stochastic hydraulic modelling of ice jam flood events. Such an approach is well suited to ice jam modelling because ice jamming is a stochastic process and can best be described in probabilistic, not deterministic, terms. The stochastic approach allows the author to currently engage in research to develop ice jam flood forecasting systems in which RIVICE serves as an integral part in providing ensembles of backwater level predictions to determine probabilities of flood level exceedances. Such forecasting systems are being developed for the Athabasca River at Fort McMurray and the Exploits River at Badger under the project “Near real-time Ice-Related Flood Hazard Assessment (RIFHA)” funded by the Canadian Space Agency and managed by C-Core, a remote sensing company in St. John’s, Newfoundland. Under the Global Water Futures research program at the University of Saskatchewan (http://gwf.usask.ca/), funded in part by a $77.8-million grant from the Canada First Excellence Research Fund, ice jam flood forecasting systems are in planning to be developed for the Porcupine River at Old Crow in the Yukon Territories and the Red and Assiniboine rivers at Winnipeg, Manitoba.

7. Conclusions

RIVICE captures many of the major river ice processes within the channel that can lead to ice jamming and subsequent flooding. Flood waters that overtop the river banks and divert around the jam are not accounted for; however, developments are underway to couple RIVICE to a two-dimensional storage cell floodplain model that will allow water to be exchanged between the river channel and the floodplain. Future studies should also concentrate on extending the model to incorporate the additional forces that support the ice jam vertically when it is grounded at the river banks or on shallower points of the river bed such as at sandbars.

The local sensitivity analysis shows varying importance of the parameters and boundary conditions to the ice jamming. Both the location of the ice jam lodgment and the position along the jam ice cover show differences in the sensitivities to backwater staging. An example is the varying sensitivities to staging along different locations of the jam, imputed by the parameters associated with the frictional and internal resistive forces and properties of the ice jam accumulation (e.g., porosity).

One limitation of RIVICE is that it can only form ice covers through the juxtapositioning of frazil-generated or rubble-transported ice. Hence, it is only suited for rivers where mean flow velocity, at freeze-up or ice jamming, exceeds 0.4 m/s, a common threshold between thermally- and mechanically-driven ice formation. Additionally, a thermal module still needs to be integrated into the model to allow for thermal thickening of the ice cover once the cover has been established along the river. This is a topic of future research and model development.