Attribution of the Climate and Land Use Change Impact on the Hydrological Processes of Athabasca River Basin, Canada

Abstract

Climate change (CC) and land use/land cover change (LULCC) are significant drivers of hydrological change, and an effective watershed management requires a detailed understanding of their individual and the combined impact. This study focused on the Athabasca River Basin (ARB), Canada, and investigated how the basin responded to their changes using the MIKE SHE-MIKE Hydro River. Our findings revealed novel insights into ARB hydrological changes, including increment in non-vegetated lands (0.26%), savannas (1.28%), forests (0.53%), and urban areas (0.02%) while grasslands (2.07%) and shrublands (0.03%) decreased. Moreover, the basin experienced rising annual minimum (1.01 °C) and maximum (0.85 °C) temperatures but declining precipitation (6.2%). The findings suggested a significant impact of CC compared to LULCC as CC caused annual reduction in streamflow (7.9%), evapotranspiration (4.8%), and recharge (6.9%). Meanwhile, LULCC reduced streamflow (0.2%) and recharge (0.4%) but increased evapotranspiration (0.1%). The study revealed spatiotemporal variability across the ARB, with temperature impacts stronger in winter and precipitation influencing other seasons.

1. Introduction

The hydrological cycle is a naturally dynamic system comprising various hydrological processes driven by climatic factors, the watershed’s physical properties, anthropogenic activities, and their mutual interactions varying spatially and temporally [1]. Consequently, changes in climatic conditions and land use/land cover have been identified as key drivers of shift in the hydrological regime [2]. Climate change (CC) alters the pattern of precipitation and temperature in terms of intensity and frequency, significantly impacting the hydrological cycle. Documented effects include flow regime alteration, changes in timing and frequency of hydrological extremes, seasonal shifts, increased snowmelt, earlier spring freshet, modified evapotranspiration, soil moisture content, groundwater resources, and shifts in groundwater–surface water interaction [3,4,5,6,7]. Meanwhile, land use dictates the response of a watershed to climatic input, and its changes are associated with alteration in hydrological processes like interception, infiltration, evapotranspiration, surface runoff, and groundwater recharge rates [8,9] which ultimately affects the spatiotemporal availability of water in the basin.

Climate and land use change affect different aspects of the hydrological cycle, and thus their combined effects are complex, often interacting to either amplify or counterbalance each other’s impacts [10,11,12]. To develop effective, basin-specific water management strategies, it is essential to assess each stressor’s individual influence on the hydrological cycle, as well as their integrated effects, within each basin [13]. While some studies have found that LULCC can exert a stronger influence on hydrological processes, most findings attribute CC as the primary driver of change [2,14,15]. However, the extent of each stressor’s impact varies significantly by catchment [16]. For instance, urbanization was the dominant influence on runoff in six catchments in southeast Queensland, Australia [17] whereas in China’s Weihe River Basin, CC had a considerably stronger impact on evapotranspiration, soil moisture, and streamflow compared to LULCC [2]. These variations underscore the need for individualized basin assessments to accurately capture localized hydrological responses, enabling more tailored water management approaches.

The impacts of these driving factors can also exhibit spatiotemporal variability. For instance, Zeng et al. [18] revealed that CC accounted for greater changes in runoff during the dry season, whereas changes induced by human activities were more dominant in the wet season in the Zhang River basin in north China. Furthermore, most studies exploring the drivers behind hydrological changes, however, have primarily focused on streamflow alone, which may not provide a comprehensive depiction of the entire hydrological cycle’s dynamics [2]. Exploring other essential hydrological processes, such as evapotranspiration, overland flow, and groundwater recharge, can provide a more comprehensive understanding of the mechanisms governing the hydrological cycle and its responses to shifts in the surrounding environment. As water management increasingly emphasizes maintaining water availability across both spatial and temporal scales, it becomes essential to evaluate the impacts of these drivers on multiple hydrological components at diverse scales.

The ARB, home to Alberta’s multi-billion-dollar oil sands industry, is vital to the province’s economy, drawing considerable attention in hydroclimatic research due to its environmental and socio-economic significance. A consensus in past studies indicates declining river flow trend in the downstream [19,20,21] while noting increased flow at the headwaters. These trends have been linked with decline in precipitation [21,22,23], rising temperatures [22], and large-scale climatic oscillation cycles [24,25,26]. Predominantly, these findings are drawn from statistical analyses such as trend analysis [20], wavelet analysis [27], and multi-variate linear regression [22]. While hydrological models like SWAT [28], VIC [29], HydroGeosphere [30], and MISBA [31] have been applied in the ARB, they are primarily utilized to study future state of the basin’s hydrology. This study presents a novel application of hydrological modeling to assess the historical impacts of CC and LULCC on the basin’s hydrology, filling a critical gap in understanding past changes and their underlying mechanism in this vital region.

