Development of an ice-jam flood forecasting modelling framework for freeze-up/winter breakup

Ice processes are crucial hydrological factors for many rivers in cold regions. The presence of ice in lotic waters alters the flow conditions, water quality, and geomorphological and ecological settings of rivers (Prowse 2001). One of the more extreme river ice processes, ice-jam formation, can create rapid water level fluctuations and severe flooding, which have caused millions of dollars of property and business losses to riverside communities. Moreover, the dynamic behaviour of ice processes can interrupt hydropower generation, disrupt hydrometric data recording, and affect navigation. Ice-jams are capable of significant erosion of river beds and banks. Furthermore, ice-jams can cause disturbances within the aquatic environment and adversely impact aquatic habitats (e.g., increased fish mortality and destroyed spawning grounds). River ice dynamics also control dissolved oxygen concentration and can change the mixing process within the water column by scouring transporting and depositing sediment.

Ice-jams often occur during spring breakup; however, winter ice-jams are also frequent in various rivers in cold regions (De Coste et al. 2022). Winter ice-jams can occur during freeze-up and stable ice cover periods and are often caused by brief thaws, rain-on-snow events, and high flow conditions along the river (Beltaos 2002). Winter ice-jams are often unanticipated and difficult to manage, as subsequent cold air temperatures freeze the flooded water in place for the entire winter period. Winter ice-jams can form relatively thick and competent ice covers and intensify the following spring breakup conditions and ice-jam formations. Many recent studies have reported an increasing trend in winter breakups throughout North America (Carr & Vuyovich 2014; Das et al. 2022a, 2022b; De Coste et al. 2022). Newton et al. (2017) reported a total of 52 mid-winter breakup events in Western Canada and Alaska since 1950. In recent years, winter breakup has become a common occurrence along the Saint John River in Maine, USA and New Brunswick, Canada, where 27 freeze-up/winter breakup events have been observed since 1981 (Das et al. 2022a, 2022b). A similar increasing trend has been observed along other rivers in the United States, such as Picataquis River in Maine (Huntington et al. 2003), the Fox River in Wisconsin, and the Grand River in Michigan (Carr & Vuyovich 2014). A mid-winter breakup has also been observed on the Klondike River, Yukon (Janowicz 2010).

Since hydroclimatic factors are the major drivers of the occurrences of winter breakup, current changing climatic trends can greatly influence the pattern and severity of this type of event in the future (Beltaos 2002; Turcotte et al. 2019; Burrell et al. 2022). Several studies suggest that a warmer climate will result in a higher frequency of winter breakup occurrences and increase the risk of winter flooding along many communities in cold regions (Das et al. 2022a, 2022b; Lamichhane et al. 2022). Therefore, it is very important to develop a forecasting framework to examine the severity of winter breakup to reduce flood impacts and better prepare for emergency measures. Although many studies have already focused on spring ice-jam flood forecasting, no study has focused on winter ice-jam flood forecasting.

Many tools have been introduced to forecast river ice-jam conditions, such as empirical approaches, logistic regression, discriminant function analysis (White 2003, 2008), and multiple linear regression (Mahabir et al. 2006, 2008). However, these tools are limited to establishing a relationship between hydrometeorological variables and providing the occurrences of ice-jamming and breakup dates rather than estimating the backwater level of the river during ice-jam formation. Moreover, these methods are site-specific and challenging to implement in different geomorphological settings. Some data-driven machine learning models have also been introduced to predict breakup conditions (dates and flows), such as the k-nearest neighbour algorithm (Sun & Trevor 2017), decision-tree models (Sun 2018), neural networks (Mahabir et al. 2008; Wang et al. 2008; Guo et al. 2018), and fuzzy-logic models (Sun & Trevor 2015; Zhao et al. 2015). The combination of one or two of these models has also proven to produce better prediction accuracy (Sun & Trevor 2017, 2018a, 2018b). Although these models produce useful results, the lack of interpreting physical processes between the inputs and outputs makes it challenging to decision-making during emergency situations. Moreover, the data-driven models highly depend on the data availability and quality. Some attempts have also been carried out by applying various deterministic models to forecast flood water levels for ice-jam, such as River2D (Brayall & Hicks 2012), HEC-RAS (Beltaos et al. 2012), and RIVICE (Lindenschmidt et al. 2019). However, properly parameterizing a deterministic model to forecast ice-jam flooding is often challenging, which requires many parameter and boundary conditions, such as streamflow forecast, ice-jam locations, and ice volume. In recent years, Lindenschmidt et al. (2019) developed a stochastic framework for real-time ice-jam flood forecasting, combining hydrological and hydraulic models to produce an ensemble of ice-jam flood water level profiles during the progression of spring breakup. This approach was then tested in various rivers in Canada, such as Lower Red River, Manitoba (Williams et al. 2021) and Saint John River, New Brunswick (Das et al. 2022a, 2022b). This prominent flood forecasting approach has been applied in this study for mid-winter ice-jam flood forecasting.

