Application of the HBV model for the future projections of water levels using dynamically downscaled global climate model data

The long-term increase in air temperature globally has modified the energy and water fluxes, and thus the hydrological cycle. Rising air temperature may accelerate the atmospheric moisture transport which could alter the climate system (Newton et al. 2014). Changes in the hydrologic cycle would occur as temperature increases, and observations along these lines have already been documented (Groisman et al. 2004; Pryor et al. 2009). Changes to seasonal air temperature and precipitation can have a significant long-term impact on the hydrology of river basins across the world, including that of Canada. The Mackenzie River Basin (MRB), the largest River Basin in Canada, is expected to be subject to the significant impact of climate change. Several studies have pointed out that the MRB is experiencing some of the greatest rises in temperature anywhere in the world (Wayland 2004), and the studies also have concluded that over the last few decades, the warming trend is going up sharply in Canada including the MRB (Shabbar et al. 1997; Zhang et al. 2000; Aziz & Burn 2006). Conversely, we expect that climatic change will modify the future water resources of MRB, especially during summer when there will be more demand for water to adequately maintain aquatic and terrestrial ecosystems, agricultural production, northern ferry operations, and hydroelectricity generation (Stewart et al. 2004; Trenberth et al. 2007).

Hydrological simulations to predict the future streamflow at the basin scale require future climate data. This future rainfall and temperature data can be obtained from the Global Climate Models (GCMs) outputs. However, the GCMs climate data resolutions are too coarse for direct use in catchment-scale hydrological modeling and need to be downscaled before the simulation process (Fowler et al. 2007). The downscaled GCM data could be a reliable way to simulate the streamflow of a river using any hydrological model like the Hydrologiska Byråns Vattenbalansavdelning (HBV) model (Driessen et al. 2010). The HBV model was developed by the Swedish Meteorological and Hydrological Institute (SMHI) in 1972, and is a widely used hydrological model. Its capability in conducting hydrological analysis related to the streamflow simulations is well known. Many studies have been reported around the world to address the potential impact of climate change on future hydrology and water resources (Driessen et al. 2010; Hirpa et al. 2019; Pokhrel et al. 2019; Chen et al. 2020), but for a large river like Makenzie River to predict the future streamflow based on climate change scenarios still remains challenging.

In this study, the Weather Research and Forecasting (WRF) model was selected to dynamically downscale the GCM data into a regional scale. The downscaled temperature and rainfall data will be the input to the hydrology model (HBV) to simulate the streamflow. The WRF model is a next-generation mesoscale numerical weather prediction system that is being used by many researchers in different parts of the world. Several studies have reported that the WRF model could reliably downscale the future temperature and rainfall data into much finer scale. It has been utilized to simulate rainfall over various regions (Chawla et al. 2018; Yang et al. 2019; Huang et al. 2020). Again, the applications of the WRF model to simulate air temperature were reported by several researchers in distinct locations (Fonseca et al. 2019; Wyszogrodzki et al. 2020). There is numerous evidence that the WRF is capable of simulating long-term climate data: Li et al. (2019) utilized the WRF model to downscale the long-term climate data and Ishida et al. (2020) used the WRF for future climate projection in Northern California. Knist et al. (2020) analyzed the projected changes of precipitation for 2038–2050 and 2088–2100 with respect to 1993–2005 derived from the WRF downscaled global climate runs based on an RCP 4.5 scenario. And Chotamonsak et al. (2011) projected climate change over Southeast Asia using the WRF model for 2045–2054. There is considerable evidence that the WRF model has been used in many climate change research studies. To simulate the streamflow using a hydrological model, finer resolution data are needed which could be obtained from the WRF model's dynamic downscaling process. However, there is not much evidence to use the WRF model data into a conceptual hydrological model to simulate the streamflow. Considering this fact, we incorporated the WRF model data into the HBV model to simulate the streamflow and to project the future water levels of the Mackenzie River.

The Mackenzie River possesses a vital role in the northern transport system. In an open water season, the ferry transport system keeps the northwest part connected to the rest of Canada. Here, the water level in the summer season is important for planning and managing the ferry operation system, which could be affected by the recent trend of climate change. The key objective of this study was to simulate the streamflow and to predict the water levels of the Mackenzie River for the 2050s based on the downscaled data of an RCM (WRF model) driven by boundary conditions of a selected GCM. The outcome of this study will provide some quantitative estimate of possible hydrologic changes, especially the water levels of the Mackenzie River in future.