Hydrological modeling is widely used to distinguish individual impacts of stressors on hydrological systems by simulating baseline conditions and comparing them with scenarios and thereby isolating CC and LULCC effects [13]. Further, the modeling approach allows for evaluating effects at a higher temporal resolution as impact at daily, seasonal, and annual time scales can be assessed [18]. In this study, we have applied the coupled MIKE SHE-MIKE Hydro River model, which can comprehensively represent the primary hydrological process of the catchment. MIKE SHE is a deterministic, physics-based, fully distributed model that can demonstrate major hydrological cycle processes like evapotranspiration, overland flow, unsaturated flow, groundwater flow, and their interaction at varying spatial scales and complexity [32]. Meanwhile, MIKE Hydro River is a one-dimensional channel flow model for the simulation of rivers using fully dynamic Saint Venant equations. Numerous studies have used it to simulate the water balance of a watershed and predict its behavior in projected climate change, anthropogenic activities, and changes in land use [10,33,34]. In addition to this, the selection of a distributed hydrological model like MIKE SHE enables the evaluation to be concentrated and compared between smaller sub-watersheds within the basin.

Previous studies in the ARB have primarily focused on quantifying observed hydrological changes using statistical methods, with limited exploration of the underlying mechanisms [20,23,27]. Meanwhile, hydrological modeling efforts have largely been scenario-based, forecasting future changes rather than attributing past changes [3,28,29]. The objective of this study is to fill this gap by employing a physically based distributed hydrological model to ascertain the distinct and cumulative impacts of CC and LULCC on hydrological dynamics of the ARB. This study is the first to utilize the MIKE SHE in the ARB to specifically quantify the individual and combined effects of CC and LULCC on hydrological alterations, offering a novel perspective on hydrology of the ARB.

2. Materials and Methods

2.1. Study Area

The Athabasca River originates at the Columbia Icefield in Jasper National Park in Alberta, Canada, and flows northeast for more than 1500 km before finally draining into Lake Athabasca. The ARB encompasses an area of approximately 155,000 km² characterized by heterogenous hydroclimatic and physiographic properties. The elevation ranges from around 3600 m above mean sea level in the upper Rocky Mountains to around 200 m in the lower region. The river is the longest free-flowing river (with no regulating structures) in Alberta [35]. The flow regime is snow-dominated, with ice phenology playing a major role in the river flow [19]. Likewise, the region experiences four distinct seasons, winter (December to February), spring (March to May), summer (June to August), and fall (September to November). The major footprint of human activities includes agriculture, urbanization, livestock production, forestry, coal mining, oil and gas, and oil sands mining [36]. Although there has been a growing use of groundwater, the surface water from Athabasca and its tributaries serves as a primary source of water supply for anthropogenic activities [19,35,36]. The annual allocated water use is less than 5% of annual river flow, the bulk of which is withdrawal by the oil sands industry concentrated below Fort McMurray (Figure 1) [21,35,36].

2.2. Hydrological–Hydrodynamic Model

The availability of required data and computational efficiency were major factors that needed careful consideration for distributed models, particularly for a larger watershed like the ARB. These aspects influenced the choice of process and spatial resolution of process representation for the MIKE SHE-MIKE Hydro model. The relevant data and parameters used for the different modules, along with their details, are provided in the following sub-sections.

2.2.1. Model Domain and Grid

The topography of the model was defined by the digital elevation model obtained from the Canadian Digital Elevation Model (CDEM) at a resolution of 0.75 arc seconds. A cell size of 2000 m was selected as a computational grid representing the compromise between required information at the finest resolution and the computational requirement for the model.

2.2.2. Climate Input

Daily gridded observed precipitation, maximum temperatures, and minimum temperature from 1950 to 2019 available at a 0.1° × 0.1° grid size [37] were obtained from Alberta Environment and Protected Areas (EPA). The Hargreaves–Samani (HS) method was adopted to calculate the daily potential evapotranspiration [38]. The applicability of the HS method was validated by comparing with potential evapotranspiration data at four climatic stations in the area (Appendix A). The modified degree day method was used for the snowmelt, which required information about the degree day melting coefficient. The initial values of the monthly degree day coefficient were selected based on [39] and were adjusted during the calibration.

2.2.3. Surface Water Component

We used Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type 3 (LULC) at a 500 m resolution to represent the LULC in 2001 and 2015. The dataset collected and processed for quality control by [40] was used in this study. Altogether, nine different land classes available in the ARB were grasslands, water bodies, shrublands, savannas, deciduous broadleaf forests, urban areas, deciduous needleleaf forests, evergreen needleleaf forests, and non-vegetated lands. Leaf area index (LAI), root depth (RD), Mannings roughness coefficient (M), and detention storage (D) were required to define the properties of each land use class. These values were initially obtained from the literature [39,41,42], and some were calibrated accordingly based on the sensitivity of the parameters.

2.2.4. River Component and Coupling with MIKE SHE

The MIKE Hydro River model was set up to represent the hydrodynamic component of the model. The details on the river network and cross-section were obtained from Alberta Environment and Protected Areas. Altogether, 84 different branches with 514 cross-sections were extracted for the model and were used for the coupling of MIKE SHE and MIKE Hydro River. The interchange between the groundwater and river link was determined by the conductivity and head differential between the river and the grid cell.

2.2.5. Unsaturated Zone Component

