Modelling groundwater quality of the Athabasca River Basin in the subarctic region using a modified SWAT model

Groundwater is a vital resource for human welfare. However, due to various factors, groundwater pollution is a paramount environmental concern. It is challenging to simulate groundwater quality dynamics with the Soil and Water Assessment Tool (SWAT) because it does not adequately model nutrient percolation processes in the soil. The objectives of this study were to extend the SWAT module to simulate groundwater quality for the parameters nitrate and Total Dissolved Solids (TDS). The results of the SWAT model for the Athabasca River Basin in Canada revealed a linear relationship between observed and calculated groundwater quality. This result achieved satisfactory values for coefficient of determination (R2), Nash-Sutcliffe efficiency (NSE), and percent bias (PBIAS). For nitrate, the model performance measures R2 ranged from 0.66–0.83 during calibration and NSE from 0.61–0.83. R2 is 0.71 during validation and NSE ranged from 0.69–0.75. Likewise, for TDS, the model performance measures R2 ranged from 0.61–0.82 during calibration and from 0.58–0.62 during validation. When coupled with soil zone and land surface processes, nitrate and TDS concentrations in groundwater can be simulated with the SWAT model. This indicated that SWAT may be helpful in evaluating adaptive management scenarios. Hence, the extended SWAT model could be a powerful tool for regional-scale modelling of nutrient loads, and to support and effective surface and groundwater management.