MRB is a high latitude continental basin with a total area of about 1.8 million square kilometers, extending from 52° to 69°N and 140° to 102°W, making it the largest river basin in Canada. The western part of MRB is mountainous and has an average elevation of about 1,000 m, whereas the interior part is a plain and the eastern part is the Canadian Shield. The map of the MRB with terrain height is shown in Figure 1. The northern part of the basin experiences harsh winter, with subzero temperature for the whole winter, with an average winter temperature ranging from −25 to −35 °C, and −50 °C is not uncommon; whereas, in the summer, the average monthly temperature ranges from 8 °C in the north to 20 °C in the south, which could be as high as 35 °C. So, the inter-annual temperature variation is substantial within this region. From 1979 to 1998, the average annual precipitation over the MRB was about 410 mm. The annual mean evapotranspiration was about 230 mm, which mainly occurs between May and October.

The Mackenzie River represents a critical corridor in Canada's Arctic transportation network, and it plays a key role in the transportation of northern communities of the Northwest Territories (NWT). In an open water season, ferry transport provides the essential transportation across the MRB. The peak daily flow of the Mackenzie River at the downstream station of the Arctic Red usually varies between 15,800 and 35,000 m³/s, lowest and highest, respectively (Yang et al. 2015) during the spring–summer season; whereas the annual mean discharge of the Mackenzie River into the Arctic Ocean is about 9,100 m³/s. From May to October, daily peak flows occur in the spring to early summer due to the snowmelt breakup of river ice and summer flows are also high. Trend analyses reveal a statistically significant (85% confident) tendency of decreasing peak flow by about 3,000 m³/s over the 1973–2011 period (Yang et al. 2015). Time series and trends analysis of the flow ratio for 1973–2011 is shown in Supplementary Appendix 1.

During the relatively short open water season, the Mackenzie River carries barge traffic for the community, including bulk fuel shipments (Prolog Canada 2011). The Mackenzie River is also a highly active transportation corridor for residents. Most people living along the river travel by snowmobile on the ice cover in winter and by personal watercraft in the summer. In other words, the water levels of the Mackenzie River play a major role in the transportation of the northern communities of the NWT in the summer.

A conceptual hydrological model called the HBV model (Bergström & Forsman 1973) was chosen to simulate the streamflow of the Mackenzie River at the selected stations. The HBV is a conceptual semi-distributed, rainfall-runoff model that has been widely used for hydrological simulations and has been applied in numerous studies for the computation of design floods, flood forecasting, and also for climate change studies (Driessen et al. 2010; Al-Safi & Sarukkalige 2017). The HBV model has been extensively used in operational hydrological forecasting and water balance studies (Abebe et al. 2010). This model typically operates on daily time steps. Input data to HBV include precipitation, temperature, and potential evapotranspiration (PET).

In this study, climate data such as precipitation and temperature were calculated as a weighted mean of climate stations in and around the basin. A Thiessen polygon method was used to determine the relative weight of different stations for mean daily temperature and precipitation. When gridded climate data are used, the areal average is obtained by this method. A similar method (Thiessen polygon) was utilized by Kobold & Brilly (2006) to prepare the areal average of the climatological data to input to the HBV model.

Hamon's (1961) method is one of the best-suited methods for calculating the PET suggested by Jianbiao et al. (2005), as he compared six different methods for calculating the PET. In this study, we also used Hamon's (1961) temperature-index PET equation to transform air temperature to potential evaporation as one of the model inputs (Equation 1 in Supplementary Appendix 2). The snow module of HBV uses a degree-day method to estimate snow accumulation and snowmelt processes. Precipitation is only considered as rainfall when the air temperature is above a threshold temperature TT (°C); otherwise, it is treated as snowfall. The amount of snowmelt, M (mm day⁻¹), and the amount of refreezing liquid water within the snowpack, R (mm day⁻¹), are computed from Equations 2 and 3, respectively (Supplementary Appendix 2). The soil moisture module of the HBV is based on three empirical parameters: β, FC, and LP, where β controls the contribution of input (I(t)) in mm, which consists of precipitation and snowmelt, the runoff response (Q_S(t)), and the contribution of (I(t) − Q_S(t)) to the soil moisture storage (S_sm(t)). Simulations of the HBV model are very sensitive to changes of CFMAX or FC when they have changed alone (Seibert 1997). Soil moisture and evapotranspiration components FC, β, and to a certain extent, LP are sensitive, whereas FC was found to be the most dominant parameter of the model, which affects both the high-flow series and the volume errors, showing sensitivity to both RMSE and BIAS objective measures. Melching et al. (1990) stated that the model structure uncertainty and parameter uncertainty may be a significant source of the combined modeling uncertainty.