The unsaturated zone was based on a two-layer water balance method. A soil map with a 1:5,000,000 spatial resolution was obtained from the Harmonized World Soil Database provided by the Food and Agriculture Organization (FAO), resulting in 22 different soil types for the study area. Saturated hydraulic conductivity, water content at field capacity, wilting point, and water content at saturation were input properties for each soil type adopted from the FAO database.

2.2.6. Saturated Zone Component

A simple linear reservoir model that represents the saturated zone by interflow and two parallel baseflow reservoirs was employed for the study. This technique, which lessens the amount of data required while boosting model execution time, has been proven especially helpful for large watersheds [43]. The time constants for the defined interflow and baseflow were calibrated during the model setup.

2.3. Model Calibration, Validation, and Sensitivity Analysis

The calibration and validation of the model were carried out against daily streamflow records for the period of thirty years from 1990 to 2019. The calibration of the model used LULC information from 2015 and a twenty-year period (2000–2019), while a ten-year period (1990–1999) was used to validate the model. The daily streamflow records from Environment Canada at three gauging stations (Athabasca, Fort McMurray, and Embarras Airport) along the Athabasca River (Figure 1) were used. However, flow data were only available for five years (2015 onwards) at Embarras Airport, so it was only used for model calibration. The performance of the model was evaluated using three statistical measures: Nash–Sutcliffe efficiency (NSE), percentage bias (PBIAS), and coefficient of determination (R²).

As the distributed physical model includes numerous parameters, reducing the number of free parameters through sensitivity analysis was critical. However, applying automatic calibration for sensitivity analysis in such a model is not desirable [44]. Thus, calibration parameters were selected that have already been consistently identified as sensitive parameters from various studies [39,43,45,46]. The selected parameters were degree day snow melting coefficient, Manning’s roughness coefficient for all land uses and riverbed, conductance between aquifer and river, and time constants for interflow and baseflow reservoirs. A local sensitivity analysis was conducted by individually varying each selected parameter within an appropriate range and evaluating the model’s performance under each scenario using three statistical metrics. Those parameters that caused significant variations in model outputs were identified as the most sensitive and critical.

2.4. Scenario Settings and Hydrological Components

The scenario analysis technique was used in this study to attribute the individual impact of CC and LULCC. The climate data from 1960 to 2019 was divided into two 30-year periods: 1960–1989 (called hereafter baseline period) and 1990–2019 (called hereafter recent period). The selection of a thirty-year period for comparison also allowed the study to be focused on evaluating the impact of CC rather than climate variability. This study used the earliest available (2001) LULC record from MODIS to represent baseline LULC conditions. The changes observed in LULC primarily stemmed from human activities, particularly those associated with population growth. An initial examination of population trends in the basin between 1981 and 2021, using census data from Canada and municipal data from Alberta, revealed significant growth, particularly in the early 2000s (Appendix B). This suggests that the LULC pattern of 2001 can serve as a reasonable representation of the vegetation characteristics from 1960 to 1989, given that more substantial changes occurred after 2001. On the other hand, the LULC map obtained from MODIS for the year 2015 (called hereafter recent LULC map), also used for calibration, was used to represent the recent period. Furthermore, an analysis of the annual trends in LULC coverage over the years revealed gradual linear changes [47], suggesting that using land use maps from 2001 and 2015 is sufficient to represent the intermediate changes during this period. The year-to-year variations are minimal, and modeling the hydrological impacts of such small LULC changes may not accurately capture their effects within the model. Overall, four scenarios based on two climate periods and the LULC map were developed for the study, as shown in

Table 1. Different scenarios with their climate period and LULC information used in the analysis.

Scenario	Climate Period	LULC	Remarks
S1	1960–1989	Baseline (2001)	Baseline
S2	1990–2019	Recent (2015)	Recent
S3	1960–1989	Recent (2015)	Hypothetical baseline
S4	1990–2019	Baseline (2001)	Hypothetical recent

. As such, two scenarios reflecting real world situations (baseline and recent) and two reflecting hypothetical scenarios were developed and simulated in the model.

The simulation results from individual scenarios were compared to identify the individual and combined impact of CC and LULCC for various selected hydrological indicators. The calculation was carried out by comparing the simulated value of hydrological indicators under different scenarios with respect to the baseline scenario’s value. For instance, percentage change due to individual stressors and their combined impact is calculated using the following equation.

$$\mathrm{Individual} \mathrm{impact} \mathrm{of} \mathrm{LULCC}={\displaystyle \frac{{\mathrm{S}}_{3\mathrm{i}}-{\mathrm{S}}_{1\mathrm{i}}}{{\mathrm{S}}_{1\mathrm{i}}}}\times 100$$

$$\mathrm{Individual} \mathrm{impact} \mathrm{of} \mathrm{CC}={\displaystyle \frac{{\mathrm{S}}_{4\mathrm{i}}-{\mathrm{S}}_{1\mathrm{i}}}{{\mathrm{S}}_{1\mathrm{i}}}}\times 100$$

$$\mathrm{Combined} \mathrm{impact} \mathrm{of} \mathrm{CC} \mathrm{and} \mathrm{LULCC}={\displaystyle \frac{{\mathrm{S}}_{2\mathrm{i}}-{\mathrm{S}}_{1\mathrm{i}}}{{\mathrm{S}}_{1\mathrm{i}}}}\times 100$$