Methods and materials
The study area. The research area is in Athabasca River Basin (ARB), which is located in the central part of the province of Alberta in the subarctic region ( Fig. 1). Athabasca River is the second largest river in Alberta, and the largest undammed river. The mean annual discharges in cubic decameters (dam 3 = 1000 cubic meters) at the points along the river are 2,790,000 dam 3 at Jasper 13,600,000 dam3 at Athabasca. The river basin in general have substantial economic contribution in Canada by providing reliable freshwater supply to the people as well as for various industries, such as oil sands mining and pulp mills 42 . The main land cover type of the ARB is the boreal forest, which shares 82% of the total land; agriculture shares 9.5%, which is in the central portion of the watershed (i.e., Pembina, Lesse Slave, McLeod and upper parts). Overall forest, agriculture, traditional oil and gas extraction, oil sand mining, and coal mining are the major industries of the ARB 41 . Mean annual precipitation of the area ranges from 300 mm in the lower portion of the river basin to > 1000 mm from the headwaters, while the mean temperature ranges 1.8-5.1 °C 43 .
The major aquifer types include near-surface sands, buried channels and valleys, Paskapoo aquifers, and bedrock aquifers, which characterize the basin with primary and/or secondary porosities 42 . In this area, groundwater within unconsolidated surficial deposits consist of a mixture of a preglacial and glacial sand and gravel aquifers. In the Athabasca River Basin, groundwater recharge is spatially variable, and groundwater contributes baseflow to rivers and streams 44,45 . This region has extensive wetlands, including groundwater-fed wetlands (fens), and one of the world's most ecologically significant wetlands designated by the Ramsar Convention. The ecology within the watershed is diverse as a result of the different natural regions. As the major river system, the ARB is influenced by a variety of climate, terrain and landscape characteristics found within its basin. The seasonality of climatic conditions is a major influence on river flow conditions. Soil and Water Assessment Tool (SWAT) description and application. The  www.nature.com/scientificreports/ complex watersheds at daily time step 46,47 . It could be used for simulating the hydrological process, nutrient concentration, runoff generation, and sediment yield in the watershed under various landuse as well as soil scenario 46,48,49 . Additionally, it is used for simulating the effects of climate change, landuse and management practices on the quantity and quality of water 50 . The SWAT has been used for simulating the N cycle in soils and aquifer at shallow depth 51,52 . In the water and soil, the N processes is dynamic. It could be added to the soil in the form of manure, bacteriological fixation, residue, and rainfalls. It may be moved out from the soil through by plant uptake, volatilization, soil erosion, denitrification, and leaching. In the SWAT model, there are various types of N pools, including inorganic forms of nitrogen, and organic forms. NO 3 − may transport with runoff, percolation or lateral flow and recharge into the aquifer to the shallow depth from the soil profile. In addition, NO 3 − can move with groundwater flow to river channels or be carried out of the shallow aquifer into the soil zone during water deficiencies. The amount of nitrate carried by water is calculated by multiplying the concentration of nitrate in the mobile water fraction by the volume of water moving in each route. The amount of NO 3 − discharge depends on NO 3 − concentration in the soil-water domain. Nitrate uptake by plants has an inverse exponential relationship with depth 47 . where NO 3conc,z refers to initial nitrate concentration (mg/kg) at the depth of z (mm). Only the parts of NO 3 − is mobile and therefore it available for the discharge through the tiles. To calculate the concentration of nitrate in the mobile water fraction, the following equation has been adopted from 53,54 . where ConcN O3,mobile refers to NO 3 − concentration in the movable water at a given layer (kg N/ha), W mobile represent the amount of water lost by runoff, percolation or side flow in the layer (mm H 2 O), θe represent fraction of porosity from which anions are excluded, SAT 1y is saturated water content in the soil layer (mm H 2 O). To  www.nature.com/scientificreports/ obtain the transport of nitrate with runoff, percolation, and lateral flow, the following generic equation has been adopted from 55 . where NO 3 − represents the NO 3 − removed by the physical transport (i.e., lateral flow, runoff, percolation) (kg ha −1 ); β NO3 is the concentration of NO 3 − in the mobile water for the top 10 mm of soil/kg ha −1 considering both surface and subsurface lateral flow in the top layer; and Q x is the water physical transport (i.e. Q surf , Q lat , Q perc ). The NO 3 − , that percolates to the shallow aquifer from the soil profile may remain in the aquifer, or it may move with groundwater flow into the main river channel or into the deep aquifer. Organic transport of N with sediment is obtained as a concentration function proposed by 56 and later applied by 57 to the separated runoff events. Estimating the daily organic N runoff loss is on the basis of concentration of organic N in the top soil layer, the sediment yield and the enrichment ratio of the organic nitrogen in sediment to organic N in soil layer 51 see Fig. 2). In the SWAT model, water quality procedures integrate essential interactions and relationships used in the QUAL2E model 58 , which includes the major interactive factors such as the nutrient cycles, benthic oxygen demand, and algae production.
TDS is all the dissolved chemicals in the water (mg/L). It can be obtained by adding all the concentrations from chemical analysis or can be measured as weight of the residue after a volume of water has been evaporated to dryness. Typically, the concentration of TDS increases with the sample depth that groundwater has traveled from the recharge area to the sample sites. High nutrient concentration in the groundwater shows contamination form different sources of pollution 59 . Groundwater module in the SWAT model. Existing SWAT model divides underground saturated aquifer into shallow and deep aquifer. The balance of water in the shallow aquifer can be described 47,60 : where, S i sh and S i−1 sh is the volume of water stored in the shallow aquifer on the day i and i − 1, respectively, Q i gw,sh is refers to the groundwater from the shallow aquifer on the day i, W i rg,sh shows the volume of recharge which is entering the shallow aquifer on the day i, W i rp is indicates the amount of water entering to the soil for the evaporation as well as transpiration on the day i, W i pump,sh is refers to the amount of water removed from the shallow aquifer in the form of pumping on the day i. Assuming that the variation in groundwater is linearly correlated to the water table, the flow of groundwater can be represented by the following equations: www.nature.com/scientificreports/ where, α gw,sh refers to the groundwater recession constant of the shallow aquifer. Groundwater from the shallow aquifer contributed to the river on day i is obtained as: where S shthr,q is a threshold value above which the stored groundwater flows to river channels, α gw is the base flow recession constant, W rg,sh is the amount of recharge entering the shallow aquifer on day i, and t is the time step. The volume of recharge entering the shallow aquifer on day i is obtained using the following equations: where W i seep is the total amount of water exiting the bottom of the soil profile on day i, and δ gw is the drainage time of the overlying geologic formations. W i rg and W i−1 rg are the amount of recharge entering the aquifers on day i and i-1, respectively. δ gw is estimated against observed data in water table level through simulating aquifer recharge.
Likewise, from the deep aquifer groundwater, which is contributed to the stream on the day, i is obtained by: where Q i gw,dp and Q i−1 gw,dp refer to flow of groundwater from the deep aquifer into the stream on day i and i-1, respectively, and α gw,dp is the groundwater recession constant of the deep aquifer.
Due to complexity of groundwater, we separate the shallow aquifer into a lower aquifer and an upper aquifer in the groundwater module of SWAT model to improve the model accuracy, similar to Shao et al. 58 . Thus, the groundwater flow in Eq. (6) can be replaced using upper and lower aquifers as follows where, Q i gw,l and Q i gw,u refers to flow of groundwater from lower and upper aquifer into the stream on the day i respectively; where as −α gw,u and −α gw,l shows groundwater recession constant of the lower and upper aquifer respectively; W i rg,l and W i rg,u shows the volume of recharge entering into the lower and upper aquifer in the day i respectively.

