Long Term Trend Analysis of River Flow and Climate in Northern Canada

: Changes in water resources within basins can signiﬁcantly impact ecosystems, agriculture, and biodiversity, among others. Basins in northern Canada have a cold climate, and the recent changes in climate can have a profound impact on water resources in these basins. Therefore, it is crucial to study long term trends in water ﬂow as well as their inﬂuential factors, such as temperature and precipitation. This study focused on analyzing long term trends in water ﬂow across the Athabasca River Basin (ARB) and Peace River Basin (PRB). Long term trends in temperature and precipitation within these basins were also studied. Water ﬂow data from 18 hydrometric stations provided by Water Survey of Canada were analyzed using the Mann-Kendall test and Sen’s slope. In addition, hybrid climate data provided by Alberta Environment and Parks at approximately 10 km spatial resolution were analyzed for the ARB and its surrounding regions during 1950–2019. Trend analysis was performed on the water ﬂow data on monthly, seasonal, and annual scales, and the results were cross-checked with trends in temperature and precipitation and land use and land cover data. The overall temperature across the basins has been increasing since 1950, while precipitation showed an insigniﬁcant decrease during this period. Winter water ﬂow in the upper ARB has been slowly and steadily increasing since 1956 because of the rising temperatures and the subsequent slow melting of snowpacks/glaciers. The warm season ﬂows in the middle and lower subregions declined up to 1981, then started to show an increasing trend. The middle and lower ARB exhibited a rapid increase in warm-season water ﬂow since 2015. A similar trend change was also observed in the PRB. The gradual increase in water ﬂow observed in the recent decades may continue by the mid-century, which is beneﬁcial for agriculture, forestry, ﬁshery, and industry. However, climate and land cover changes may alter the trend of water ﬂow in the future; therefore, it is important to have a proper management plan for water usage in the next decades.


Introduction
Climate dynamics can significantly affect hydrological regimes in river basins across the world. Changes in water availability in watersheds have profound impacts on ecosystems, agriculture, forestry, fishery, industry, and hazard management [1]. Therefore, continuous monitoring of water flow and the monitoring of its influential factors, such as precipitation and temperature, are crucial. Hydrological modelling for large basins is very difficult due to hydroclimatic complexity linked to changes in land cover, topography, soil drainage capacity, climate, and anthropogenic activities. In cold weather basins, water flow during winter is usually low, while during spring and summer, precipitation and snowmelt significantly increase water flow, which may result in flooding [2]. Watersheds located in high latitudes in Canada are subject to rapid hydroclimatic changes that directly influence the water availability in the regions. In the Athabasca River Basin (ARB), the changing hydrology impacts the natural ecosystem in national parks and reserves, the livelihood of the population, the agricultural and mining industries, as well as the oil and gas industry in the lower ARB [1]. The Peace-Athabasca Delta (PAD), which falls at the confluence of the Peace and Athabasca Rivers with Lake Athabasca, is home to a diverse ecosystem and is a breeding environment for wildlife [2]. Fossil fuel excavation and extraction within the ARB has utilized a growing amount of the Athabasca River, about 4.5% of river's annual flow in 2012 [3].
Many researchers have studied trends in water flow across the ARB. A general decline in the summer water flow has been reported across Canadian rivers due to climatic and atmospheric factors [4,5]. The ARB has demonstrated an earlier spring peak water flow, a higher winter flow [6][7][8], and around a 20% reduction in the summer flow [9]. These observations were supported by studying the trends in the runoff, temperature, and precipitation along the ARB and their impact on water flow [10]. Historical data have also been used to analyze and explain different hydrological phenomena using modelling [11]. Hwang et al. (2018) [12] studied the water balance in the ARB to explain the increasing downstream flow during water deficit conditions using a 3D model for the surface and subsurface flow and evapotranspiration. They also integrated land cover data to investigate the role of different types of vegetation in the water balance.
Cold weather basins are generally very vulnerable to climate change. For example, the highly glacierized and mountainous Upper Indus Basin was shown to be highly affected by climate change [13]. The findings of Shah et al. [13] imply that other similar basins may also be affected by climate change, such as the ARB, since it is also mountainous and highly glacierized in the upper ARB and extends to prairies in the lower ARB. In addition to climate, land use and land cover (LULC) also contribute to the hydrological behavior of watersheds. With the changing climate, new analyses of the ARB hydrology with newly available data are required to have an up-to-date understanding of the current hydrological conditions and to make informed projections and decision making. Remote sensing data for 2001-2020 show a warming trend in the Canadian province of Alberta, where the ARB lies [14]. A recent LULC map was developed for the ARB and was used to detect changes in vegetation and human activities within the watershed, which can have a significant impact on the Athabasca River water flow [12,15]. A combined analysis of water flow, temperature, precipitation, and land cover can help in proper water resource planning, both currently and in the future.
In this study, hydrometric and hybrid climate datasets were analyzed to obtain a current hydrological understanding of the ARB. To the authors' best knowledge, this is the first time that such analyses have been performed spatially for over 18 hydrometric stations, at pixel and subregion levels for climate, and temporally at monthly, seasonal, and annual scales during 1950-2020. These findings are comprehensive and provide more in-depth understanding of climate and water flow in the ARB and its north-western basin, the Peace River Basin (PRB). Bawden et al. [16] applied the Mann-Kendall (MK) test and Sen's slope to hydrometric time series data for the ARB and surrounding regions for three different periods: 1976-2010 (35 years), 1971-2010 (40 years), and 1966-2010 (45 years). They investigated the trends of monthly water flow, annual mean and peaks, and warm season mean between 1976 and 2010. They compared the trends in the ARB with the surrounding watersheds using data from the Water Survey of Canada (WSC) and Reference Hydrometric Basin Network (RHBN) gauging stations. They also analyzed the trends in temperature and precipitation in the ARB using a 10 km resolution dataset [17]. In the present study, on one hand, eight different periods were investigated at monthly, seasonal, and annual scales: 1991-2020 (30 years), 1986-2020 (35 years), 1981-2020 (40 years), 1976-2020 (45 years), 1971-2020 (50 years), 1966-2020 (55 years), 1961-2020 (60 years), and 1956-2020 (65 years). From the statistical point of view, estimated trends for longer periods have greater power than shorter periods. However, more stations could be analyzed for recent years because many stations have started to operate recently, and some had major data gaps in the past. Furthermore, the potential regime shifts in water flow were investigated through a windowing strategy. On the other hand, the hybrid climate data, obtained from multiple climate datasets and prepared for the ARB and its surrounding regions, were recently released for the period 1950-2019, providing continuous and reliable data since 1950, which also cover the mountainous regions [18]. The long historical and continuous hybrid climate data allow the rigorous investigation of climate impact on water flow both in time and space, which were utilized in the present study. Moreover, the possible influence of recent land cover changes (2001-2020) on waterflow dynamics is also discussed here.
Note that the detailed water flow analysis results are presented for three main stations in the ARB located at three subregions of the ARB with long historical and continuous data, and the results for the rest of the stations are presented in the Supplementary Materials. The main contributions of this work are summarized below.
(1) Long-term and short-term trend analyses of hydrometric datasets from 18 gauging stations in ARB and PRB at monthly, seasonal, and annual scales. These analyses show the direction of historical and recent water flow trends as well as the detection of months and seasons when the flow was decreasing or increasing. (2) Estimated monthly hydrographs for the water flow datasets for 30-year-long moving windows since 1956 to find whether there is any regime shifts in the mean of water flow for each station. (3) Estimated trends for temperature and precipitation at both pixel and subregion levels to investigate their possible relations with waterflow dynamics. (4) Discussed the recent trends in land cover types across the subregions and their potential interactions with water flow dynamics.