where S_1i, S_2i, S_3i, and S_4i are the simulated values of the ith hydrological indicators under baseline, recent, hypothetical baseline, and hypothetical recent scenarios, respectively.

To depict the overall impact on the hydrological system, the analysis has been carried out for three different hydrological processes; streamflow, evapotranspiration, and groundwater recharge for annual and monthly time scales. Furthermore, as the streamflow is widely regarded as a master variable capable of reflecting different aspects of the hydrological cycle, we have analyzed not only the annual average magnitude of streamflow but also its annual extremes feature. Therefore, in addition to annual and monthly average values, we have selected four additional streamflow indicators from the indicators of hydrological alteration (IHAs), which include annual 7-day maximum, annual 7-day minimum, date of annual maximum, and date of annual minimum flow. The first two indicators show the maximum and minimum flow observed for seven consecutive days in a water year and reflect the extreme conditions in terms of high and low flow in the basin. The latter two indicators represent the timing of extreme flows as the date of occurrence of annual maximum flow and annual minimum flow for each water year. In order to address the spatial variability that may be present between areas with high human activities and with low human activities, this study has divided the basin into upper and lower regions, with Fort McMurray as the boundary (Figure 1). This allowed results to be evaluated for the upper region and lower region separately in addition to the entire watershed.

3. Results

3.1. Evolution of Land Use Pattern in the Basin

The distribution of LULC within the basin over two selected periods of 2001 and 2015 is presented in Figure 2. In general, it can be observed that savannas and forests of various types dominate over three-quarters of the total area in the basin. The areas of savannas were found have increased in 2015 compared to 2001. The distribution of urban areas demonstrates that human settlements are concentrated in a few urbanized areas, which have increased from 0.04% in 2001 to about 0.06% of the basin’s area in 2015. This corresponds to an almost 50% increment in urban areas. Likewise, areas of evergreen forests and non-vegetated lands have also increased, while deciduous forest, grassland, and shrubland coverage decreased in these periods.

The basin’s overall landscape exhibits three distinct properties that align with the natural region division of Alberta [48]. The uppermost areas are characterized by the Rocky Mountains that remain mostly untouched by human activities, resulting in little to no change in land use patterns. The upper central part lies in the foothills region, with a gentle slope consisting of coniferous and mixed wood forests. The boreal forests cover the central and lower parts of the basin and play a pivotal environmental role in maintaining the natural biodiversity of the region [49,50]. The study found that while the upper region exhibited changes in LULC similar to those across the entire watershed, the lower region showed distinct patterns. From 2001 to 2015, the coverage of grasslands, savannas, and non-vegetated lands increased by about 5.5%, 4.0%, and 0.3%, respectively, in the lower region (

Table 2. Percentage coverage of LULC class in 2001 and 2015 for the entire watershed, upper region, and lower region.

LULC Class	Entire		Upper		Lower
	2001	2015	2001	2015	2001	2015
Water Bodies	1.76	1.75	1.96	1.94	0.57	0.61
Grasslands	12.26	10.19	13.92	10.54	2.83	8.36
Shrublands	0.04	0.01	0.01	0.01	0.15	0.04
Savannas	44.88	46.16	42.13	42.92	60.41	64.27
Evergreen Broadleaf Forests	4.68	6.35	4.91	7.17	3.35	1.69
Deciduous Broadleaf Forests	35.60	33.85	36.21	35.68	32.51	23.87
Evergreen Needleleaf Forests	0.02	0.63	0.02	0.61	0.03	0.72
Non-Vegetated Lands	0.73	0.99	0.81	1.08	0.14	0.44
Urban Areas	0.04	0.06	0.04	0.07	0.01	0.01

). This is at the expense of broadleaf forests, whose coverage area was reduced by more than 10.0%. The oil and gas industries are mainly located in the lower region of the basin, leading to intense urbanization in these areas. Consequently, the lower region is very dynamic in terms of landscape changes. In addition, the basin is frequently subjected to wildfires, which have also contributed to alterations in the land use pattern of the area [51].

3.2. Climate Change Observed in the Basin

CC in the basin exhibits a trend of reduced precipitation while increasing maximum and minimum temperatures (Figure 3). The average annual precipitation in the basin during the baseline period was 536 mm, and it reduced to 503 mm (−6.2%) during the recent period. Meanwhile, the average maximum temperature in the basin increased by 0.85 °C from its baseline value of 6.03 °C. The change in the average minimum temperature was even more pronounced, with a rise of almost 1.00 °C observed (increased from −5.35 °C to −4.34 °C) during the two periods.

The spatial distributions of observed climate variables in the baseline and recent periods, along with changes between these two periods, are shown in Figure 4. It is evident that the upper part of the basin (southwest) receives higher precipitation than the lower part. Although there is an overall decrease in precipitation, some scattered areas in the lower basin have experienced increased precipitation. The central basin is considerably warmer than other parts, with higher maximum and minimum temperatures. Surprisingly, the lower part of the basin, despite being at a lower elevation, is colder and experiences temperatures like those of the mountainous areas upstream. The comparison of temperature during the two periods reveals that the lower part of the basin is warming at a faster rate than other regions. Changes in the maximum and the minimum temperature in the lower part of the basin are up to 1.40 °C and 2.00 °C, respectively, in the considered period. This temperature increase aligns with findings by [23] who observed a consistent uptrend in maximum (0.23 to 0.28 °C per decade) and minimum (0.17 to 0.45 °C per decade) temperatures from 1950 to 2019 in the ARB.