SWAT model setup.
Delineation processes were done in ArcGIS. Based on soil types and landuse classes, the hydrological response unit (HRU) was defined, based on similar soil, similar landuse, and slope types. About 1370 HRUs were identified for the basins. Weather data, such as, temperature, rainfall, wind speed, humidity, and radiation observation were obtained from Alberta Environment and Parks. Groundwater well and water quality data (nitrogen species and TDS) at 1,300 different monitoring stations was obtained from Alberta Environment and Parks as shown in Fig. 1. Data availability is one of the most prominent factors which affect the model accuracy. Therefore, there were two monitoring stations with > 15 samples for the selected parameters between 2004 and 2016 namely IOR-KRL-03 (ss) and IOR-KRL-04 (ss) monitoring stations. Both groundwater monitoring wells are completed in surficial sand aquifers at 3 m and 11 m respectively.
Model performance metrics and uncertainty analysis. Measured groundwater quality data for NO 3 − and TDS were used for model calibration and validation. To evaluate model performance, coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE), and percent bias (PBIAS) were used 61 . The descriptions of NSE, R 2 , and PBIAS can be found in 62 as follows: gw,dp = Q i−1 gw,dp e −α gw,dp + W i rg,dp (1 − e −α gw,dp )  (2013)(2014)(2015). In order to perform parameter sensitivity and uncertainty analysis, we used Sequential Uncertainty Fitting2 (SUFI-2) algorithm in the SWAT-CUP 63 . The program generated various parameter sets for these selected parameters from a specified range of values using the Latin-Hypercube sampling technique (Table 1). To identify the total predictive uncertainty band of the simulated results, the SWATCUP was run several times, each one identifying a narrow parameter range for each parameter listed in Table 1, up to the point where reasonable goodness-of-fit statistics values were attained.

Results
Sensitivity and uncertainty analysis. Sensitivity test analysis is to identify the most sensitive parameters that govern NO 3 − and TDS in groundwater ( Table 1). The range of selected sensitive parameters are presented in Table 1. Some parameters are basin scale and others are basin-based parameters. The analysis of model sensitivity can be processed to find out the relative response of the SWAT model to the changes in relative value of specific model parameters. Hence, some parameters are sensitive to control the whole system processes as the most significant parameters. Here, the most sensitive parameters governing the groundwater parameters were assessed on the basis of the values obtained during primary model calibration. The CDN.bs, r__GWNO 3 .gw, HLIFE_NGW.gw, and NPERCO.bsn were ranked as the most sensitive parameters among the NO 3 − parameters. While NSETLR1.1wq and N_UPDIS.bsn were the sensitive parameters from the N parameters, that influenced the concentration. On the other hand, the parameters listed in Table 1 were the most sensitive for TDS, groundwater conditions, and surface runoff. The results of the sensitive analysis confirmed that measured input parameters have a substantial influence on the model prediction.
Model Calibration and Validation performance. Model calibration and validation are done after the sensitivity analysis. To compare the model performance, the simulated outputs were compared with observed value. Therefore, the daily NO 3 − and TDS from the groundwater monitoring stations IOR-KRL-03 (ss) and  Table 2 summarizes the performance statistics of the model for the daily nutrient concentration simulations for two groundwater monitoring stations. Table 2 shows satisfactory to very good for both stations with an averaged R 2 of 0.74 for nitrate during calibration and 0.71 during validation. This confirmed that the model was able to capture the concentration of nutrients after model modification (Table 2 and Fig. 3). Therefore, the overall model performance of the new SWAT module for the daily nutrient concentration simulations was in acceptable range of model calibration and validation in the ARB. In contrast, a lower model performance for TDS simulation was observed at both monitoring stations, for which the value of NSE is found to be 0.58 during validation at IOR-KRL-03 (ss) and 0.48 during calibration at IOR-KRL-04 (ss).