The main purpose of this study is to develop a framework for forecasting mid-winter ice-jam flooding. The specific objectives are (i) to analyze mid-winter breakup along the Saint John River from Dickey, Maine, USA to Grand Falls, New Brunswick, Canada; and (ii) to develop a stochastic framework for mid-winter ice-jam flood forecasting.

Hydroclimatic data (e.g. streamflow and temperature) were analyzed and used to examine the winter breakup along the river. Daily mean streamflow data were obtained from the hydraulic gauge stations at Fort Kent (#ID: USGS01014000), Edmundston (#ID: WSC01AD004), and Grand Falls (#ID: WSC01AF002). The gauge station at Fort Kent is operated by the U.S. Geological Survey (USGS) (https://waterdata.usgs.gov/nwis) and Edmundston and Grand Falls are operated by Water Survey Canada (WSC) (https://wateroffice.ec.gc.ca/) (Figure 1). It should be noted that ice-related data are occasionally unavailable because ice events often damage the gauge station and stop data recordings during the events. Observed mean daily air temperature data were obtained from meteorological stations (Climate ID # 8101303) at Edmundston, New Brunswick, Canada.

The meteorological forcing data (wind speed, air temperature, specific humidity, incoming longwave radiation, barometric pressure, and incoming shortwave radiation) required for MESH simulations were obtained from the Regional Deterministic Prediction System (RDPS, Caron et al. 2015), which has spatial and temporal resolutions of 10 km and 1 h, respectively. RDPS is a part of the Global Environmental Multiscale (GEM) numerical weather prediction model (Côté et al. 1998), developed and operationally run by ECCC, and provides short-term forecasts (3.5 days). For model calibration and validation purposes, meteorological forcing data were archived from the first day of the forecast horizons (0–23 h) from the RDPS system. However, the precipitation data were taken from the Canadian Precipitation Analysis (CaPA) (Mahfouf et al. 2007) which are available daily at 6-h time intervals at a spatial resolution of 10 km. CaPA is a reliable precipitation product for the Canadian domain (Boluwade et al. 2018).

Hydrological forecasting for 2018 was performed using the meteorological products from the Global Deterministic Prediction System (GDPS, Buehner et al. 2015). Like RDPS, GDPS is also part of the GEM numerical weather prediction model and is a global operational forecasting system that provides deterministic atmospheric variable predictions with a 10-day lead time (Bélair et al. 2009; Charron et al. 2012). The forecasts are produced twice a day (00 UTC and 12 UTC) at 3-h temporal resolution and have a spatial resolution of approximately 25 km. Previous studies (e.g., Lindenschmidt et al. 2019; Rokaya et al. 2020) have found GDPS to provide reasonable results. All the forcing datasets (CaPA, RDPS and GDPS) were downloaded from the Canadian Surface Prediction Archive (CaSPAr) platform (Mai et al. 2020) (https://caspar-data.ca/caspar).

Historical winter breakup information was obtained from the Cold Regions Research and Engineering Laboratory (CRREL) ice database (https://icejam.sec.usace.army.mil/ords/icejam/r/icejam/home). The information from this database such as breakup dates, jam locations, and a brief description of the breakup and ice-jams from this database was used to examine the winter breakup phenomenon and modelling setup. Moreover, the study site's visual history of river ice processes is recorded by the Department of Environment and Local Government, New Brunswick. The winter breakup and associated ice-jam location data were analyzed from these visual history records. Furthermore, the progress of ice cover formation during the freeze-up was captured using Sentinel-1 satellite imagery from the sentinel hub through the EO browser (https://www.sentinel-hub.com/explore/eobrowser/). EO browser is an open-source database that provides an image every 7 days, allowing us to capture the ice cover progress along the river.

MESH (Modélisation Environnementale Communautaire – Surface Hydrology) is a semi-distributed hydrological land surface model developed by Environment and Climate Change Canada (ECCC) (Pietroniro et al. 2007). The hydrological processes represented in MESH consist of three components: (i) vertical exchange of water between soil, plants and atmosphere within a grid cell using CLASS (Verseghy 1991; Verseghy et al. 1993), (ii) generation of surface and subsurface runoff within a grid cell, and (iii) routing of lateral fluxes within the stream network using WATROUTE (Kouwen 2016).

In the soil profile, the moisture content and its vertical movement in the soil layers are calculated by solving Richard's equation. As a result, CLASS estimates overland runoff, interflow, and baseflow for each computational grid cell. The horizontal flow for each Grouped Response Unit (GRU) consists of overland flow, interflow, and baseflow. Overland flow is calculated using Manning's approximation of the kinematic wave propagation. Interflow, which accounts for the downward hill slope flow, is estimated from the bulk saturation of each soil layer, calculated at each time step. Baseflow is treated as any water that percolates out of the bottom of a soil column in a GRU and is immediately added to two interconnected hydrological reservoirs that slowly release water towards streams. The total amount of surface runoff is then calculated for each grid cell by the areal average of GRUs within that grid cell. More details on MESH, including the recent advances, is available from Wheater et al. (2022).

The stochastic modelling approach can simulate hundreds of probable scenarios under many conditions using a set of random variables. In this study, the RIVICE hydrodynamic model was embedded into the Parameter Estimation Software (PEST) program to perform the Monte-Carlo (MOCA) method. The RIVICE model is able to simulate various river ice processes, such as frazil ice generation, transport, juxtaposition and hanging dam formation, shoving and ice cover melting and ablation using a set of river ice and hydraulic parameters and boundary conditions (EC 2013b; Lindenschmidt 2017). The boundary conditions include upstream discharge (Q), downstream water level elevation (W), location of the ice-jam toe (x), and volume of inflowing ice (V_ice). Key parameters include porosity and thickness of rubble pans (PS and ST, respectively), porosity and thickness of the ice cover at its front (PC and FT, respectively), thickness of the intact ice cover (h) downstream from the ice-jam, depositional and erosional velocity thresholds (v_d and v_e, respectively), ice and river bed roughness (n_8m and n_bed, respectively), and longitudinal to lateral and longitudinal to vertical stress distributions within the ice-jam cover (K1TAN and K2, respectively).

ECCC's Green Kenue software (EnSim Hydrologic, 2014) was used to prepare the drainage database required for MESH. For this, the Digital Elevation Model (DEM) was retrieved from the hydrologically adjusted elevation of MERIT Hydro database (Yamazaki et al. 2019). The landcover and vegetation parameters were retrieved from the Commission for Environmental Cooperation (CEC) North American landcover (30 m resolution) and CLASS manual (Verseghy, 2009), whereas the soil texture information for different soil depths were retrieved from the Unified North American Soil Map (UNASM) (LIU et al. 2014). The model was set up with the spatial grid resolution of 0.125° resulting in 373 grid cells and having a drainage area of 41,000 km² with the outlet at the Mactaquac Dam station (long./lat.: −66.672, 45.951). Further details on the model setup are available from Budhathoki et al. (2022).

OSTRICH (Optimization Software Toolkit for Research Involving Computational Heuristics) (Matott 2005), an open sourced auto-calibration optimization software, was used for model calibration. Ten parameters of six dominant GRU's were calibrated for the period of 2002–2010 (model spin-up from Oct 2002 to Oct 2003) using the Nash–Sutcliffe efficiency (NSE) and its logarithmic value (LogNSE) as objective functions.

The modelling framework still requires many assumptions and is based on historical observation data, which may create some uncertainty in the results. For example, the locations of the ice-jam toe locations may not always be in the right stretch; therefore, it may create some level of overestimation or underestimation in the results. To overcome this limitation, real-time monitoring data or field observation can be used to select the probable toe location during freeze-up. Moreover, the study only simulates ice-jam scenarios; therefore, the framework can only produce accurate water level elevation when there is an ice-jam along the river. The ice volume was used from the simulation results; however, more robust study and field data are required to estimate the total volume of ice cover for a freeze-up or winter periods. Furthermore, the hydrological model did not consider the ice processes along the river; therefore, there is some uncertainty in the forecasting flow results.

In this study, a first attempt was carried out to develop a framework to forecast ice-jam flooding along the Saint John River. The stochastic modelling framework was developed to forecast the severity of mid-winter breakup along the Saint John River. The study combines the hydrological and hydraulic models and multi-sources of data such as satellite, gauge, and historical studies to develop the stochastic framework. The MESH hydrological model provides the forecasted flows, and river ice modelling results provide another important boundary condition, the volume of ice cover. Overall modelling results were applied to a real case study or mid-winter breakup in December 2018. The framework was successfully able to model the ice-jam scenarios and backwater level elevation along the river.

This stochastic modelling framework can be used to forecast mid-winter breakup severity in an operational context. This can help to make appropriate management decisions and emergency measures before the breakup. The stochastic outcomes allow the assessment of hundreds of probable scenarios that could be used in decision-making during emergency situations. The percentile results from these hundreds of forecasts could be used to determine the level of ice-jam flood risk for the study area. For example, if most of the simulated water level profiles between 50th and 95th percentile exceed the historical flood water level elevation, it could be marked as a high risk/emergency. While the 50th percentile water level profile is below the historic flood level, it can be indicated as a low-risk situation. However, further studies are required to make the framework more robust in changing climate. To do this, real-time time freeze-up observation and hydrological forecasting model for each river basin are needed to obtain the forecasting modelling data.

The authors thank the Global Water Future program at the Global Institute for Water Security, University of Saskatchewan for their funding support of this research.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.