The seasonal distribution of the precipitation shows that the basin receives higher precipitation during the summer period, with almost 50% of annual precipitation falling between June and August (Figure 5). The average monthly precipitation in July (97 mm) was highest during the analysis period. Meanwhile, the precipitation during the late fall to early spring is lower and the least is generally observed in February (20 mm). The comparison of the two periods reveals that precipitation in the basin decreases throughout the year except during April (13.68% increment) and June (5.09% increment). The change is particularly significant during the winter, with a reduction of up to 21% (4.75 mm) observed in February. Although the precipitation is still higher in summer and early fall months, they also experience a reduction ranging from almost 3.75 mm in July to 9.15 mm in September.

However, sub-watershed analysis reveals that the precipitation changes in the lower part of the basin do not necessarily reflect the overall basin’s changes. This is particularly visible during the late spring, late winter, and spring months, where the changes in lower region precipitation were opposite to that of the entire watershed and upper regions. For instance, the precipitation during October reduced slightly by 1 mm on average throughout the basin, but it declined significantly by almost 20 mm in the lower region. At the same time, precipitation increased by about 2 mm in the upper region, further highlighting the spatial variability of CC in the basin.

The seasonal changes in the two periods highlighted that temperature is increasing throughout the year, with only October experiencing a reduction (−0.74 °C and −0.06 °C in maximum and minimum temperature). The warming tendency was more profound in winter, with a rise of more than 3 °C observed during January for both maximum and minimum temperatures. Although the annual change in minimum temperature was higher than the maximum, summer months experienced a greater change in maximum temperature.

3.3. Performance of the Hydrological Model

The result of model calibration and validation against daily streamflow at three gauging stations is shown in Figure 6. The model performed well in capturing both high and low flows of the streamflow variability. During the calibration and validation periods, the NSE value for all three stations was greater than 0.70, indicating a very good performance [52]. The coefficient of determination was also above the acceptable range, ranging from 0.71 to 0.77 during calibration and 0.72 to 0.77 during validation. There was a slight overestimation bias at the Embarras station (negative PBIAS), while the Athabasca and Fort McMurray stations showed underestimation bias ranging from 0.34% to 11.48%. Overall, the model performed satisfactorily, and the model performance was comparable to other models previously applied in the basin [3,28,29].

Among the selected parameters, the time constants for both interflow and baseflow reservoirs were found to be highly sensitive. These time constants directly govern the contribution towards streamflow in the form of interflow and baseflow in the linear reservoir method, resulting in higher sensitivity. Meanwhile, Manning’s coefficient, degree day melting coefficient, and conductance between river and aquifer were found to be less sensitive compared to the saturated zone’s parameter. The final calibrated values of parameters have been provided in Appendix C.

3.4. Impacts of CC and LULCC on Hydrological Components