Excessive water from precipitation and snowmelt is transformed by the runoff response function to Q_S(t) (Equation 4 in Supplementary Appendix 2). Evapotranspiration E_a is calculated using Equation 5 (Supplementary Appendix 2). The runoff response function and the model structure of the HBV are shown in Supplementary Appendix 3. The model consists of two conceptual storages (tanks) that simultaneously re-distribute the generated runoff in terms of quick and slow responses, respectively. Q_u(t) is obtained from Equation 6, and then, Q_l(t) and Q_sim(t) are calculated using Equations 7 and 8, respectively (Supplementary Appendix 2).

The Monte Carlo procedure entails performing a large number of simulations, and each individual simulation was allocated random parameter values from pre-defined parameter ranges. The other two HBV model parameters, refreezing coefficient (CFR) and water holding coefficient (CWH), are set at 0.05 and 0.1, respectively, as recommended, instead of being calibrated. Optimal values are fine-tuned using an objective function based on the Nash–Sutcliffe (NS) coefficient (Nash & Sutcliffe 1970; Equation 9 in Supplementary Appendix 2) that essentially minimizes differences between the simulated runoff of the HBV with observed streamflow. The calibrated HBV parameters are then validated using data independent of the calibration period.

The Australian National University SPLINe (ANUSPLIN) data are daily observational climate data produced by Natural Resources Canada and were available at 300 arc-second (10 km) spatial resolution over Canada from 1950 to 2015. The ANUSPLIN likely is the best-gridded dataset available for Canada and has been used as the source data to compare climate products (Eum et al. 2014; Wong et al. 2017) and to evaluate the accuracy of regional climate models (Eum et al. 2012) for Canada.

Pacific Climate Impacts Consortium (PCIC) offers statistically downscaled daily Canada-wide climate scenarios, at a gridded resolution of 300 arc-seconds for the simulated period of 1950–2100. The PCIC downscaled datasets are based on GCMs' projections using the BCSD statistical downscaling method (Bias Correction with Constructed Analogues and Quantile mapping, Version 2; BCCAQv2) derived from 24 CMIP5 GCMs for two emissions scenarios (RCP 4.5 and RCP 8.5). For each scenario, the PCIC dataset provides daily temperatures (maximum and minimum) and total precipitation at a 10 km × 10 km resolution for all of Canada. Climate projections are provided for two future 30-year periods (2021–2050 and 2051–2080) and a baseline period (1976–2005). Studies showed that the high-resolution gridded data from PCIC are suitable for climate change research (Islam et al. 2017). Data are available at https://www.pacificclimate.org/data/statistically-downscaled-climate-scenarios.

Water Survey Canada's Hydrometric Database is HYDAT. These data are collected and compiled by the Water Survey of Canada's eight regional offices. HYDAT is a relational database that contains the actual computed data for the listed stations. These data include daily and monthly means of flow, water levels, and sediment concentrations (for sediment sites). For some sites, peaks and extremes are also recorded. HYDAT have been extensively used for streamflow and climate change impact studies (Burn et al. 2010; Tan & Gan 2015). In this study, daily flow data were collected from 1974 to 2004 at the Fort Simpson station and the Arctic Red River station.

Tables relating to the stage and discharge data at several gauging stations for the Mackenzie River were obtained from the Water Survey of Canada, Environment Canada. Rating curves were developed from the stage–discharge data that will be applied to convert simulated discharge values to river stages value in meters. For the Fort Simpson station of the Mackenzie River, a table was obtained for a river stage ranging from 1 to 10.28 m at 0.01 m increments. At the Arctic Red River station, the corresponding table of a stage ranging from 1.3 to 9.50 m at 0.01 m increments was obtained. Supplementary Appendix 4 illustrates the stage–discharge relationships at the Fort Simpson and the Arctic Red River Stations.