Nitrate (NO 3 − ).
Nitrate calibration processes happened in a nitrogen percolation coefficient of 0.5, with a denitrification threshold of water content and exponential rate of denitrification coefficient of 2.5 (Table 1).   (Fig. 3). Mineralization was adjusted by minimizing the default values of the rate factor for CMN to 0.000131. To control the depth distribution of nitrogen uptake, and to slow down simulated kinetics, the N-UPDIS was increased from the default value of 20 to 28 and therefore a large mass of NO 3 − was removed from the upper layers as per the report by 64,65 . The range of nitrogen settling in reservoirs was kept constant during the year as per 66 , the range was set greater than the default values from the Sava River Basin, to efficiently simulate the substantial retention in major wetlands because wetland specific retention is not applied in the extended SWAT model.
The nitrification rate decreased from the upper to lower Athabasca River Basin, which mirrors the rainfall distribution, with lower rainfall leading to lower soil saturation and lower nitrification rates. The maximum values of nitrate in groundwater of the ARB were 0.62 mg/L at IOR-KRL-03(ss) and 0.1 mg/L at IOR-KRL-04(ss) stations respectively, whereas the minimum values were 0.02 mg/L at the IOR-KRL-03(ss) station and 0.01 mg/L at the IOR-KRL-04(ss) station. The nitrification calibration resulted to satisfactory predictions of NO 3 − daily concentration in the two observation stations in the calibration and validation dataset. Ae per the comparison of NO 3 − between simulated and observed, better model performances were obtained in the model evaluation dataset. The simulated daily NO 3 − concentration agreed with the observed data and were acceptable range (Fig. 3). The percentage BIAS obtained from the observed and calculated daily loads also was ranked as acceptable to very good for all the selected stations ( Fig. 3 and Table 2).
Total dissolved solids (TDS). The TDS module for groundwater quality used the TDS concentrations from the selected groundwater monitoring stations in the ARB. The calibrated TDS parameters with specific ranges were presented in Table 1. Daily TDS concentrations were calibrated by adjusting the parameters related to groundwater. The model performance evaluation criteria reported by 62 were used for the daily nutrient simulation as a guideline in evaluating the model performance for the daily TDS concentrations. Figure 4 shows that for the whole period of simulation, R 2 , NSE and PBIAS values were found to be good to very good during calibration while they were satisfactory to good during the validation (Table 2). Generally, the overall model performance for the daily TDS concentration simulations in groundwater of the ARB shows that the model could capture the observed concentrations. When comparing the calculated and observed concentrations (Fig. 4), the simulated TDS concentrations were also acceptable as both show similar trends. Yet, local inconsistencies were noticed in the river basin. The highest percentage of overestimation and underestimation of TDS in the www.nature.com/scientificreports/ calibration dataset occurs for high concentration observations and probably reflects the SWAT model's representation of high level concentrations 67 , which may cause errors in the process of calculating TDS in the SWAT model 68 . The other possible sources of overestimation and underestimation of the model simulations originate from uncertainties in the input data and observed data. For example, the frequency of data collection varied from once or twice per month to once every few years. In contrast, the simulations used a daily time step. The low frequency of sampling might miss the peak or valley of nutrient concentrations. Total Dissolved Solids is an indicator of the availability of total dissolved salt and other constituents that affect groundwater quality. High TDS concentrations in groundwater can occur naturally, and may also indicate a downward movement of leachate into groundwater 65,[69][70][71] . The average TDS concentrations in groundwater in the study area are 143 mg/L and 1421 mg/L at IOR-KRL-03(ss) and IOR-KRL-03(ss) stations, respectively. The TDS of groundwater samples ranges from 120-182 mg/L at the IOR-KRL-03(ss) station, while the TDS ranges from 1380-1500 mg/L at the IOR-KRL-04(ss) station (Fig. 4). High concentrations of TDS in groundwater reduce the palpability of water for drinking and may cause gastrointestinal pain and emetic effects in humans 67 . TDS is an important indicator for evaluating the quality of groundwater. High levels of TDS typically occur for hard water and might require groundwater treatment to decrease concentrations below 500 mg/L 67 . The performance of model results were investigated by using the analysis of scatter plot between observed and simulated estimates at each groundwater monitoring stations (Fig. 5). The slope at the selected stations were considerably far from zero. This showed that the model prediction accuracy was enhanced by extended groundwater model. Thus, the results of the 1:1 fitting line confirmed that the extended groundwater SWAT module was effective prediction of nutrient concentration in the groundwater. The scatters are closer to the 1:1 line for the entire study time at two monitoring stations, albeit at some points relatively far from the fitting line probably due to limited availability of observed data. In general, it is worth to conclude the extended SWAT groundwater module shows better efficiency for simulating nutrient concentrations in the groundwater.
Groundwater quality analysis and nutrient concentrations. Analysis of observed against simulated nitrate and TDS concentrations in ARB groundwater helps to identify sources of errors across two groundwater monitoring stations within the river basin. Observed nitrate concentrations vary significantly across the area (Fig. 3). Concentrations of NO 3 − in groundwater have been highly affected by emissions of both point and nonpoint sources in different watersheds across the world 71 . To assess pollution sources and quantify the loads of NO 3 − entering the whole river basin, NO 3 − load to groundwater has been simulated at various hydrogeological On the other hand, at station at IORKRL-04(ss) the longterm daily average NO 3 − concentrations recorded a maximum value of 0.1 mg/L during the year 2015, and a minimum value of 0.013 mg/L in most observation years. Generally, the concentrations of nitrate were well captured in all the groundwater monitoring stations, although some overestimations were observed at some points at IOR-KRL-04(ss) monitoring station (Fig. 3) and some underestimations were found at IOR-KRL-03(ss) and IOR-KRL-04(ss) monitoring stations. The probable reason for model overestimation is the data uncertainty such rainfall and snow distribution, in which the rainfall and snow distribution caused in nutrient concentrations. The Underestimation was probably subjected to the uncertainty of the input data and measured data. Furthermore, the observed NO 3 − data is not sufficient for continuous daily NO 3 − concentration for the entirely considered period for this study. Therefore, observed data could be sources for the errors. The simulation of NO 3 − removal by soil denitrification could be improved by varying the river basin parameters. However, the existing SWAT model could not represent perfectly the seasonal difference in nitrate concentrations. After carefully extending the model, the findings of this study highlight the need to improve spatial representation of nitrate concentrations in groundwater and in the parameters that influence nitrate concentrations at the watershed scale.