The analysis of the annual change in hydrological indicators clearly shows that CC significantly influenced the hydrological cycle, surpassing the impact of LULCC (

Table 3. Average annual value of streamflow indicators, evapotranspiration, and groundwater recharge for all simulated scenarios with the contribution of CC, LULCC, and total change for the upper, lower, and entire watershed. Note. The values in parentheses represent percentage contributions.

Variables	Region	Mean Value in Each Scenario				Absolute Contribution of Each Factor (%)
		S1	S2	S3	S4	CC	LULCC	Combined
Annual streamflow (m3/s)	Entire	675.0	620.9	673.6	622.1	−52.9 (7.9%)	−1.5 (0.2%)	−54.1 (8.0%)
	Upper	574.9	527.6	573.5	528.7	−46.2 (8.0%)	−1.4 (0.2%)	−47.3 (8.2%)
	Lower	107.5	99.5	107.5	99.4	−8.1 (7.6%)	−0.1 (0.01%)	−8.0 (7.4%)
Annual 7-day minimum flow (m3/s)	Entire	169.3	162.4	168.6	161.5	−7.8 (4.6%)	−0.6 (0.4%)	−6.8 (4.0%)
	Upper	143.1	135.9	142.5	136.4	−6.7 (4.7%)	−0.5 (0.4%)	−7.2 (5.0%)
	Lower	11.7	12.6	12.0	13.1	1.4 (12.4%)	0.3 (2.5%)	0.9 (7.6%)
Annual 7-day maximum flow (m3/s)	Entire	2442.1	2042.7	2431.3	2088.4	−353.6 (14.5%)	−10.8 (0.5%)	−399.4 (16.4%)
	Upper	2160.3	1826.4	2151.1	1831.6	−328.7 (15.2%)	−9.2 (0.4%)	−333.9 (15.4%)
	Lower	511.3	426.1	510.0	425.4	−86.0 (16.8%)	−1.3 (0.3%)	−85.3 (16.6%)
Annual date of minimum flow (day)	Entire	63.0	58	64	57	−6.0 (9.5%)	1.0 (1.6%)	−5.0 (7.9%)
	Upper	63.0	58	64	57	−6.0 (9.5%)	1.0 (1.6%)	−5.0 (7.9%)
	Lower	172.0	173	173	171	−1.0 (0.6%)	1.0 (0.6%)	1.0(0.6%)
Annual date of maximum flow (day)	Entire	184.0	186.0	186.0	185	1.0 (0.5%)	2.0 (1.1%)	2.0 (1.1%)
	Upper	182.0	183.0	183.0	180	−2.0 (−1.1%)	1.0 (0.5%)	1.0 (0.5%)
	Lower	190.0	187.0	191.0	186	−4.0 (2.1%)	1.0 (0.5%)	−3.0 (1.6%)
Annual evapotranspiration (mm)	Entire	384	365.6	384.3	365.4	−18.6 (4.8%)	0.2 (0.0%)	−18.4 (4.8%)
	Upper	390.5	371.5	391.3	370.8	−19.7 (5.0%)	0.7 (0.2%)	−19.0 (4.9%)
	Lower	346.6	331.8	344.1	334.2	−12.4 (3.6%)	−2.5 (0.7%)	−14.8 (4.3%)
Annual groundwater recharge (mm)	Entire	144.3	133.8	143.7	134.3	−10.0 (6.9%)	−0.6 (0.4%)	−10.5 (7.3%)
	Upper	148.3	136.5	147.4	137.3	−11.0 (7.4%)	−0.9 (0.6%)	−11.8 (8.0%)
	Lower	121.4	118.2	122.5	117.1	−4.3 (3.6%)	1.1 (0.9%)	−3.2 (2.7%)

). For annual streamflow, CC caused an almost 40 times greater reduction than LULCC, with an 8% decrease versus a mere 0.2% decline. The upper and lower regions exhibited a more substantial impact of CC, with a higher impact observed in the former. On the other hand, seasonal variances in CC’s impact on streamflow were observed (Figure 7). Spring months (March and April) experienced increased flow due to CC, while reductions were noted in other months, declining up to 17% in May. The changes due to CC were similar throughout the basin area. Meanwhile, LULCC had a clear semiannual pattern as it decreased from January to June (0.11 to 4.62%) while it increased in the latter half of calendar year (0.37 to 1.80%). The spatial variation in LULCC impact was evident, with the lower region showing inconsistent seasonal streamflow changes compared to the upper region.

Similar patterns are observed in annual extremes of streamflow. CC has contributed to a significant reduction in both the annual 7-day maximum (14.5%) and minimum (4.6%) flows, while LULCC had minimal impact. However, minimum flow in the lower region showed an opposite pattern, increasing due to CC and LULCC by approximately 1%. Notably, the reduction in minimum flow in the lower region due to LULCC was proportionally higher compared to other LULCC-induced changes. In the case of extreme timings, low flow was found to be shifting earlier due to CC (6 days) while it was delayed by one day due to LULCC. As a result, the low flow date was found to be shifted by five days ahead due to both stressors. Likewise, for maximum flow timing date was delayed by CC (1 day) and LULCC (2 days). However, spatial variation was observed in the case of CC’s impact on high flow timings, as it was found to be delayed in upper (2 days) and lower (4 days) regions.

In contrast to streamflow, CC and LULCC had counteracting impacts on annual evapotranspiration. CC reduced evapotranspiration for the entire watershed by about 4.9%, while LULCC had a slight increasing effect of 0.1%. The sub-watershed analysis, however, shows that the changes observed are not consistent throughout the basin. Evapotranspiration has been found to be increasing in the upper region (0.2%) but declining in the lower region (0.7%) due to LULCC. The seasonal pattern of evapotranspiration follows a similar trend to precipitation, with higher rates during summer and lower rates during winter. A notable observation is the overall increase in evapotranspiration under CC from late fall to early spring, with the most significant percentage increment occurring in January, approximately 32%. On the other hand, LULCC caused a reduction in evapotranspiration in winter and spring months, while it increased in summer and fall. Sub-watershed analysis highlighted differences between the upper and lower regions’ responses to LULCC impact, with evapotranspiration increasing throughout the year in the upper region and decreasing throughout the year in the lower region.

Both CC (6.9%) and LULCC (0.4%) caused a reduction of groundwater recharge for the entire watershed. The reducing impact of CC was consistent throughout the basin. As seen for other hydrological processes, the lower region had varying responses to LULCC as recharge was found to have increased by 0.9%. Seasonal changes showed reduced groundwater recharge due to CC during summer and fall in the upper region, which was consistent with the entire watershed. At the same time, CC contributed to increased recharge during late fall to early winter (12% to 14% per month) in the entire watershed. However, as with the case of other hydrological components, the lower region’s response to CC showed a few differences with the upper region. For instance, the recharge in the lower region mostly increased in summer as opposed to declining summer recharge in the upper region due to CC. Meanwhile, LULCC resulted in decreased recharge (0.2% to 0.8% per month) throughout the year in the upper region, which was similar for the entire watershed. In contrast, the recharge rate in the lower region increased due to LULCC except during winter.

4. Discussion