The second-generation Canadian Earth System Model (CanESM2) consists of the physical coupled atmosphere-ocean model CanCM4 coupled to a terrestrial carbon model (CTEM) and an ocean carbon model (CMOC). The raw data collected from the Canadian GCM (CanESM2) model has spatial resolutions 2.81° by 2.81° (310 km × 310 km). The CanESM2 is a widely used GCM, and it provided reliable results for many climate change studies (Gebrechorkos et al. 2019). To evaluate the performance of the Canadian GCM to simulate the local hydrology, the CanESM2 model data were chosen. A Regional Climate Model is necessary to downscale the GCM model data into a regional scale before using it in any hydrological models. Figure 2 illustrates the typical GCM, Regional Climate Model (RCM), and local hydrology parameters.

The calibrated and validated WRF model setup (Pervin & Gan 2020) was used to simulate temperature, rainfall, and other climate data for May, June, July, August, September, and October (MJJASO) over the MRB for the historical period (1979 to 2005), and the WRF simulated historical data were compared with the reference ANUSPLIN data. Figures 3 and 4 show the WRF simulated temperature and rainfall data (respectively) compared with the reference data.

For downscaling the future data, the CanESM2 Representative Concentration Pathway, RCP 4.5 and RCP 8.5 climate scenarios (IPCC 2013) were chosen which are called the ‘Low Carbon’ and ‘High Carbon’ scenarios, respectively. The WRF model then ran to downscale the temperature and rainfall data for the future period 2041–2070 for both RCP 4.5 and RCP 8.5 climate scenarios. The output temperature and rainfall data were 30 km × 30 km in spatial resolutions at 6 hourly time steps. Figures 5 and 6 illustrate the future data and the anomaly with respect to the historical data for the RCP 4.5 climate scenario, and Figures 7 and 8 represent the same but for the RCP 8.5 climate scenario. The downscaled summer (MJJASO) temperature and rainfall for the 2050s (2041–2070) under RCP 4.5 and RCP 8.5 were then converted to daily time steps, and PET data were calculated. These future temperature, rainfall, and PET data will be used as the input to the HBV hydrological model to simulate the streamflow for the future period at the selected stations of the Mackenzie River to assess the future water levels.

HBV was calibrated against streamflow of seven river gauging sites separately, starting first from sub-basin #1, then sub-basin #2, until sub-basin #7, based on simulated streamflow versus observed streamflow for each individual sub-basin. This resulted in seven sets of optimized HBV parameters for the seven sub-basins of the MRB. Streamflow of the Peace River has been regulated after the W.A.C. Bennet Dam was built. Naturalized streamflow for stations 07HA001 and 07KC001 were obtained from the Government of Alberta. An upstream portion of the Athabasca River basin was modeled at station 07BE001. Station 07NB001 on the Slave River included the flow from the Athabasca River and the Peace River. The Liard river streamflow was modeled at station 10ED002 which, in addition to the outflow from the Great Slave Lake, contributed to the streamflow at the Fort Simpson at station 10GC001. The Fort Simpson station was the most upstream site where the observed streamflow of the Mackenzie River was used for calibration. The HBV model is not suitable to model the hydrology of the Mackenzie Delta because it is a wetland. Therefore, the most downstream station selected for the Mackenzie River was the Arctic Red River at station 10LC014. Based on the availability of observed input data (streamflow, precipitation, and temperature), the calibration period was selected as 1 January 1974 to 31 October 1997 and the validation period from 1 April 1998 to 31 December 2004. The first year (1 January 1973 to 31 December 1973) data were used as spin-up in the HBV and, therefore, were not included in the analysis. The statistics and NS coefficient for both the calibration and validation period are illustrated in Table 1.