Concluding remarks
Groundwater is a precious natural resource supporting the existence of life on earth. However, various factors, such as soil properties, crop growth, industrial wastes, and agricultural management, can influence its quality. Development of industries and agriculture increase water use, which puts pressure on groundwater quality and may, in turn, influence the ecosystem of the Athabasca River Basin. Particularly, the use of fertilizers, manure, and industrial wastes may contribute to the pollution of groundwater and connected surface water environments. It is necessary to perform optimal management of groundwater resources in the basin.
State-of-the-art groundwater quality modelling at the ARB is recognized as a vital component of groundwater management, where further in-depth studies may be required to offer valuable insights related to groundwater condition and nutrient processes at the various spatio-temporal scales (i.e. site and region, and daily, monthly and annually). This can better support nutrient monitoring networks for river basin management, and enhance understanding of changing nutrient concentrations in the ARB. The SWAT model results could support the development of indicators for groundwater quality parameters, and also support integrated surface water and groundwater management. However, the accuracy of SWAT predictions is limited by data availability and structure of the model. This may result in errors while modeling processes in the river basin. For example, the highest percentage of overestimation and underestimation for TDS in the calibration dataset probably reflected the poor representation of SWAT model to high concentrations, as well as uncertainties of the input data and observed data, which in turn may cause errors in estimating groundwater TDS in the SWAT model. With continued monitoring of the groundwater network at regular frequencies, high quality input data would be available for better calibration and validation of the model. Furthermore, additional data would improve the representation of processes and pathways controlling groundwater pollution and therefore allow evaluation of the effect of land and water management practices to support the implementation of the best management practices.
In this study, we extended the existing SWAT model to improve modelling of groundwater quality (i.e. NO 3 − and TDS). A systematic calibration and validation of the SWAT model has been performed to compare the observed groundwater NO 3 − and TDS concentrations in the ARB. The results revealed that the new groundwater quality model in the SWAT is able to capture the daily nutrient concentrations in groundwater. The simulated results agree with the observed data with satisfactory and good performance for both groundwater monitoring stations as per the model performance measures. Thus, the process-based hydrologic groundwater quality model is an effective tool in simulating the groundwater quality dynamics (NO 3 − and TDS) for sustainable groundwater and surface water management in the river basin.