4.1. Impact of CC and LULCC

CC has been found to be highly dominant in inducing the observed hydrological changes in the basin. This finding of a higher impact of CC aligns with various previous results from other regions around the globe [53,54,55,56]. For instance, a study by Ahmed et al. [54] found that CC caused an almost 98% increment in annual runoff while LULCC only contributed a 2% rise at the Yangtze River Basin, China. Similarly, Tan et al. [15] reported that CC had a higher impact on both annual streamflow and evaporation in the Johor River Basin, Malaysia. In our study, the primary influence of CC on streamflow has predominantly been caused by decreased precipitation, resulting in less water availability for other hydrological processes. This led to a reduction of flow in rivers and streams. Likewise, less water availability for evapotranspiration might have resulted in the decline of annual evapotranspiration. These findings are consistent with other studies in the ARB highlighting precipitation as the major driver of change [20,23,27]. The study by Bawden et al. [20] used a multiple linear regression (MLR) model to analyze the relationship between runoff, precipitation, and temperature in the ARB. They found that warm season precipitation trends explained over twice as much variation in runoff as temperature trends, with R² values of 62% and 30%, respectively. They also identified a significant correlation between precipitation and streamflow trends during the warm season (March–October) from 1976 to 2010, at a 10% significance level.

The seasonal impact of CC, however, shows that the impact is quite complex and dynamic in nature. The findings suggest that the governing climate variable varies with season. During the late winter and early spring, it has been seen that the temperature drives the hydrological processes as snowmelt forms a significant portion of the streamflow. In contrast, precipitation is the driving factor during the other periods of the year. This seasonal variability explains the increase in streamflow due to higher snowmelt from warming temperatures in March and April. These months are generally when spring freshet, a key event of cold region hydrology, is initiated in the watershed. The increasing annual 7-day minimum flow observed in the lower region can also be linked to this effect. As the lower region is warming at a higher rate, its impact has been directly reflected in increased low flows, which usually occur during winter and early spring. Likewise, the diminishing precipitation during the relatively warm season has direct implications for the streamflow generation, resulting in reduced flow. Furthermore, the warming tendency during late winter and early spring has likely contributed to increased rain-on-snow events. These events involve rainfall directly onto the snowpack, creating more overland flow, resulting in greater overland flow and consequently augmenting water availability in the river. Rising winter and spring temperatures also have cascading effects on summer at a lower elevation. The diminished snowpack due to snowpack sublimation in winter leads to reduced water storage, ultimately limiting the supply to the river in the summer.

The seasonal distribution and the governing mechanism of evapotranspiration are similar to those of streamflow. The reduced precipitation during late spring and fall has likely resulted in the reduction of ET during these months. However, during the winter and early spring, the increased temperature has accelerated the snowmelt, which in turn has allowed for increased evapotranspiration during this period. A similar finding of a possible increment in winter evapotranspiration due to rising temperatures in the future was also reported by Farjad et al. [10] in the Elbow watershed, Canada.

This study has shown that the impacts of LULCC in the hydrological processes have been manifested through reduced streamflow, alteration of peak flows, increased evapotranspiration, and reduction of groundwater recharge. During the specified period, there was an observed increase in the coverage of savannas and forests, resulting in a notable expansion of canopy areas conducive to evapotranspiration. This phenomenon likely contributed to the overall rise in annual evapotranspiration rates, while concurrently leading to a reduction in both annual streamflow and recharge. Despite minor increments in non-vegetated and urban lands, these changes were insufficient to counteract the prevailing trends. Notably, in the lower region, a significant decrease in deciduous forests was noted, partially explaining the increased evapotranspiration recharge in that area. Moreover, the increment in grasslands, non-vegetated lands, and developed areas introduced variable responses in terms of evapotranspiration and recharge dynamics. The landscape of the ARB is still undergoing changes at a decent rate [57]. The upper region of the basin experienced a reduction in agricultural lands and increased rural and urban settlements from 2011 to 2022. These changes are expected to continue in the future, which might have a significant impact on the basin’s hydrology. Although the annual impact might not be significant, it can still significantly alter the timing and frequency of extreme flows. Moreover, changes in vegetation properties can have implications for the seasonal distribution of streamflow as well as other important hydrological processes.

The integrated impact of both stressors has resulted in seasonal disturbance in the hydrological cycle. These induced seasonal shifts have significant implications for aquatic ecosystems and water availability for various industrial activities, such as utilities, mining, and bitumen extraction in the lower part of the basin. The flow alteration in the ARB has cascading impacts on the downstream Peace-Athabasca Delta (PAD), which provides habitat for 45 mammal species, 214 bird species, and 20 fish species [58]. The flow shifts in the Athabasca River, alongside other nearby rivers, influence the formation of perched basins and their periodic flooding, processes vital for the survival of species such as muskrats, waterfowl, and various other wildlife. The lower Athabasca River supports year-round spawning, rearing, and feeding habitats for various fish species [35], with any changes in flow potentially impacting their distribution. For instance, a study by Morales-Marin et al. [59] found that reduced streamflow results in lower water velocities and depths, creating sub-optimal conditions for the Athabasca rainbow trout (ART), a species classified as “at risk”.