The HBV model was automatically calibrated using streamflow data from 1 January 1974 to 31 October 1997 for each of the seven selected sub-basins separately. Göran (1997) introduced a simple automatic calibration routine for the HBV model. The optimization is made for one parameter at a time, while the others are kept constant. This one-dimensional optimization is repeated in a loop for all parameters. Results for the Fort Simpson station 10GC001 (61.9°N, 121.4°W) and the Arctic Red River station 10LC014 (67.5°N, 133.7°W) are discussed below, since these two stations are the main interest of the study. Both the climate station data and the ANUSPLIN data were used in the calibration, which yielded different optimized parameters as the ANUSPLIN is gridded data and not adjusted for changes in elevation. The coefficient of determination and the NS coefficient for the Fort Simpson station based on the climate station data are 0.86 and 0.82 for the calibration period, respectively. The corresponding statistics using the ANUSPLIN as the calibration data were marginally better, 0.89 and 0.88, respectively. The R² and NS for the validation runs at the Fort Simpson using the climate station (ANUSPLIN) data were 0.85 (0.85) and 0.77 (0.81), respectively. At the Arctic Red River station, the climate station data had goodness-of-fit statistics ranging between 0.84 and 0.87 for the calibration and validation runs. The ANUSPLIN data resulted in an R² of 0.89 and NS of 0.89 for the calibration run, and an R² of 0.87 and NS of 0.83 in the validation run, respectively.

Given that HBV's simulated streamflow driven by both the climate station and the ANUSPLIN input data generally matches well with the observed flow at the validation stage, which means that the HBV is well-calibrated, its optimized model parameters are physically meaningful, and therefore, it can be used to simulate the future streamflow of the Mackenzie River. The HBV accurately simulated the spring snowmelt which has a sharp peak compared with the winter flow, but it marginally under-simulated the falling limb of the hydrograph. The simulated base flow in the winter is under-stimulated compared with the observed flow. This is expected as the observed flow includes the effect of regulated streamflow which have higher winter releases. The HBV model was used to simulate historical water levels using the WRF downscaled data and the PCIC data, then the modeled water levels were compared with the observed water levels for both stations. The streamflow simulation using the WRF downscaled CanESM2 historical data are comparable with the observed streamflow, for both the Fort Simpson and the Arctic Red River stations, the R² values were found to be 0.85 and 0.86, respectively. Figure 9(a) and 9(b) is presented at a weekly time-scale averaged over the historical period in terms of water levels of the Mackenzie River.

After the HBV model was calibrated and validated with the HYDAT streamflow data of Environment Canada, WRF outputs (rainfall and temperature) were forced with the HBV to simulate the future streamflow of the MRB. Since the WRF over-simulates the precipitation and slightly under-simulates the temperature over the MRB, the quantile–quantile bias correction method was applied to the WRF downscaled CanESM2 data with respect to the ANUSPLIN data, before being used to drive the HBV model. The streamflow hydrographs of MRB simulated using the WRF downscaled CanESM2 historical data were compared with the observed historical streamflow at the Fort Simpson station and the Arctic Red River Station; they were further compared with the simulated streamflow using the PCIC data for the historical period at the same locations.

Then, the future water levels at the selected stations were simulated using the HBV model; the future climate scenarios RCP 4.5 and RCP 8.5 of CanESM2 data downscaled by the WRF were input to derive the streamflow hydrographs of MRB for the 2050s. Figure 10(a) and 10(b) shows the future water levels in the 2050s of CanESM2 RCP 4.5 and PCIC 2050s at the Fort Simpson and the Arctic Red River stations; in comparison with the observed historical water levels, it is projected that the water levels for the 2050s over the MJJASO period will be significantly decreasing at the Fort Simpson station; whereas the Arctic Red River will be experiencing reduced peak flow by the 2050s.

Seasonal (MJJASO) streamflow hydrographs of the MRB simulated for the 2050s expressed in terms of water levels of the Mackenzie River for RCP 8.5 are presented in Figure 11(a) and 11(b) at weekly time-scale averaged over 30-year periods for the Fort Simpson and the Arctic Red River stations, respectively. The corresponding observed water levels over the 1974–2004 period are also plotted to demonstrate the projected changes to the water levels of the Mackenzie at the Fort Simpson River and the Arctic Red River stations in the 2050s for RCP 8.5 climate scenarios. Substantial lowering of water levels at the Fort Simpson station is projected by the 2050s using the WRF downscaled CanESM2 RCP 8.5 data; whereas the Arctic Red River will be experiencing reduced peak flow which is also projected to shift about 2 weeks earlier than the historical peak flow.