Study Region
The Athabasca River Basin (ARB) is one of the largest watersheds in Alberta, with a drainage area of approximately 138,000 km 2 . The Athabasca River extends from the Athabasca Glacier, which extends from the Columbia Glacier in the Rocky Mountains and flows North-East into Lake Athabasca, then into the Arctic Ocean through the Mackenzie River [17]. The topography and ecology vary along the ARB, where the upper reach of the river flows through the Rocky Mountains with alpine, sub-alpine, and montane ecoregions. The middle reach of the Athabasca River then flows through low-populated foothills and prairie lands with boreal forests, where the river is surrounded by forests, coal mines, limestone quarries, and croplands. The McLeod River, Pembina River, and Lesser Slave River flow into the Athabasca River within the middle reach. The lower reach of the river starts at the urban centre of Fort McMurray. The lower ARB is the national centre for the oil industry, the third-largest oil deposit in the world. The Clearwater River and other small tributaries flow into the Athabasca River within the lower reach. Fort McMurray is highly populated (~60,500 in 2011 and~66,500 in 2016) relative to other towns in the ARB, with an expanding residential and industrial footprint due to the influx of workers in the energy sector [19].
For the ARB, the annual mean daily maximum temperature was 6.2 • C, and the annual mean daily minimum temperature was −5.2 • C on average during the period 1950-2019. Additionally, the average annual precipitation was 492 mm for the ARB. Around 75% of the precipitation in the ARB is in the form of rain, while the rest is snow [16]. Water flow in the Athabasca River is characterized by low winter flows, followed by a rapid increase due to snowmelt in April. Peak flow occurs between May and August, depending on the location of the gauging station along the river. Figure 1 shows the location and topography of the ARB within the provinces of Alberta and Saskatchewan, as well as the locations and ID numbers of the gauging stations in and around the ARB that were used in this study. The details of the gauging stations are shown in Table 1. According to multiple criteria, such as topography, land cover, climate, and biodiversity, the ARB is divided into three subregions, namely the upper (the southern mountainous region), middle, and lower ARB [20], separated by dashed lines in Figure 1 (Hatfield Consultants 2008, see the Regional Aquatics Monitoring Program (RAMP): http://www.ramp-alberta.org/river/ geography/basin+landscape.aspx, accessed on 2 April 2022). The focus of this study was on the ARB; however, for comparison purposes, three stations with long and continuous historical discharge records were also chosen in the PRB (see p1, p2, and p3 in Figure 1). topography of the ARB within the provinces of Alberta and Saskatchewan, as well as the locations and ID numbers of the gauging stations in and around the ARB that were used in this study. The details of the gauging stations are shown in Table 1. According to multiple criteria, such as topography, land cover, climate, and biodiversity, the ARB is divided into three subregions, namely the upper (the southern mountainous region), middle, and lower ARB [20], separated by dashed lines in Figure 1 (Hatfield Consultants 2008, see the Regional Aquatics Monitoring Program (RAMP): http://www.ramp-alberta.org/river/geography/basin+landscape.aspx, accessed on 2 April 2022). The focus of this study was on the ARB; however, for comparison purposes, three stations with long and continuous historical discharge records were also chosen in the PRB (see p1, p2, and p3 in Figure 1).

Figure 1.
Map of the ARB and its location and topography, and the ID and locations of the waterflow gauging stations used in this study. Hydrometric stations marked with closed circles were considered to represent the subregion in which they lie because they are located along the Athabasca River and provided long-duration and continuous flow measurements. Note that l2 was the only station that provided continuous historical records since 1956 as compared to other stations in the lower ARB; see the Supplementary Materials. Stations p1, p2, and p3 are in the PRB. Table 1. Gauging station attributes in the ARB and surrounding area. Letters u, m, and l in the first column refer to the upper, middle, and lower ARB, respectively. Additionally, letter p in the first column refers to the PRB.   A high-resolution (0.1 × 0.1 degree or~10 km within ARB) hybrid gridded climate dataset for 1950-2019 was generated from a framework called the REFerence Reliability Evaluation System (REFRES) that systematically generates a performance-based hybrid climate dataset from multiple climate datasets for a region. In fact, Eum and Gupta [18] recommended the use of the hybrid climate dataset as reference for hydrological model forcing and statistical downscaling because it is a good proxy for the historical climate conditions. REFRES is composed of three modules, performance measures, ranking, and data generation.

Water Flow Data
The water flow data used in this study were acquired from the Water Survey of Canada (WSC). The WSC collects real-time hydrometric data (water flow and water level) from gauging stations across Canada. The parameters used in this study, which were extracted from the WSC data, were the monthly averages, annual averages, cold season averages, and open warm-season averages. The months included in the calculation of the cold and warm seasons varied across the ARB based on the observed flow conditions. In the upper ARB, the cold season with relatively higher water flow included the months November-April, while the warm season was May-October. The Middle and Lower cold and warm seasons were November-March and April-October, respectively. The WSC gauging stations were set up at various times since 1914, leading to a varying range of recorded data. Some stations also had gaps in recording that ranged between months to years, and other stations did not record in the winter months (seasonal operation). The gauging stations selected in this study were chosen based on data availability. Chosen stations were those that had continuous data, where the gaps in any parameter were not more than five years in total or three consecutive years in any gauging station, according to the 1981-2010 Canadian Climate Normals Report (https://climate.weather.gc.ca/climate_normals/, accessed on 2 April 2022). Furthermore, the chosen stations recorded data over the 12 months without seasonal discontinuity in the winter season, and finally, the chosen stations had a period of record that extended at least 30 years prior to 2020.
In addition to the WSC water flow data, another dataset was included in the study that covered the same region, even though it has a shorter period with more missing values than the WSC. The Regional Aquatics Monitoring Program (RAMP) was an industry funded program between 1997 and 2015. Its intent was to monitor the aquatic environments in the Wood Buffalo Area and study the impact of the oil sands projects. The RAMP program collected water flow data across the lower ARB and complemented the existing WSC data.

Land Use and Land Cover (LULC) Data
The Canada Centre for Remote Sensing (CCRS) has generated a land cover map of Canada for the base year 2010, as the first of a planned series of maps to be updated every five years or more frequently. This land cover dataset is also the Canadian contribution to the 30 m spatial resolution 2010 Land Cover Map of North America, produced by Mexican, American, and Canadian government institutions under a collaboration called the North American Land Change Monitoring System (NALCMS) [21]. The NALCMS is a collaborative framework involving CCRS' parent organization, Natural Resources Canada (NRCan), the United States Geological Survey (USGS), and Mexican institutions, including the National Institute of Statistics and Geography (INEGI), the National Forestry Commission (CONAFOR), and the National Commission for Knowledge and Use of Biodiversity (CONABIO). The 2015 Land Use and Land Cover (LULC) classes were defined using the Land Cover Classification System (LCCS) standard developed by the Food and Agriculture Organization (FAO) of the United Nations. This map is based on 2015 Landsat satellite imagery for Canada at 30 m spatial resolution (the map is publicly available through the link: (http://www.cec.org/north-american-environmental-atlas/land-cover-30m-2015-landsat-and-rapideye/), accessed on 23 November 2021).

Climate Trend Analysis Methods
Linear regression models the relationship between dependent variables (e.g., temperatures) and independent variables (e.g., times). The least-squares approach is usually employed in linear regression to estimate the slope (m) and intercept of the linear fit line through minimizing the residual series. Considering the normality and independence of the residual series, one can determine the statistical significance of the estimated slope via the test statistics t =m/E, which follows the student's t-distribution, where E is the standard error of the estimated slopem. In testing the statistical significance of the estimated slope, the null hypothesis is usually defined as no trend (m = 0), and the alternative hypothesis is defined as trend exists (m = 0).

Water Flow Trend Analysis Methods
The Mann-Kendall (MK) method is a non-parametric test that is common in detecting trends in hydrometric data [22,23]. It is based on the null hypothesis that there is no trend in the data, i.e., the measurements acquired over time are independent and identically distributed. The MK test can detect trends in data at certain significance levels (90%, 95%, and 99% in this study). The trend slopes were calculated using the Sen approach [24]. Sen slope is a non-parametric slope estimator that determines the magnitude and direction of the trend in hydrometric variables. The mathematical formulation of the MK test and Sen slope is outlined below. The MK test statistic S is calculated by: and where n is the number of measurements, and x j and x k are data values in time series i and j. The variance is calculated by: when n >10, the standard normal test statistic Z is calculated by: Hydrology 2022, 9,197 7 of 21 Positive Z values indicate an increasing trend and vice versa. The Sen slope β is calculated by:

Climate Trends
Climate is one of the main drivers of vegetation coverage and water flow. However, there is usually a time lag between these parameters, and their impact on each other depends on several other factors, such as urbanization, harvesting activities, soil types, topography, etc. To show the spatial variation of climate over the ARB and its surroundings, the mean of daily maximum and minimum temperatures was calculated between 1950 and 2019, as well as the mean of yearly-accumulated precipitation within that period. The results are shown in Figure 2. Furthermore, the overall rates of change in temperature and precipitation were calculated during that period, and the climate change map is also illustrated in Figure 2. Temperature has been increasing over the region, particularly in the middle and lower ARB. Additionally, precipitation has been decreasing within the ARB since 1950, particularly in the upper ARB.
Positive Z values indicate an increasing trend and vice versa. The Sen slope is calculated by:

Climate Trends
Climate is one of the main drivers of vegetation coverage and water flow. However, there is usually a time lag between these parameters, and their impact on each other depends on several other factors, such as urbanization, harvesting activities, soil types, topography, etc. To show the spatial variation of climate over the ARB and its surroundings, the mean of daily maximum and minimum temperatures was calculated between 1950 and 2019, as well as the mean of yearly-accumulated precipitation within that period. The results are shown in Figure 2. Furthermore, the overall rates of change in temperature and precipitation were calculated during that period, and the climate change map is also illustrated in Figure 2. Temperature has been increasing over the region, particularly in the middle and lower ARB. Additionally, precipitation has been decreasing within the ARB since 1950, particularly in the upper ARB. The MK trend analysis and Sen's slope of the climate data also demonstrated a decline in the precipitation across the ARB while temperatures were increasing, as shown in Table 2. Precipitation in the upper ARB has been decreasing by 11.28 mm per decade at a 90% significance level, since 1950. The precipitation in the middle and lower ARB has also been declining but with a significance level of less than 90%. As it can be seen in the bottom right panel of Figure 2, the increase and decrease in the precipitation rate were approximately equally distributed within the middle and lower ARB, while the upper ARB showed a more significant decline in the precipitation rate. The temperature, however, has been increasing across the ARB with over 99% significance. The maximum temperature has been uniformly increasing in the upper, middle, and lower ARB with rates of 0.23, 0.27, and 0.28 • C per decade. The minimum temperature demonstrated an increasing tapered rate along the ARB, where the upper, middle, and lower subregions increased by 0.17, 0.33, and 0.45 • C per decade.   Table 2, further confirming the trend results presented here. The climate time series for the lower, middle, and upper ARB are shown in Figure 3, Figure 4, and Figure 5, respectively. The average temperature was simply obtained by averaging the maximum and minimum temperatures. Note that "Annual Max Temp" on the y-axis is the mean of daily maximum temperatures in each year. Likewise, "Annual Min Temp" is the mean of daily minimum temperature in each year. To compare the past and current overall changes in temperature and precipitation, the time series were segmented from 1950 to 1990 (1950-1989) and from 1990 to 2020 (1990-2019). The breakpoint was manually chosen so that we could estimate and compare the climate trend for the 40 years prior to 1990 with the trend for the recent 30 years. The regression trend results for each segment are also illustrated in green and cyan in these figures.

Water Flow Trends
Three major stations were used to represent each subregion, as they were the m downstream stations in each subregion, on the Athabasca River itself, with complete tasets. Stations 07AD002 (at Jasper, u3), 07BE001 (at Athabasca, m6), and 07DA001 (bel Fort McMurray, l2) were used as a proxy to represent the upper, middle, and lower A

Water Flow Trends
Three major stations were used to represent each subregion, as they were the most downstream stations in each subregion, on the Athabasca River itself, with complete datasets. Stations 07AD002 (at Jasper, u3), 07BE001 (at Athabasca, m6), and 07DA001 (below Fort McMurray, l2) were used as a proxy to represent the upper, middle, and lower ARB, respectively. Note that l2 was the only station in the lower ARB that had continuous and long-term measurement records, although it may not be a good proxy for the lower ARB as it is not located near the outlet of the lower ARB. The flow measured at l2 is the joint flow contribution from the Athabasca River and Clearwater River in the eastern part of the lower ARB. Station 07DD001 (at Embarras Airport, l3) lies on the Athabasca River at the river delta, near the outlet of the lower ARB. However, the water flow data at this station were only available for six years (1972,1973,1975,2015,2018, and 2019). Therefore, station 07DD001 was not included in the MK trend analysis but was used for comparing the water flow for these specific years. Three stations on the Peace River were also included in the MK test to compare with the Athabasca River. Stations 07GE001 (p1), 07GE001 (p2), and 07KC001 (p3) represented the upper, middle, and lower Peace River subregions, respectively.
The results of the water flow MK trend analysis and Sen's slope calculation are shown in Table 3 for the three main stations in the ARB. Table 3 shows the trend slope and significance level for each parameter and each station at each test period. The slope and significance for the remaining stations are shown in Tables S1-S8 in the Supplementary Materials. The results show that the cold season average water flow in the upper ARB has increased with a significance level of 95%. This trend was detected in all time periods between the past 30 and 65 years. The open warm season exhibited a variable trend across the warm months. The water flow in May has significantly increased since 1981 in the middle ARB, while the August water flow showed a statistically significant decline. It is worth noting that the water flow decline in August was no longer statistically significant after 1981. The overall annual average flow has decreased since 1956. However, it started to increase starting from 1981, although with a significance level of less than 90%. The map of the RAMP water flow gauging stations in the lower ARB and the water flow time series are shown in Figures S3-S7. The MK trends of the RAMP stations are shown in Tables S9-S15. Some RAMP stations investigated in this study duplicated WSC stations between 1997 and 2017, where the stations were physically the same. These stations were 07DA006-S38, 07DA008-S07, 07DA018-S39, 07DB001-S26, 07DC001-S27, and 07DD001-S46. Tables 4 and 5 show the total number of stations that demonstrated a significant (≥90% significance) increasing or decreasing trend in the water flow data. Table 4 was developed for the ARB stations, except for the station Athabasca River at Embarras Airport, due to the lack of sufficient water flow data for robust trend analysis. Table 5 was developed for the three stations along the Peace River. The arrows adjacent to each number indicate whether it is increasing or decreasing. The numbers in Tables 4 and 5 align with the observed trends in Table 3. As shown in Table 4, the number of stations with increasing water flow in the warm season was higher in more recent study periods. When including more years in the past, mainly before 1986, more declining trends appeared. This was also observed in the Peace River stations, as shown in Table 5. Jan  Jan The increasing trend in the open warm season flow has been observed in the middle and lower ARB. Figure 6 shows the chronological change in the open warm season average every year between 1991 and 2020. The year 2015 appears to be an inflection point, after which the water flow at Athabasca (07BE001) and Fort McMurray (07DA001) started to increase every year. That behavior was not present in the upper ARB, where the water flow remained relatively constant. The increasing trend was reflected in the MK trend analysis shown in Table 3. However, its significance was lower than 90%. Analyzing 30-year windows of data that move every five years for the ARB since 1956 showed that the definition of cold and warm seasons can vary according to the climate and hydrology of the region. Possible regime shifts can be detected when analyzing the 30-year moving windows of data, which can reveal the periods when significant changes occurred. Figure 7 shows the average hydrographs of each of the 30-year windows, as well as the cold and warm season averages and the annual average for the three major stations along the Athabasca River. The remaining stations are shown in Figures  S1-S2 in the supplementary information. The hydrographs show that the month of April did not exhibit significant increases in water flow over the cold season months in all upper ARB stations. In the middle and lower subregions, the April water flow was significantly higher than the March water flow due to factors, such as warmer climate and flow contribution from other tributaries. For this reason, the calculation of the cold season and warm season averages for the ARB in this study was adjusted to reflect this behavior. The cold season was November to April in the ARB, while the rest of the year was considered open warm season. This was true in the upper ARB only, where the water flow in April was consistent with freezing conditions. The middle and lower ARB had high water flow measurements in April, which led to moving April to the warm season in these two subregions. Analyzing 30-year windows of data that move every five years for the ARB since 1956 showed that the definition of cold and warm seasons can vary according to the climate and hydrology of the region. Possible regime shifts can be detected when analyzing the 30-year moving windows of data, which can reveal the periods when significant changes occurred. Figure 7 shows the average hydrographs of each of the 30-year windows, as well as the cold and warm season averages and the annual average for the three major stations along the Athabasca River. The remaining stations are shown in Figures S1 and S2 in the Supplementary Information. The hydrographs show that the month of April did not exhibit significant increases in water flow over the cold season months in all upper ARB stations. In the middle and lower subregions, the April water flow was significantly higher than the March water flow due to factors, such as warmer climate and flow contribution from other tributaries. For this reason, the calculation of the cold season and warm season averages for the ARB in this study was adjusted to reflect this behavior. The cold season was November to April in the ARB, while the rest of the year was considered open warm season. This was true in the upper ARB only, where the water flow in April was consistent with freezing conditions. The middle and lower ARB had high water flow measurements in April, which led to moving April to the warm season in these two subregions.
Stations 07GE001, 07GJ001, and 07KC001 are located on the Peace River with longterm and almost continuous historical records and in comparable locations to stations 07AD002, 07BE001, and 07DA001 along the Athabasca River, i.e., they are considered as a proxy for subregions of the Peace River Basin. Figure S2 includes the hydrographs of stations 07GE001, 07GJ001, and 07KC001. The hydrographs show that the Peace River has similar flowrates to the Athabasca River in the upper and middle subregions, but its lower subregion flowrate is higher. The Peace River had a consistent November-to-March cold season in its upper, middle, and lower reaches. Another observation was that the water flow in the cold months was significantly higher in the Peace River than in the Athabasca River.
The contribution of each of the upper, middle, and lower ARB subregions to the outlet water flow was calculated as percentages and shown in Figure 8. The only downstream WSC gauge station near the outlet of the Athabasca River, and with recorded water discharge data, was the Athabasca River Near Embarras Airport station (07DD001). The other stations used for representing the subregion contributions were the Athabasca River stations at Hinton (07AD002), Athabasca (07BE001), and below Fort McMurray (07DA001). Due to the presence of several data recording gaps, the only years where all four stations had complete annual recorded data were 1972,1973,1975,2015,2018, and 2019. These years were used to calculate the average water flow for each month and seasonal average and were plotted with the outlet water flow at Embarras Airport being 100%. Tables 6 and 7 show the average flow percentages split according to the decade of available data.  Stations 07GE001, 07GJ001, and 07KC001 are located on the Peace River with longterm and almost continuous historical records and in comparable locations to stations 07AD002, 07BE001, and 07DA001 along the Athabasca River, i.e., they are considered as a proxy for subregions of the Peace River Basin. Figure S2 includes the hydrographs of stations 07GE001, 07GJ001, and 07KC001. The hydrographs show that the Peace River has similar flowrates to the Athabasca River in the upper and middle subregions, but its lower subregion flowrate is higher. The Peace River had a consistent November-to-March cold season in its upper, middle, and lower reaches. Another observation was that the water flow in the cold months was significantly higher in the Peace River than in the Athabasca River.
The contribution of each of the upper, middle, and lower ARB subregions to the outlet water flow was calculated as percentages and shown in Figure 8. The only downstream WSC gauge station near the outlet of the Athabasca River, and with recorded water discharge data, was the Athabasca River Near Embarras Airport station (07DD001). The other stations used for representing the subregion contributions were the Athabasca River stations at Hinton (07AD002), Athabasca (07BE001), and below Fort McMurray (07DA001). Due to the presence of several data recording gaps, the only years where all four stations had complete annual recorded data were 1972,1973,1975,2015,2018, and 2019. These years were used to calculate the average water flow for each month and seasonal average and were plotted with the outlet water flow at Embarras Airport being 100%. Tables 6 and 7 show the average flow percentages split according to the decade of available data.  The contribution of each subregion to the downstream flowrate could be affected by physical characteristics such as topography, drainage patterns, and soil properties, as well as the land cover in the corresponding watershed. The LULC map shown in Figure 9 illustrates the prevalent land cover across the ARB in 2015. It is important to note that the LULC map is only relevant to the recent hydrology of the ARB, as historical LULC data were not available. The upper ARB is mainly mountainous, as it lies within the Rocky Mountains, with some forested cover. The middle ARB is covered with a mixture of grassland and forested areas. The lower ARB is also covered with grassland and forests, but the distribution is more uniform where the western side is more forested than the eastern  The contribution of each subregion to the downstream flowrate could be affected by physical characteristics such as topography, drainage patterns, and soil properties, as well as the land cover in the corresponding watershed. The LULC map shown in Figure 9 illustrates the prevalent land cover across the ARB in 2015. It is important to note that the LULC map is only relevant to the recent hydrology of the ARB, as historical LULC data were not available. The upper ARB is mainly mountainous, as it lies within the Rocky Mountains, with some forested cover. The middle ARB is covered with a mixture of grassland and forested areas. The lower ARB is also covered with grassland and forests, but the distribution is more uniform where the western side is more forested than the eastern side. The lower ARB is also home to the larger residential and industrial communities at Fort McMurray.

Discussion
In this study, a trend analysis was conducted on the water flow and climate data in the ARB region between 1956 and 2020. To the best of the knowledge of the authors, this study is the only one that included the recent records of the ARB (2010-2020). The study period of 1956-2020 for the water flow data provided a trade-off between a comprehensive historical analysis and the completeness of the datasets. Previous studies covered inconsistent study periods, which varied from 50 years [16,25] to more than 90 years [3,10]. In this study, the selection criteria of gauge stations were based on the presence of data gaps. It was found that records prior to 1956 were segmented and not suitable for robust trend analysis. In addition to the WSC data, the lower ARB water flow records from the

Discussion
In this study, a trend analysis was conducted on the water flow and climate data in the ARB region between 1956 and 2020. To the best of the knowledge of the authors, this study is the only one that included the recent records of the ARB (2010-2020). The study period of 1956-2020 for the water flow data provided a trade-off between a comprehensive historical analysis and the completeness of the datasets. Previous studies covered inconsistent study periods, which varied from 50 years [16,25] to more than 90 years [3,10]. In this study, the selection criteria of gauge stations were based on the presence of data gaps. It was found that records prior to 1956 were segmented and not suitable for robust trend analysis. In addition to the WSC data, the lower ARB water flow records from the RAMP dataset were also analyzed (see Supplementary Materials). While the RAMP data were also mostly seasonal and segmented, their results agreed with the results of the WSC data employed in this study.
The individual station water flow trend analysis demonstrated that the discharge in the Athabasca River has been steadily declining for the past 60 years, which is consistent with the literature [3,26]. The climate trends were also consistent with the hydrological behavior of the Athabasca River, where the precipitation has been declining across the ARB, with the most significant rate of −11.28 mm/decade at the upper ARB. The temperatures have also been steadily rising with a rate of 0.2-0.45 • C/decade across the ARB, leading to changes in the regional hydrology in terms of snowmelt, evapotranspiration, and groundwater table variation. The overall temperature within the ARB increased by about 2 • C during the past sixty years, particularly near Jasper. Previous studies on the climate trends between 1976 and 2010 were unable to detect significant trends in temperature [10,16], which could be due to several factors, such as pacific decadal oscillation, change in wet-dry cycles, and change in LULC [10]. The distribution of cooling and warming regions across ARB during warm seasons between 1976 and 2010 was approximately the same (see Figure 5 in [16]). In this study, when the period of record was expanded between 1950 and 2019 for climate data, the MK test revealed >99% significant increasing trends in temperature across the ARB, and most of the ARB showed significant warming trends during this period (see Figure 2). As Figures 3-5 show, the temperature had an increasing trend with a faster rate during 1950-1989 than during 1990-2019, particularly with a negative rate during 1990-2019 for the upper ARB, though it was not statistically significant.
The trend results in this study agreed with earlier studies, where a significant decrease in the precipitation trend and a significant increase in the temperature trend were estimated over ARB [16,26]. However, the present study provided more in-depth trend results both spatially and temporally compared to earlier studies [16,26]. For example, herein, it was found that the overall precipitation in the upper ARB was declining in the more significantly than the middle and lower ARB during the past 70 years. In [26], climate trends for ARB were estimated for two regions, northern and southern ARB, and declining precipitation and increasing temperature since 1960 were observed for both regions. In the present research, however, ARB was divided into three subregions, where each subregion is relatively homogeneous in terms of topography, climate, land cover, soil type, and biodiversity [20]. Table 2 and Figures 2-5 showed that the most significant declining precipitation was in the upper ARB, while the temperature was increasing for all the three subregions at a 99% confidence level since 1950.
The annual water flowrates for the three stations in Table 3 for the period 1961-2020 were approximately the same as the estimated flowrates shown in Table 3 in [26]. Other studies also estimated declining annual flowrates for longer historical periods [27]. However, in the present study, the flowrates for each month as well as cold and open warm seasons were also estimated, and it was observed that the flow in the upper ARB (07AD002, u3) was increasing during the cold months and cold season for all periods, 30-year to 60-year, particularly in March. In [27], it was reported that January and February for Hinton (u3) had an increasing flowrate, while July, August, and September showed a declining flowrate during 1955-2011, which is similar to the long-term trend analysis in this study (see Table 3). This shows that the gradual warming in the upper ARB had a profound impact on the melting of glaciers during cold months and particularly in early spring, contributing to water flow. While the flowrate had a declining trend for 07BE001 (m6) and 07DA001 (l2) during longer periods (e.g., 60-year), it was increasing significantly during the cold season in the recent period, 1991-2020. One reason is the water flow increase in the upper ARB during the cold season (minimum temperate rate of~0.4 • C/decade), which significantly contributed to the total basin flow [20]. Tables 6 and 7 show that the contribution of water flow from Hinton to the middle and lower Athabasca River was significantly higher in the recent years (2015, 2018, and 2019) compared to earlier years (1972, 1973, and 1975). Some studies [28,29] suggest that deep purple algae, likely because of forest fires, caused the glaciers to absorb more solar energy, which, together with winter warming, caused the rapid melting of glaciers over the last few decades. A significant positive flowrate in the tributaries of the Athabasca River was another major factor of flow increase in m6 and l2 during the recent decades. For example, Tables S1-S4 show a significant flowrate increase at a 99% confidence level since 1981 for 07CD001 (l1), a major tributary of l2. Figure 2 demonstrates a gradual warming trend and a declining precipitation rate near 07GE001 (p1) and 07GJ001 (p2), which are, respectively, stations at Wapiti River and Smoky River, one of the main tributaries of the Peace River. Tables S1-S3 show increasing trends in the cold season and the annual water flow at p1 and p2, whereas, for longer periods, the annual water flow had a declining trend for p1 and p2 though not statistically significant, which was similar to the results of the Athabasca River (see Tables S5-S7). The gradual melting of glaciers near p1 is likely one of the main contributing factors in water flow increase over the past few decades. Figure 7, Figures S1 and S2 demonstrate that the average water flow in the open warm season during 1981-2010 was relatively lower than other 30-year periods, particularly for stations 07GJ001 (p2), 07BE001 (m6), and 07DA001 (l2), where the water flow started to exhibit a positive gradient after 1981 (see Tables S1-S8). These results suggest the non-stationary behavior of water flow time series in the ARB and PRB, which is in line with the findings of Wu et al. [30].
Satellite remote sensing imagery acquired for 2001-2019 showed that grasslands and non-vegetated lands increased, while evergreen needleleaf forests decreased in the past two decades in the upper ARB [31], which is potentially a contributing factor to the high-altitude regional warming and water flow increase. Furthermore, an increase in water flow during 1991-2019 (see Table 3) is likely due to the recent changes in land cover of the ARB subregions. For example, from the same satellite imagery, it was observed that evergreen needleleaf forests in the middle ARB have increased by 4% per decade since 2001 [31]. Grasslands have also increased by 3% per decade since 2001 in the lower ARB, with a declining trend in forest cover [31]. Figures 3-5 show that the annual precipitation had a decreasing trend for upper and lower ARB, while it had an increasing trend for the middle ARB during 1991-2019, though not statistically significant at 90%. These gradual changes in precipitation might have had a significant impact on the water flow and forest cover of middle ARB, which was increasing, and forest cover of upper and lower ARB, which was decreasing during 2001-2019, as observed from satellite imagery. Deforestation by human and wildfires were some other factors of declining trends in forest covers [15,29,31].
When including the most recent water flow data, it was observed that the water flow had been steadily increasing between 2015 and 2020 in the middle and lower ARB. For some of the stations, the flow had an increasing trend using the MK test, which was less than 90% significant. The MK trends since 1952 demonstrated that, while the water flow was decreasing, the number of stations with significant declining slopes was also decreasing. At the same time, the stations with increasing significant trends were increasing in recent times. Several reasons could contribute to this behavior, including changes in land cover, groundwater, and human activities within the ARB. The 2015 LULC map displayed in Figure 9 shows an abundance of grassland and peatland in the middle ARB, which can influence the water balance in the watershed. A recent study presented a surface-subsurface flow and evapotranspiration model for the ARB to study the effect of the distribution of forests and peatland on the water balance in the ARB [12]. They reported that the peatlands in the middle and lower ARB reduce evapotranspiration and increase water availability, which aligns with the findings of this study.
The contributions of the upper, middle, and lower subregions to the total discharge of the Athabasca River, shown in Figure 8, reflect the factors affecting the river water flow. The upper subregion water flow behavior, captured by stations 07AD002 and 07BE001, demonstrates significant seasonal variability, which indicates that the water flow is mainly driven by climatic factors, such as temperature and precipitation, leading to snowmelt. The middle and lower ARB demonstrated less seasonal variation due to the year-long contributions of different tributaries as well as the interaction with groundwater, which diluted the effect of the upper subregion seasonal variation. The abundance of peatland and grassland across the middle and lower ARB may also contribute to the lower seasonal variability due to lowering the impact of evapotranspiration [12]. The long-term warming and long-term precipitation decrease, as observed in Figures 3-5, may potentially result in the reduction of water flow in the future, noting that glaciers in the upper ARB will disappear in the future due to their persistently negative mass balances [28].
The use of the MK trend test should be applied with care to any datasets characterized by the occurrence of seasonal regimes. Climate and water flow data are typically seasonal, with cyclic variations between warm and cold seasons. In this study, seasonality was eliminated by considering each month as a separate parameter, with the time steps being the years. In this case, each parameter represented the average of a specific month every year, which removes cyclic seasonal trends but retains other trends. Additionally, the correction of the MK test for serial correlation and regime shift has been a debatable topic in the literature [7,16]. A false trend can be detected in small datasets with small trend slopes due to linear serial correlations in the data. The Trend-Free Pre-Whitening (TFPW) correction can be applied to the dataset to remove the effect of serial correlation by distinguishing between the true trend and the serial correlation [32]. On the other hand, the TFPW correction was not recommended due to the significant power loss it causes, and it is not needed for sample sizes that exceed 50, as serial correlation was found to have a negligible effect on the rejection rate of the MK test in this case [33]. Natural variability in hydrologic data can cause regime shifts or a lack of continuity, which compromises the accuracy of the trend analysis that may be performed on the data [34]. Therefore, the data must be tested for the presence of regime shifts [35]. However, in the ARB, the water flow data did not exhibit signs of regime shifts and, therefore, was considered to reflect actual trends due to continuous climatic changes in the region [16].

Conclusions
In this study, trend analysis was conducted for the water flow and climate data in the ARB and PRB using the MK test and least-squares method. The data used in this study were collected from multiple locations across the entire ARB and PRB, with a period of record between 1956 and 2020 for water flow data and between 1950 and 2019 for climate data. The monthly, seasonal, and annual averaged trend analysis results for various historical periods were demonstrated. Possible interactions between climate, LULC, and water flow dynamics of the basins were also discussed. The findings of the study showed that, while the water flow trends in the past 60 years had a declining trend in the middle and lower subregions, the flowrate was increasing during the cold and open warm seasons over the past 30-40 years, which is likely because of the flowrate increase in the upper ARB and the recent climate and land cover changes in these subregions. The analysis of a longer climate record also confirmed an increasing trend in temperature and a decline in precipitation across the basins. Further research is recommended to analyze the contribution of the land cover and human activities to the water availability as well as the interaction between surface and groundwater in the different subregions of the ARB and PRB. It is also recommended to conduct the same analysis in the next five to ten years to confirm if the detected water flow increase will continue to occur. Finally, a comprehensive analysis of water availability across the entire ARB and PRB is recommended for future research, where a better understanding can be obtained regarding the contribution of climate, snowmelt, groundwater, and LULC to the Athabasca River and Peace River flows.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/hydrology9110197/s1, Figure S1: Hydrographs (left panels) and seasonal averages with standard deviation (right panels) of ARB stations; Figure S2: Hydrographs (left panels) and seasonal averages with standard deviation (right panels) of ARB stations (continued); Figure S3: Map of the RAMP and WSC gauging stations in the lower ARB. Some stations were replicated in the RAMP and WSC data under different IDs: 07DA006-S38, 07DA008-S07, 07DA018-S39, 07DB001-S26, 07DC001-S27, and 07DD001-S46; Figure S4: Hydrographs of the WSC stations between 1965 and 2020; Figure S5: Hydrographs of the RAMP stations between 1996 and 2017; Figure S6: Hydrographs of the RAMP stations between 1996 and 2017; Figure S7: Average flow percentages of hydrometric stations at Fort McMurray, Athabasca River below Eymundson Creek, and Embarras Airport for the years 2012, 2013, and 2014; Table S1: MK water flow slope and significance for a period of 30 years (1991-2020); Table S2: MK water flow slope and significance for a period of 35 years (1986-2020); Table S3: MK water flow slope and significance for a period of 40 years (1981-2020); Table S4: MK water flow slope and significance for a period of 45 years (1976-2020); Table S5: MK water flow slope and significance for a period of 50 years (1971-2020); Table S6: MK water flow slope and significance for a period of 55 years (1966-2020); Table S7: MK water flow slope and significance for a period of 60 years (1961-2020); Table S8: MK water flow slope and significance for a period of 65 years (1956-2020); Table S9: MK trend analysis of the RAMP water flow data for the time period of 1996-2017. Cells highlighted in red represent insufficient data for a robust analysis. Slope and significance (ˆ= 90%), (* = 95%), and (** = 99%); Table S10: MK trend analysis of the lower ARB water flow for the period of 1991-2020. Cells highlighted in red represent insufficient data for a robust analysis. Slope and significance (ˆ= 90%), (* = 95%), and (** = 99%); Table S11: MK trend analysis of the lower ARB water flow for the period of 1986-2020. Cells highlighted in red represent insufficient data for a robust analysis. Slope and significance (ˆ= 90%), (* = 95%), and (** = 99%); Table S12: MK trend analysis of the lower ARB water flow for the period of 1981-2020. Cells highlighted in red represent insufficient data for a robust analysis. Slope and significance (ˆ= 90%), (* = 95%), and (** = 99%); Table S13: MK trend analysis of the lower ARB water flow for the time period of 1976-2020. Cells highlighted in red represent insufficient data for a robust analysis. Slope and significance (ˆ= 90%), (* = 95%), and (** = 99%); Table S14: MK trend analysis of the lower ARB water flow for the period of 1971-2020.