Water consumption from the Athabasca River has surged significantly in recent decades [19], posing potential challenges for its users. This is particularly critical for the oil sands industry, which depends heavily on the river, accounting for over 70% of the total water usage [35]. To address this, a framework has been implemented to regulate cumulative water withdrawals, ensuring a balance between human and environmental demands [35]. The framework establishes weekly management triggers and withdrawal limits based on the seasonal flow variations. However, as climate and land use changes increasingly affect the flow regime, updates to the framework may be necessary to balance social, environmental, and economic priorities in the future. Different climate impact assessment studies carried out in the ARB have already reported that the flow regime is likely to experience shifts in the future in the form of increased annual flow, higher seasonal fluctuations, shifts in hydrography, and altered peak flow frequency, snow coverage, glacier contribution, and more which might result in spatiotemporal heterogeneity in water availability leading to water security challenges [3,28,29,60,61]. With anthropogenic activities expected to further alter the landscape properties, their integrated impact might lead to more uncertainty in the spatiotemporal water availability in the basin.

4.2. Uncertainties

Hydrological models have been widely recognized as significant contributors to uncertainties when simulating hydrologic processes. These uncertainties can arise from multiple sources, such as input data, parameter selection, the degree of process simplification, calibration, and structural uncertainty [62,63]. Furthermore, studies have observed that there are additional uncertainties associated with modeling at high altitudes [64] as well as the inadequate accounting of river ice cover dynamics [29] in most hydrological models. The hybrid climate datasets utilized in this study were developed using the REFerence Reliability Evaluation System (REFRES), which produces high-resolution hybrid climate data by mathematically ranking multiple datasets [37]. While these hybrid datasets offer extensive spatial coverage, they may not be as reliable as ground-based observational data. The authors also highlighted that the reliability of the climate data diminishes at higher altitudes, potentially due to factors such as topographic variability, orographic effects, and the limited availability of ground-based measurements.

Additionally, the selection of only two LULC maps, each representing one period of analysis, has essentially limited the temporal representation of the watershed’s vegetation property. This limitation suggests the possibility of underestimating or overestimating LULC features, leading to uncertainties in the simulated outcomes [2]. The model calibration was carried out by fine tuning only a few selected parameters, and ignoring other parameters may also induce uncertainties in the modeling results. Moreover, the calibration and validation of the model against daily streamflow are likely to introduce uncertainties in the model performance inherently. The calibration against daily streamflow may particularly induce uncertainties in a cold climate region like the ARB. During the winter season, significant ice coverage occurs in the river, which is not sufficiently addressed in MIKE SHE, potentially contributing to uncertainty. Additionally, the presence of ice coverage may have impacted the accuracy of the streamflow data used for calibration.

Future studies can address these uncertainties by adopting a more robust modeling framework. Employing physics-based methods to better capture snowmelt and thaw processes, incorporating river ice dynamics, and providing a more detailed representation of the underlying saturated zone could enhance model performance. Additionally, calibrating the model using multiple variables, such as soil moisture and groundwater levels, would improve its reliability and accuracy.

5. Conclusions

This study assessed the individual and integrated effects of CC and LULCC on the hydrological system of the ARB using the MIKE SHE-MIKE Hydro River model. The results showed that MIKE SHE can adequately represent hydrological conditions of cold climate regions like ARB and the model could be applied to other similar climatic regions with adequate calibration and validation. The findings clearly indicate that there have been considerable changes in climatic conditions and LULCC patterns in the ARB. These changes have altered the basin’s hydrological system in which the CC impact emerged as a dominant factor of change.

The evaluation of CC revealed a warming and drying trend in the region, with the lower part of the basin witnessing a greater rise in temperature. The annual maximum and minimum temperatures have increased by 0.85 °C and 1.01 °C, while precipitation has decreased by 33 mm in the recent thirty-year period compared to the earlier baseline from 1960 to 1989. The findings suggest that CC alone would have been responsible for reducing annual flow along with its extremes, annual actual evapotranspiration, annual recharge, and annual baseflow. The hydrological response of the basin to climatic factors was found to be dynamic as temperature fluctuations had greater implications during winter, while precipitation fluctuations governed climatic impacts in other seasons. This also led to seasonal shifts in the streamflow, overland flow evapotranspiration, recharge, and baseflow of the basin.

At the same time, the basin was subjected to increasing anthropogenic activities that led to an increment in non-vegetated lands, forested areas, and developed areas at the expense of grasslands and shrublands. The modeling result suggests that these changes induced an increment in annual evapotranspiration and recharge. Sub-watershed analysis in the upper and lower regions highlighted the heterogeneity in the basin’s response to these stressors, further complicated by seasonal variations in the dominant stressor. The lower region had different LULCC patterns, and this was also reflected in its hydrological response. These findings emphasize the complexity of the basin’s hydrological system. The study demonstrated that both CC and LULCC have affected the hydrological processes of the Athabasca River Basin in various ways. Neglecting either of these factors in assessing the basin’s current or future hydrological state could lead to an overestimation or underestimation of the hydrological response. Consequently, it is crucial to consider the potential changes in the hydrological system due to both stressors when formulating water management practices and policies for the foreseeable future.