Since the Fort Simpson station shows a significant change in the future water levels, we use the boxplots to statistically compare the results. In Figure 12, the boxplots represent the average weekly water levels of MJJASO for each time period (historical and 2050s for RCP 4.5 and RCP 8.5). The observed historical mean water level at the Fort Simpson station matches with the WRF downscaled CanESM2 data-derived HBV simulated mean historical water levels. In the 2050s, the mean water level will be lowered by 300 cm with respect to the historical period; on the other hand, the peak flow and the high water events will be lowered significantly for both RCP 4.5 and RCP 8.5 scenarios.

At the Fort Simpson Station, mean flow is projected to decrease by about 5% under both RCP 4.5 and RCP 8.5 scenarios in the 2050s, whereas the peak flow is likely to reduce by about 12 and 9% for RCP 4.5 and RCP 8.5 scenarios, respectively, compared with the historical streamflow. However, the projected minimum flow indicates that it will likely increase a little in the 2050s.

Under a warmer climate, more evaporation loss is expected in the summer; this increase of evaporation could offset the projected increase in precipitation, giving rise to an overall decrease in streamflow and lower water levels in the 2050s. The projected lower streamflow and water levels could significantly affect the future hydrology and thus the navigability of the Mackenzie River which is the backbone of water transportation in the western portion of the NWT. Navigation problems of the Mackenzie River include a short shipping season (beginning of June to mid-October) and lower water levels. Given that the summer water levels of the Mackenzie are projected to decrease in the 2050s, navigation problems related to low water levels are expected to increase in future because safe transit through the Mackenzie River depends on its water levels.

To project the future streamflow and to assess the changes in water levels of the Mackenzie River, a conceptual hydrological model (HBV) was used in this study. In order to simulate realistic streamflow, the Canadian GCM CanESM2 model data were downscaled to 30 km × 30 km resolution using a regional climate model called the WRF model. Then, the WRF downscaled climate data were utilized as input data to the HBV model. Two river stations, namely the Fort Simpson and the Arctic Red River station, were selected, which are the most upstream and downstream stations of the Mackenzie River. Based on the availability of observed input data (streamflow, precipitation, and temperature), the calibration period was selected as 1 January 1974 to 31 October 1997 and the validation period from 1 April 1998 to 31 December 2004. The R² and NS for the validation runs at the Fort Simpson station using the climate station data were 0.85 and 0.77, respectively. At the Arctic Red River station, the climate station data had goodness-of-fit statistics ranging between 0.84 and 0.87 for the calibration and validation runs. As the input data to simulate historical streamflow, the downscaled GCM CanESM2 data by the WRF model, which are 30 km × 30 km in spatial resolution, were forced with the calibrated and validated HBV model. The simulated historical flow with the WRF data was compared with the measured station data and also with the simulated streamflow with the PCIC data; in both cases, the results were comparable. For the future streamflow simulations at the selected stations, the WRF downscaled CanESM2 data for the 2041–2070 period under both future climate scenarios RCP 4.5 and RCP 8.5 were input to the HBV model. The WRF model outputs were bias-corrected with respect to the ANUSPLIN data. Again, we also applied the PCIC data for the same period (2041–2070) to the HBV model to simulate the future streamflow in order to compare the results using the WRF model outputs. The simulated discharges were converted to stage data using the rating curve of that station. To assess the changes in future water levels of both locations, we compared the future water level data with the historical water level data. The observations are as follows:

• (a)

  Under the projected increase in air temperature and precipitation of RCP 4.5 and RCP 8.5 climate scenarios of the CanESM2 downscaled by the WRF, the streamflow of the Mackenzie at the Fort Simpson and the Arctic Red River stations shows a general decrease of flow in future. Under a warmer future climate, more evaporation loss is expected in the summer which could offset the projected increase in summer precipitation, resulting in an overall decrease in water levels in the 2050s.

• (b)

  The average annual maximum water levels are also expected to decrease significantly under both scenarios of future periods.

• (c)

  The projected lower streamflow and water levels could significantly affect the future hydrology and thus the navigability of the Mackenzie River which is the backbone of water transportation in the western portion of the NWT. Navigation problems of the Mackenzie River include a short shipping season (beginning of June to mid-October) and lower water levels. Given that the summer water levels of the Mackenzie are projected to decrease in the 2050s, problems related to low water levels are expected to increase in future.

We thank Transport Canada for funding this study; Environment Canada, NCEP, and ECMWF for the use of their data source; WestGrid supercomputing program and the wrf-model.org.

The authors declare no potential conflict of interests.

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