Application of an Integrated SWAT – MODFLOW Model to Evaluate Potential Impacts of Climate Change and Water Withdrawals on Groundwater – Surface Water Interactions in West-Central Alberta

It has become imperative that surface and groundwater resources be managed as a holistic system. This study applies a coupled groundwater–surface water (GW–SW) model, SWAT–MODFLOW, to study the hydrogeological conditions and the potential impacts of climate change and groundwater withdrawals on GW–SW interactions at a regional scale in western Canada. Model components were calibrated and validated using monthly river flow and hydraulic head data for the 1986–2007 period. Downscaled climate projections from five General Circulation Models (GCMs), under the RCP 8.5, for the 2010–2034 period, were incorporated into the calibrated model. The results demonstrated that GW–SW exchange in the upstream areas had the most pronounced fluctuation between the wet and dry months under historical conditions. While climate change was revealed to have a negligible impact in the GW–SW exchange pattern for the 2010–2034 period, the addition of pumping 21 wells at a rate of 4680 m3/d per well to support hypothetical high-volume water use by the energy sector significantly impacted the exchange pattern. The results showed that the total average discharge into the rivers was only slightly reduced from 1294 m3/d to 1174 m3/d; however, localized flowrate differences varied from under 5 m3/d to over 3000 m3/d in 320 of the 405 river cells. The combined potential impact is that intensive groundwater use may have more immediate effects on river flow than those of climate change, which has important implications for water resources management and for energy supply in the future.


Introduction
Proper water management is a key global challenge.Indeed, there are already indicators that anthropogenic influences are impacting the flow rates of rivers and the winter recovery of glaciers [1].The energy, agricultural and domestic sectors are found worldwide, and even as relative water withdrawals in these sectors change, the overall global volume of extracted water is projected to increase.This becomes especially clear as factors such as population growth and improvements to the quality of life come into play [2].When the water needs of each sector are combined, it becomes clear that water must be managed systematically in order to ensure a sustainable supply for the future.Therefore, understanding how both natural and anthropogenic factors affect a watershed is key to the proper management of this resource.

Surface Hydrology
The Little Smoky watershed of Alberta (Figure 1) covers 11,494 km 2 and is a tributary of the Peace River, which eventually flows into the Mackenzie River and the Arctic Ocean.The hydrology of the Little Smoky watershed is heavily influenced by the seasons.As it is principally fed by the Rocky Mountains, more discharge is recorded in the river in the spring and summer months, which is the result of increased precipitation and the yearly spring snowmelt.Similarly, a pronounced decrease in discharge is consistently recorded in the autumn and winter months, coinciding with freezing and lower precipitation levels.Like much of Alberta, the area is subjected to periodic droughts and floods [25], which can heavily influence local hydrologic flow patterns and, as a consequence, water use.The downstream portion of the Little Smoky River watershed flattens into forest and prairie land, much of which is dominated by agricultural and upstream oil and gas activity, with the forestry industry also having a presence.The location of the Little Smoky River in west-central Alberta means that it experiences a generally semi-arid climate, and the watershed also incorporates a wide range of elevations (from ~475 m to ~1400 m above sea level) which brings with it several soil and land use types.

Geology and Hydrogeology
Given the size of the study area, the depths and properties of the underlying rock formations vary considerably.In the study area, similar to the rest of Alberta's Rocky Mountains, various formations outcrop or come close to outcropping, with a general trend toward the northeast [26].Quaternary-Neogene sediments cover the majority of the southern near-surface area of the study region, with various members of the Paskapoo Formation directly underlying it.However, in the downstream (northern) area of the Little Smoky River watershed, older formations outcrop, such as the Scollard, Battle and Wapiti Formations.Each of the bedrock formations has been evaluated for aquifer potential from the perspective of the entire WCSB [27]; however, locally the Paskapoo Formation has largely been the focus of groundwater source wells.As the Paskapoo Formation is also the shallowest bedrock unit in the area of study, it has exchange between SW and GW [28].Deposited in the Paleogene, the Paskapoo Formation is a highly heterogeneous fluvial deposit [22] that appears as the shallowest bedrock layer in most of the Little Smoky River watershed.Even within the study area, this formation exhibits a high range in thickness, following a thickening trend toward the foothills [29].Three primary units, or members, within the formation have been identified in the study area based on lithological differences, with two aquifers separated by an aquitard [22,30].The porosity and hydraulic conductivity observed by Hughes et al. [22] in these members share a linear relationship, with higher values for both porosity and hydraulic conductivity recorded in the aquifers.In the study completed by Hughes et al. [22], average conductivities were 1.3 × 10 −6 m/s and 4.4 × 10 −9 m/s for the sandstone and mudstone/shale units respectively, primarily obtained via air permeameter testing.

SWAT Model Development
ArcSWAT 2012 was used to build a SWAT model of the study area.A Digital Elevation Model (DEM) [31] was then imported, and used to delineate the watershed into 28 sub-basins (Figure 2) and to replicate the water network of the study area (primarily the Little Smoky River).The model calculates the direction and accumulation of surface water flow based on low points in elevation and represents, in each sub-basin, a tributary stream that eventually joins into the Little Smoky River.Upon defining the main outlet (or outlets) of the watershed, the boundaries and hydrologic parameters of the entire basin are then delineated.For this study, the primary outlet of the watershed is located at the northernmost point of the modelled region.

Geology and Hydrogeology
Given the size of the study area, the depths and properties of the underlying rock formations vary considerably.In the study area, similar to the rest of Alberta's Rocky Mountains, various formations outcrop or come close to outcropping, with a general trend toward the northeast [26].Quaternary-Neogene sediments cover the majority of the southern near-surface area of the study region, with various members of the Paskapoo Formation directly underlying it.However, in the downstream (northern) area of the Little Smoky River watershed, older formations outcrop, such as the Scollard, Battle and Wapiti Formations.Each of the bedrock formations has been evaluated for aquifer potential from the perspective of the entire WCSB [27]; however, locally the Paskapoo Formation has largely been the focus of groundwater source wells.As the Paskapoo Formation is also the shallowest bedrock unit in the area of study, it has exchange between SW and GW [28].Deposited in the Paleogene, the Paskapoo Formation is a highly heterogeneous fluvial deposit [22] that appears as the shallowest bedrock layer in most of the Little Smoky River watershed.Even within the study area, this formation exhibits a high range in thickness, following a thickening trend toward the foothills [29].Three primary units, or members, within the formation have been identified in the study area based on lithological differences, with two aquifers separated by an aquitard [22,30].The porosity and hydraulic conductivity observed by Hughes et al. [22] in these members share a linear relationship, with higher values for both porosity and hydraulic conductivity recorded in the aquifers.In the study completed by Hughes et al. [22], average conductivities were 1.3 × 10 −6 m/s and 4.4 × 10 −9 m/s for the sandstone and mudstone/shale units respectively, primarily obtained via air permeameter testing.

SWAT Model Development
ArcSWAT 2012 was used to build a SWAT model of the study area.A Digital Elevation Model (DEM) [31] was then imported, and used to delineate the watershed into 28 sub-basins (Figure 2) and to replicate the water network of the study area (primarily the Little Smoky River).The model calculates the direction and accumulation of surface water flow based on low points in elevation and represents, in each sub-basin, a tributary stream that eventually joins into the Little Smoky River.Upon defining the main outlet (or outlets) of the watershed, the boundaries and hydrologic parameters of the entire basin are then delineated.For this study, the primary outlet of the watershed is located at the northernmost point of the modelled region.Soil and land use data for the study were imported via component-based GIS layers so that single polygons could host many soil and land use types [32].Based on available soil and land use maps (see Table A2), 20 types of soil and 13 land use types were defined in the area.The various land use types of the study area correspond loosely with the elevations at which they are found, with southern sub-basins dominated by evergreen forests, and agricultural activity being more pronounced as ground elevations level out in the northern portion of the watershed.After data for soil and land use were included, we chose a multiple-slope model to better represent the topographic variation found in the field.We opted for three slope classes: 0-5% incline, 5-10% incline, and over 10% incline.With the streams, sub-basins, soils, land use and slopes in place, Hydrological Response Units (HRUs), which are sub-sections of the model that further break the spatial units up into pockets with similar properties to represent local soil, land use, slope characteristics of a watershed [4], could then be delineated.Since they are defined in this manner, HRUs may not comply with the boundaries between sub-basins.Instead, SWAT calculations are performed for each HRU, and then the results are aggregated to represent values at sub-basin outlets, thus producing outputs for each sub-basin and every time step.
Historical weather data, including precipitation, temperature, relative humidity, wind speed, and solar radiation, were collected from different sources (Table A2), and were incorporated into the model for the time period of 1983-2007 (25 total run years).In the SWAT model, snowfall is simulated based on total precipitation and temperature data.

MODFLOW Model Development
The MODFLOW model was constructed using the Visual MODFLOW Classic Interface.To allow for the drying and re-wetting of grid cells over a transient simulation, the MODFLOW-NWT (Newtonian Solver) engine was used [33], similar to the models used by Guzman et al. [12] and Bailey et al. [4].The model area for this component was 44,955 km 2 , with a length of 195 km in the x-direction and 230 km in the y-direction.Discretization for this model was intentionally kept coarse to ease the computational burden once coupled with SWAT.A 100 × 100 grid was used, with each grid cell having dimensions of 1.95 km × 2.30 km.While this does result in relatively low spatial resolution in the x-y directions, the MODFLOW component remains relatively complex compared to those used in previous SWAT-MODFLOW models by virtue of its seven layers, and the heterogeneous Soil and land use data for the study were imported via component-based GIS layers so that single polygons could host many soil and land use types [32].Based on available soil and land use maps (see Table A2), 20 types of soil and 13 land use types were defined in the area.The various land use types of the study area correspond loosely with the elevations at which they are found, with southern sub-basins dominated by evergreen forests, and agricultural activity being more pronounced as ground elevations level out in the northern portion of the watershed.After data for soil and land use were included, we chose a multiple-slope model to better represent the topographic variation found in the field.We opted for three slope classes: 0-5% incline, 5-10% incline, and over 10% incline.With the streams, sub-basins, soils, land use and slopes in place, Hydrological Response Units (HRUs), which are sub-sections of the model that further break the spatial units up into pockets with similar properties to represent local soil, land use, slope characteristics of a watershed [4], could then be delineated.Since they are defined in this manner, HRUs may not comply with the boundaries between sub-basins.Instead, SWAT calculations are performed for each HRU, and then the results are aggregated to represent values at sub-basin outlets, thus producing outputs for each sub-basin and every time step.
Historical weather data, including precipitation, temperature, relative humidity, wind speed, and solar radiation, were collected from different sources (Table A2), and were incorporated into the model for the time period of 1983-2007 (25 total run years).In the SWAT model, snowfall is simulated based on total precipitation and temperature data.

MODFLOW Model Development
The MODFLOW model was constructed using the Visual MODFLOW Classic Interface.To allow for the drying and re-wetting of grid cells over a transient simulation, the MODFLOW-NWT (Newtonian Solver) engine was used [33], similar to the models used by Guzman et al. [12] and Bailey et al. [4].The model area for this component was 44,955 km 2 , with a length of 195 km in the x-direction and 230 km in the y-direction.Discretization for this model was intentionally kept coarse to ease the computational burden once coupled with SWAT.A 100 × 100 grid was used, with each grid cell having dimensions of 1.95 km × 2.30 km.While this does result in relatively low spatial resolution in the x-y directions, the MODFLOW component remains relatively complex compared to those used in previous SWAT-MODFLOW models by virtue of its seven layers, and the heterogeneous hydraulic conductivity found within these layers (thus maintaining variability in soil hydraulic conductivity (K) in all three dimensions).
Each of the MODFLOW layers correspond to a geological unit or sub-unit with the formation elevation data (in meters above sea level) determined from unpublished data from the Alberta Geological Survey and additional well log information from IHS AccuMap ® software.Although the formation elevations are highly variable, a uniform model bottom of 1250 m below sea level was included to keep the lower boundary of the model well away from surface water features, and to allow for possible water flow in areas where deeper bedrock formations outcrop to the surface at the northern end of the model (see Figure 3).
Water 2019, 11, x FOR PEER REVIEW 6 of 30 hydraulic conductivity found within these layers (thus maintaining variability in soil hydraulic conductivity (K) in all three dimensions).Each of the MODFLOW layers correspond to a geological unit or sub-unit with the formation elevation data (in meters above sea level) determined from unpublished data from the Alberta Geological Survey and additional well log information from IHS AccuMap ® software.Although the formation elevations are highly variable, a uniform model bottom of 1250 m below sea level was included to keep the lower boundary of the model well away from surface water features, and to allow for possible water flow in areas where deeper bedrock formations outcrop to the surface at the northern end of the model (see Figure 3).Hydraulic conductivity ranges for each unit are based on values determined for the study area [22], and estimates based on known lithology.Each layer in the model has multiple conductivity values, which is supported by published values regarding the hydrostratigraphy of formations in the study region [22].The table of the MODFLOW input data can be viewed in the Appendix A (Table A2).
Three types of boundary conditions were imposed upon the MODFLOW model: constant heads (CHD), rivers (RIV), and recharge (RCH).The CHD data used was implemented based on observed hydraulic heads in nearby observation or legacy wells.In order to better simulate the variation in observed heads, linear gradients were used for the majority of CHD inputs, with hydraulic heads increasing or decreasing in a linear fashion.In addition, heads in areas with lower surface water elevation (northern/downstream portion of the model) were assumed to be at or near surface level for near-surface formations such as the Paskapoo.
RIV and RCH data were both taken directly from the corresponding SWAT model, so that the values remained constant once the two models were coupled.The complete process of transferring SWAT river data to MODFLOW is given in the Appendix A (Workflow 1).Further detail on the MODFLOW model construction can be found in Appendix B.

SWAT-MODFLOW Integration
To integrate SWAT and MODFLOW, SWAT HRUs must be disaggregated, after which they are called DHRUs.This separates a given HRU into polygons, which allows for it to be geo-located, facilitating the connection between SWAT and MODFLOW [34].DHRUs pass location and SW flow data to MODFLOW grid cells on a daily time step, which is then returned to SWAT as calculated hydraulic heads and GW-SW interaction flow rates [4].For the 40 HRUs in this model, there are 31047 DHRUs.After the disaggregation of the HRUs, linking files must be created that tell the SWAT-MODFLOW code how many HRUs, DHRUs and MODFLOW grid cells exist in the project.Only then can the DHRUs be linked to their corresponding MODFLOW grid cells for data to be transferred.The entire process of creating/exporting the correct linking files and coupling the SWAT-MODFLOW Hydraulic conductivity ranges for each unit are based on values determined for the study area [22], and estimates based on known lithology.Each layer in the model has multiple conductivity values, which is supported by published values regarding the hydrostratigraphy of formations in the study region [22].The table of the MODFLOW input data can be viewed in the Appendix A (Table A2).
Three types of boundary conditions were imposed upon the MODFLOW model: constant heads (CHD), rivers (RIV), and recharge (RCH).The CHD data used was implemented based on observed hydraulic heads in nearby observation or legacy wells.In order to better simulate the variation in observed heads, linear gradients were used for the majority of CHD inputs, with hydraulic heads increasing or decreasing in a linear fashion.In addition, heads in areas with lower surface water elevation (northern/downstream portion of the model) were assumed to be at or near surface level for near-surface formations such as the Paskapoo.
RIV and RCH data were both taken directly from the corresponding SWAT model, so that the values remained constant once the two models were coupled.The complete process of transferring SWAT river data to MODFLOW is given in the Appendix A (Workflow 1).Further detail on the MODFLOW model construction can be found in Appendix B.

SWAT-MODFLOW Integration
To integrate SWAT and MODFLOW, SWAT HRUs must be disaggregated, after which they are called DHRUs.This separates a given HRU into polygons, which allows for it to be geo-located, facilitating the connection between SWAT and MODFLOW [34].DHRUs pass location and SW flow data to MODFLOW grid cells on a daily time step, which is then returned to SWAT as calculated hydraulic heads and GW-SW interaction flow rates [4].For the 40 HRUs in this model, there are 31047 DHRUs.After the disaggregation of the HRUs, linking files must be created that tell the SWAT-MODFLOW code how many HRUs, DHRUs and MODFLOW grid cells exist in the project.Only then can the DHRUs be linked to their corresponding MODFLOW grid cells for data to be transferred.
The entire process of creating/exporting the correct linking files and coupling the SWAT-MODFLOW model is explained in the tutorial made available by Park and Bailey [34].Further information on executing the SWAT-MODFLOW model can be found in Appendix B.

SWAT Calibration
Calibration guidelines are described within the documentation for the SWAT-CUP software [35] and were implemented in this study.Within the SWAT software, the model runs successfully for the simulated duration of 25 years.For the optimal run, a warm-up period of three years was used, bringing the first data year to 1986.The model was then calibrated by using observed monthly stream flow data at two hydrometric stations within the study area (Outlet Gauge and Tributary Gauge, Figure 3).In this study, we used the Sequential Uncertainty Fitting program (SUFI2) [3].We performed sensitivity analyses to determine which parameters most affect the streamflow simulation.The calibration procedure starts by providing a large range that is defined by user.However, this initial range is physically meaningful, and is based on input data that is initially obtained from different sources.For example, soil-related parameters such as soil hydraulic conductivity (K) and soil available water capacity (AWC) are attributed to the soil map, and they are initially obtained by taking soil samples and analysis of soil profiles in separate studies by the data providers.These parameters, that are called default parameters, provide a base for the modeler to make decisions about the initial parameter range.The parameter range narrows down through an iterative procedure, until the best performing parameter range is obtained.In each iteration, the SUFI2 algorithm performs Latin Hypercube sampling for user defined parameter ranges and creates multiple parameter set samples.Each parameter set sample is fed into the SWAT model, and the possible outputs are simulated within a 95% confidence interval (95PPU).The SUFI2 program compares the simulated discharge to the observed data using a variety of statistics, including R 2 and bR 2 .The R 2 , or coefficient of determination, relates how well the simulated data trend replicates the observed trend.The primary factor observed in this study, however, is bR 2 , which takes both trend and closeness to observed results into account, and is calculated by multiplying R 2 by b, or the slope of the regression line between observed and simulated data [36].Ranging from 0 to 1, a higher value indicates a better-calibrated model.However, since it is impossible to reach a bR 2 value of 1 other than by over-calibrating, values are considered acceptable once they start to plateau after many iterations (in the range of 0.5-0.7).Also important are the p-factor and r-factor of the results.The p-factor of an iteration, measured on a scale of 0 to 1, considers how well the best simulation stays within the 95% confidence interval that is due to the parameters' uncertainty range.As such, a p-factor closer to 1 is ideal.The r-factor represents the range of simulated uncertainty.As a well-calibrated model should simulate data within a range close to the observed results, a smaller r-factor is ideal, with a result under 1 considered to be passable [3,37].The best performing parameter set in each iteration is considered to guide the user to narrow the parameter range in subsequent iterations until a desirable performance of p-factor and r-factor is obtained.
Within the study area, two hydrometric stations are present, which are situated in sub-basins 1 (station 119_1) and 24 (station 118_24) and are used as calibration points (Figure 2).These hydrometric stations yield river discharge data on a monthly basis, which are used to calibrate against the simulated river discharge calculated in SWAT at the outlet of previously mentioned sub-basins.In total, we evaluated 45 parameters of the SWAT model (see Table A4).Due to the proximity of sub-basin 24 to the Rocky Mountains, it was treated as its own system, without further regionalization of the parameters for other sub-basins.Snow-related parameters were also added as a variable, due to their potential impact on the hydrologic budget of the model (see Table A4).Between iterations, widening and narrowing ranges of certain parameters within a physically meaningful range resulted in changes to the accuracy of calibration.The parameter ranges used within this study were based on those in the SWAT input-output documentation [38], as well as those of the soils database for the project (see Table A2 for data source).When the simulated values for river discharges are close to matching the observed values for each stream gauge (yielding bR 2 , p-factor and r-factor values within the desirable range discussed above) for the entire simulation period, the SWAT model is considered to be well-calibrated.As more recent historical data are generally accepted to be more reliable (and was of greater interest due to the planned forecasting scenarios), the SWAT model was calibrated for the 1996-2007 period, with model validation occurring for the 1986-1995 period.

MODFLOW Calibration
The MODFLOW model was calibrated using measured GW water level data from 377 wells available from the Alberta Water Well Information Database and two wells from the Alberta Groundwater Observation Network (GOWN) that have time series data (Figure 4).The project setup and initial runs were carried out under steady-state (SS) conditions, to establish an initial condition for subsequent runs.This was done so that the initial heads file could be used to run transient simulations.After this initial run, variable recharge data were added to the model, with the corresponding SWAT model as its source.These recharge data were subdivided by SWAT sub-basin, and were added to the top layer of the MODFLOW model in the corresponding locations (see Figure 3).As expected, the computational time increased as more detailed recharge data were introduced to the model.When monthly values were introduced for each sub-basin and run under transient conditions, the total run-time was just over thirty minutes.
The water levels of each calibration well were converted to hydraulic head data by subtracting them from the corresponding ground elevations from the DEM, with screen elevations assumed to be at the bottom of each well.A total of 377 water wells were available within the MODFLOW modelled area, with some falling outside the Little Smoky River watershed to aid in representing the regional groundwater conditions.Due to the number of data points, calibration was limited to well measurements taken in the last five years of the simulation period (January 2003-December 2007).
Water 2019, 11, x FOR PEER REVIEW 8 of 30 the desirable range discussed above) for the entire simulation period, the SWAT model is considered to be well-calibrated.As more recent historical data are generally accepted to be more reliable (and was of greater interest due to the planned forecasting scenarios), the SWAT model was calibrated for the 1996-2007 period, with model validation occurring for the 1986-1995 period.

MODFLOW Calibration
The MODFLOW model was calibrated using measured GW water level data from 377 wells available from the Alberta Water Well Information Database and two wells from the Alberta Groundwater Observation Network (GOWN) that have time series data (Figure 4).The project setup and initial runs were carried out under steady-state (SS) conditions, to establish an initial condition for subsequent runs.This was done so that the initial heads file could be used to run transient simulations.After this initial run, variable recharge data were added to the model, with the corresponding SWAT model as its source.These recharge data were subdivided by SWAT sub-basin, and were added to the top layer of the MODFLOW model in the corresponding locations (see Figure 3).As expected, the computational time increased as more detailed recharge data were introduced to the model.When monthly values were introduced for each sub-basin and run under transient conditions, the total run-time was just over thirty minutes.
The water levels of each calibration well were converted to hydraulic head data by subtracting them from the corresponding ground elevations from the DEM, with screen elevations assumed to be at the bottom of each well.A total of 377 water wells were available within the MODFLOW modelled area, with some falling outside the Little Smoky River watershed to aid in representing the regional groundwater conditions.Due to the number of data points, calibration was limited to well measurements taken in the last five years of the simulation period (January 2003-December 2007).In addition to the observation wells mentioned above, the two river monitoring gauges within the study area (Outlet Gauge and Tributary Gauge) were also input as continuous calibration points (data for the entire simulation period), bringing the total to 379.Although the flow rates at each point are commonly used, monthly river stage was used instead for this component (provided by Alberta Environment and Parks; http://environment.alberta.ca/forecasting/FAQ/near_real_time.html),due to In addition to the observation wells mentioned above, the two river monitoring gauges within the study area (Outlet Gauge and Tributary Gauge) were also input as continuous calibration points (data for the entire simulation period), bringing the total to 379.Although the flow rates at each point are commonly used, monthly river stage was used instead for this component (provided by Alberta Environment and Parks; http://environment.alberta.ca/forecasting/FAQ/near_real_time.html),due to the MODFLOW primary calibration parameter being hydraulic head.Further information can be found in Appendix B.
At the onset of MODFLOW calibration, the thin, pinching-out layers in the northern end of the model (seen in Figure 3) caused problems for vertical flow.Since the simulated water could not flow down through the cells effectively, this resulted in areas where hydraulic heads would build to be higher than what can be observed (or is expected to be realistic).Introducing Constant Head boundary conditions (CHB) on the perimeter of the MODFLOW modelled area helped to bring the simulated hydraulic heads down to a more realistic range.These CHBs were based primarily on the ground elevation, and on the hydraulic head map of Brinsky [39] for deeper formations.Care was taken to not include CHBs in the model area that would be directly coupled with SWAT, allowing the simulated hydraulic head to fluctuate as needed.
Apart from imposing CHBs proximal to the study area, changing the hydraulic conductivity of each layer was another necessary factor in achieving a successful calibration.While remaining faithful to the realistic ranges outlined by literature [22,29], each MODFLOW layer was given multiple hydraulic conductivity values so as to simulate geological units that were heterogeneous in two dimensions (x and y).These layers could not be vertically heterogeneous due to each being a single grid cell thick, regardless of the actual thickness of the cells.Although calibration initially began with single-value recharge data for each sub-basin, this was subsequently changed to yearly, and finally monthly data.

Scenario Analysis
Two potentially significant contributors to hydrological change in this watershed are water well pumping for industrial use, and climate change.For this purpose, climate data simulated by five widely-used General Circulation Models (GCMs) were used (Table 1), which are climate change tools based on physical processes [40].Although GCMs simulate processes associated with climate and weather changes at a large spatial scale (some with resolutions of roughly 10 km), their output data must be downscaled for application to smaller-scale watersheds and local climates [41].A widely used method to accomplish this task is by establishing an empirical link between the selected GCM(s) and observed data for the historical period from the study area, known as the change factor approach [41].By using this method, the formulations obtained for the historical period can then be applied for future projections.The strength of GCMs lies in their ability to forecast scenarios with different assumptions for the future, known as Representative Concentration Pathways, or RCPs [42].Based on the projected emissions of greenhouse gases, the most widely used of these RCPs are RCP2.6 (the most environmentally optimistic scenario) and RCP8.5 (the heaviest greenhouse gas emission scenario).Each produces distinct model outputs and dedicated climate change studies use these to obtain the widest range of possible climate outcomes, especially when propagating climate change into groundwater models [43].
A key objective of this study was to evaluate the ability of SWAT-MODFLOW to simulate a watershed under the effects of both natural and anthropogenic influences.The results for three scenarios and a base scenario are presented, with the first focusing solely on the effects of pumping 21 simulated wells in the MODFLOW component.These wells were proximal to or within the study area (found along the same trends as the water wells used for calibration), pumping at a rate of 468 m 3 /d for six months of each simulation year.Although the durations of single energy operations are generally short-lived relative to the simulation period, the volumetric rate used here is consistent with one provided by an Albertan energy company with hydraulic stimulation operations near the study area.The second coupled scenario will include climate change and no pumping, averaging the results for five distinct GCMs (Can_ESM2, CCSM4, CNRM_CM5, CSIRO-MK3, and MIROC5) from the Pacific Climate Impacts Consortium (PCIC, https://www.pacificclimate.org/data)under the RCP8.5 Water 2019, 11, 110 10 of 28 (high carbon emissions) scenario (Table 1).The third applies both the climate change scenario described above and increased the pumping rate of the wells used in the first scenario by an order of magnitude to 4680 m 3 /d.The ten-fold increase in the pumping rate at the existing wells is meant to act as a way of simulating the installation of 10× as many wells, because the wells are well-distributed in the geographic region of the model where groundwater withdrawals for industry are most likely to occur.Although this scenario is unlikely to play out in the future (to this extent), it was included to maximize the stress on the hydrologic system to determine which extreme (climate change or heavy pumping) would have a greater influence on the system.These GCM models (Table 1) were downscaled to match the observed interval of the reference period  to establish an empirical relationship between them and the model.This relationship is then used to project data for future scenarios, resulting in a simulation period of 25 years (2010-2034), similar to the historical model.The pumping wells were primarily situated in the downstream portion of the MODFLOW component, where energy operations are more common.The boundary conditions and physical parameters of the MODFLOW model remained unchanged apart from the increased pumping rates, as the subsurface properties are assumed to be constant on a scale of decades.Furthermore, the precipitation is replaced by the new climate change-affected values of the SWAT model.

SWAT
With a preliminary SWAT run prior to calibration, the simulated results overestimated river peak flows when compared to the observed values (Figure 5a).For calibration purposes, a total of 21 physical parameters sensitive to river discharge were initially selected from literature [3].The parameters were further regionalized to capture a higher degree of spatial variability due to land use, soil, and hydrologic conditions for a total of 45 sensitive parameters.From these, a one-at-a-time, and then global sensitivity analysis was performed, evaluating model reactions to changes in single parameters and later to all parameter changes [45].The addition of elevation bands (PLAPS and TLAPS, allowing multiple temperature and precipitation values to capture the effects of orographic changes within a sub-basin) improved the overall calibration performance, particularly in Sub-basin 24, which was treated as its own system.Allowing for more water capacity in the soil (SOL_AWC), and less runoff (curve number; CN2) also decreased the overestimated peak flows.The final values were kept within the range provided by the soil input data, with a maximum water capacity of 0.333 mm H 2 O and a curve number decrease of 20%.To further constrain this peak flow reduction, the soil hydraulic conductivity values were increased as well.In addition to the above, snow-related parameters (.sno) were altered so that the snowfall and snowmelt simulations occurred in the correct temperature range (−5 to 5 • C).A complete list of the parameters used and their initial and final value ranges, can be viewed in the Appendix A (Figure A4).
Water 2019, 11, x FOR PEER REVIEW 11 of 30 mm H2O and a curve number decrease of 20%.To further constrain this peak flow reduction, the soil hydraulic conductivity values were increased as well.In addition to the above, snow-related parameters (.sno) were altered so that the snowfall and snowmelt simulations occurred in the correct temperature range (−5 to 5 °C ).A complete list of the parameters used and their initial and final value ranges, can be viewed in the Appendix A (Figure A4).In Figure 5b, the calibrated results can be seen for the outlet gauge.The data are shown on a monthly basis, and clearly outline the contrast between the peak flows of the summer and the low flows of the winter months.By changing the parameters as described above and in Section 2.3.1 and 2.3.2, the simulated overprediction of peak stream discharge was mitigated.Validation of the model was performed by comparing the simulation to observed data for the 1986-1995 period.Overall, the model follows very similar trends, with the peak flows of some years being slightly overpredicted, and the low flows of the winter months falling to zero for most simulation years.Apart from not being the primary focus of calibration, the discrepancies between the simulated and observed curves for the validated timespan of the SWAT model may be due to the short period of calibration/validation.In general, model calibration in both the tributary and watershed outlets ensured proper apportioning of precipitation and soil water into surface runoff, actual evapotranspiration, and groundwater recharge.This improved model performance compared to the pre-calibration model (Table 2, Figure 5a).Overall, 48% of the observed river flow data were captured by the simulated 95PPU and the average r-factor was about 0.75 at the watershed scale.The final summary statistics for the SWAT calibration and validation process at each stream gauge are given below (Table 2).The lower performance at the tributary gauge may be due to inadequate quality and quantity of input data, which is not representative at the study's level of detail.The spatial and temporal quality and quantity of geospatial maps and climate time series are often not able to perform ideally for each individual sub-basin (including upstream tributaries with small drainage In Figure 5b, the calibrated results can be seen for the outlet gauge.The data are shown on a monthly basis, and clearly outline the contrast between the peak flows of the summer and the low flows of the winter months.By changing the parameters as described above and in Sections 2.3.1 and 2.3.2, the simulated overprediction of peak stream discharge was mitigated. Validation of the model was performed by comparing the simulation to observed data for the 1986-1995 period.Overall, the model follows very similar trends, with the peak flows of some years being slightly overpredicted, and the low flows of the winter months falling to zero for most simulation years.Apart from not being the primary focus of calibration, the discrepancies between the simulated and observed curves for the validated timespan of the SWAT model may be due to the short period of calibration/validation.In general, model calibration in both the tributary and watershed outlets ensured proper apportioning of precipitation and soil water into surface runoff, actual evapotranspiration, and groundwater recharge.This improved model performance compared to the pre-calibration model (Table 2, Figure 5a).Overall, 48% of the observed river flow data were captured by the simulated 95PPU and the average r-factor was about 0.75 at the watershed scale.The final summary statistics for the SWAT calibration and validation process at each stream gauge are given below (Table 2).The lower performance at the tributary gauge may be due to inadequate quality and quantity of input data, which is not representative at the study's level of detail.The spatial and temporal quality and quantity of geospatial maps and climate time series are often not able to perform ideally for each individual sub-basin (including upstream tributaries with small drainage areas).This area of uncertainty has been addressed by earlier studies [3,44].Due to the focus on applying new scenarios to a coupled SWAT-MODFLOW model, the summary statistics below were taken to be satisfactory for the purpose of coupling with the MODFLOW model.

MODFLOW
As outlined in Section 2.3.2, the MODFLOW hydraulic head output was calibrated using 379 points, consisting of two groundwater observation wells, water wells and the two stream gauges used to calibrate the SWAT model.Throughout the calibration process, it was found that the most sensitive parameter of the model was hydraulic conductivity.While maintaining values that were realistic representative [22], the hydraulic conductivity values were modified to create heterogeneous layers representative of the geological setting.
Calibration progress was measured by comparing the simulated hydraulic heads to observed hydraulic heads (Figure 6), with several commonly used summary statistics for each time step including the normalized RMS, the correlation coefficient, the standard error and the residual mean.For most of these statistics, lower values indicate a more accurate calibration (with the exception of the correlation coefficient, where values closer to 1 indicate a better fit).As more detailed recharge data were added from the SWAT model to the MODFLOW model, the results of the calibration for each time step improved, demonstrating the positive reaction of the model to more variable recharge.The best performance statistics we obtained for the entire region, after many trials of manual parameter change, were 0.96, 6.5, and −21.3 for R 2 , NRMS, and Rm, respectively.This proved that, even though these specific recharge values would be replaced by those of the SWAT model, the MODFLOW model performs well under the same conditions as are used in the SWAT model, and the Normalized root mean square (RMS) error remains relatively constant between 6% and 7%.Although areas).This area of uncertainty has been addressed by earlier studies [3,44].Due to the focus on applying new scenarios to a coupled SWAT-MODFLOW model, the summary statistics below were taken to be satisfactory for the purpose of coupling with the MODFLOW model.As outlined in Section 2.3.2, the MODFLOW hydraulic head output was calibrated using 379 points, consisting of two groundwater observation wells, water wells and the two stream gauges used to calibrate the SWAT model.Throughout the calibration process, it was found that the most sensitive parameter of the model was hydraulic conductivity.While maintaining values that were realistic representative [22], the hydraulic conductivity values were modified to create heterogeneous layers representative of the geological setting.
Calibration progress was measured by comparing the simulated hydraulic heads to observed hydraulic heads (Figure 6), with several commonly used summary statistics for each time step including the normalized RMS, the correlation coefficient, the standard error and the residual mean.For most of these statistics, lower values indicate a more accurate calibration (with the exception of the correlation coefficient, where values closer to 1 indicate a better fit).As more detailed recharge data were added from the SWAT model to the MODFLOW model, the results of the calibration for each time step improved, demonstrating the positive reaction of the model to more variable recharge.The best performance statistics we obtained for the entire region, after many trials of manual parameter change, were 0.96, 6.5, and −21.3 for R 2 , NRMS, and Rm, respectively.This proved that, even though these specific recharge values would be replaced by those of the SWAT model, the MODFLOW model performs well under the same conditions as are used in the SWAT model, and the Normalized root mean square (RMS) error remains relatively constant between 6% and 7%.

SWAT-MODFLOW
Using the calibrated component models, a coupled SWAT-MODFLOW model was created using a similar method as in Bailey et al. [4].The median SWAT-simulated 95PPU results were used, out of 100 runs from the best parameter ranges, to couple to the MODFLOW model.To ensure accuracy, the river flow results for the SWAT-MODFLOW model were compared to observed data at two hydrometric stations.We found that the calibrated SWAT model consistently under-predicted river discharge in low-flow periods (winter months), likely due to the limited ability of the model to take GW discharge into account, with flux values dropping to 0 m 3 /s in most cases (Figure 5).However, river flux values of above 0 m 3 /s were simulated after use of the SWAT-MODFLOW model, indicating influence from the GW system and necessity of incorporating the MODFLOW component.The simulated streamflow of the coupled model had similar accuracy to that of the SWAT model by itself at both stream gauges (Figure 7), with the R 2 statistic at the outlet and tributary gauges for this run of SWAT-MODFLOW were 0.54 and 0.484 respectively.Figure 6.Calibration graph for MODFLOW performance, comparing simulated hydraulic head values to the observed hydraulic head values at each observation well.Similar graphs can be displayed at any time step throughout the simulation period.

SWAT-MODFLOW
Using the calibrated component models, a coupled SWAT-MODFLOW model was created using a similar method as in Bailey et al. [4].The median SWAT-simulated 95PPU results were used, out of 100 runs from the best parameter ranges, to couple to the MODFLOW model.To ensure accuracy, the river flow results for the SWAT-MODFLOW model were compared to observed data at two hydrometric stations.We found that the calibrated SWAT model consistently under-predicted river discharge in low-flow periods (winter months), likely due to the limited ability of the model to take GW discharge into account, with flux values dropping to 0 m 3 /s in most cases (Figure 5).However, river flux values of above 0 m 3 /s were simulated after use of the SWAT-MODFLOW model, indicating influence from the GW system and necessity of incorporating the MODFLOW component.The simulated streamflow of the coupled model had similar accuracy to that of the SWAT model by itself at both stream gauges (Figure 7), with the R 2 statistic at the outlet and tributary gauges for this run of SWAT-MODFLOW were 0.54 and 0.484 respectively.

Figure 7.
Comparison of the monthly simulated vs observed streamflow for the watershed outlet (Outlet Gauge in Figure 3).
The outputs produced by the coupled model fall within the acceptable range of values for the study area, as the simulated river flow values produced by the coupled SWAT-MODFLOW model faithfully reproduce the results of the calibrated SWAT model (Figure 6).With regards to the GW-SW exchange flow rates, even MODFLOW grid cells with discharge values of over 50,000 m 3 /d convert to m 3 /s values of less than 1, which are minor in comparison to the river discharges in the area, particularly at the outlet gauge.
A key output of the SWAT-MODFLOW model is detailed GW-SW interaction for each MODFLOW river cell in the watershed (see Table A5).This output is yielded at each time step (daily), with flux values recorded with respect to the GW system.As such, positive values indicate recharge into the GW system, and negative values indicate discharge from groundwater to rivers [34].For nearly every river cell in the study area, the magnitude of simulated discharge varied considerably.These variations were based on a multitude of factors, but primarily owe to the hydraulic head found at each of these river cells.The hydraulic conductivity of the aquifers can also affect the simulated discharge, especially in areas of variable conductivity (Figure 8).
The simulated GW-SW exchange flux values varied spatially, but were primarily negative, indicating discharge from the aquifers.Positive values reached into the hundreds of cubic meters per day (to a maximum of 794 m 3 /d per grid cell), whereas negative values, in many grid cells and time steps, reached into the tens of thousands of cubic meters per day (to −62,479 m 3 /d per grid cell).This indicates a significant interaction occurring between the GW and SW systems of this watershed and The outputs produced by the coupled model fall within the acceptable range of values for the study area, as the simulated river flow values produced by the coupled SWAT-MODFLOW model faithfully reproduce the results of the calibrated SWAT model (Figure 6).With regards to the GW-SW exchange flow rates, even MODFLOW grid cells with discharge values of over 50,000 m 3 /d convert to m 3 /s values of less than 1, which are minor in comparison to the river discharges in the area, particularly at the outlet gauge.
A key output of the SWAT-MODFLOW model is detailed GW-SW interaction for each MODFLOW river cell in the watershed (see Table A5).This output is yielded at each time step (daily), with flux values recorded with respect to the GW system.As such, positive values indicate recharge into the GW system, and negative values indicate discharge from groundwater to rivers [34].For nearly every river cell in the study area, the magnitude of simulated discharge varied considerably.These variations were based on a multitude of factors, but primarily owe to the hydraulic head found at each of these river cells.The hydraulic conductivity of the aquifers can also affect the simulated discharge, especially in areas of variable conductivity (Figure 8).
The simulated GW-SW exchange flux values varied spatially, but were primarily negative, indicating discharge from the aquifers.Positive values reached into the hundreds of cubic meters per day (to a maximum of 794 m 3 /d per grid cell), whereas negative values, in many grid cells and time steps, reached into the tens of thousands of cubic meters per day (to −62,479 m 3 /d per grid cell).This indicates a significant interaction occurring between the GW and SW systems of this watershed and will prove useful for providing dynamic data to future GW-SW projections using the workflows applied here.
The simulated GW-SW exchange data (Figure 8) show that overall exchange trends remain fairly constant throughout the simulation period, while the southern (upstream) portion of the river showed the most pronounced fluctuation between the wet and dry months, with lower relative recharge from the river occurring during January, and more river seepage (i.e., less aquifer discharge) taking place in June.In addition, the tributary found within Sub-basin 13 of the study area consistently had the highest discharge rates observed.The reason for these simulated patterns may be due to the regional bedrock strata: more formations pinch out and outcrop to the north (see Figure 3), so the changes in flow patterns may be partially caused by the hydraulic conductivity of layers that are close to the surface.
Water 2019, 11, x FOR PEER REVIEW 14 of 30 will prove useful for providing dynamic data to future GW-SW projections using the workflows applied here.The simulated GW-SW exchange data (Figure 8) show that overall exchange trends remain fairly constant throughout the simulation period, while the southern (upstream) portion of the river showed the most pronounced fluctuation between the wet and dry months, with lower relative recharge from the river occurring during January, and more river seepage (i.e., less aquifer discharge) taking place in June.In addition, the tributary found within Sub-basin 13 of the study area consistently had the highest discharge rates observed.The reason for these simulated patterns may be due to the regional bedrock strata: more formations pinch out and outcrop to the north (see Figure 3), so the changes in flow patterns may be partially caused by the hydraulic conductivity of layers that are close to the surface.

Model Application for Examining GW Pumping and Climate Change Scenarios
As is discussed in Section 2.2, 21 wells were included in the SWAT-MODFLOW model (Figure 2a), with pumping rates of 468 m 3 /d.Screens for each well were assigned to one of the two aquifers in the study area: either the Paskapoo or the Wapiti formation, each already containing water wells.The pumping rate was taken from an example water extraction well drilled by an operator in the

Model Application for Examining GW Pumping and Climate Change Scenarios
As is discussed in Section 2.2, 21 wells were included in the SWAT-MODFLOW model (Figure 2a), with pumping rates of 468 m 3 /d.Screens for each well were assigned to one of the two aquifers in the study area: either the Paskapoo or the Wapiti formation, each already containing water wells.The pumping rate was taken from an example water extraction well drilled by an operator in the study area, and was active in the model for half of each year after the 3-year warm-up period, beginning at the start of each simulation year.
As expected, hydraulic heads in the regions with wells decrease, with a difference of up to 5 m of drawdown observed in comparison with the SWAT-MODFLOW model without pumping wells.In addition, the influence of water well pumping had an effect on the observed cell-by-cell GW-SW flow rates.On average, the aquifer discharge rate to the river decreased, with the simulated flow rate changing from −1294 m 3 /d to −1261 m 3 /d.This discharge rate decreased in 329 of the 405 river cells, and the maximum observed decrease in discharge was 938 m 3 /d.
The GW-SW interaction change observed is likely due to the fact that as water is extracted from the subsurface, less of the remaining water is available to provide baseflow to rivers.Depending on the rate and schedule of pumping (including the number of pumping wells), these results show that the effects of pumping are immediately observable, and may affect the discharge rates of associated rivers.However, at the tested pumping rate of 468 m 3 /d, the discharge decreases in terms of the overall discharge rates (~150 m 3 /s) of the Little Smoky River are relatively insignificant.This suggests that the bedrock aquifers within the study area could support water well pumping at rates similar to those tested.
A climate change run over the 2010-2034 period was the next to be applied.The coupled SWAT-MODFLOW model was run five separate times, each with a different downscaled GCM dataset (see Section 2.4).The results of each run were then averaged to yield one harmonized output, with the GW-SW exchange rates shown in Figure 9. Immediately noticeable in the figure is that the long-term average annual historical results for GW-SW interaction patterns look similar to the results of the climate change scenarios.This is indeed the case, as the overall average historical flow rate was −1294.04m 3 /d, while that of the averaged climate change scenarios was −1294.14m 3 /d (both indicating overall discharge into the river).These results indicated that the effects of climate change on this watershed over the 2010-2034 period are negligible.A possible reason for the lack of influence of the climate change scenario is that the simulation period was quite short in terms of overall climate trends, with little change occurring in the predicted CO 2 concentration over this timespan for the RCP8.5 scenario, as well as all other scenarios (~375 to ~475 ppm; 41).Although the simulation used in this study for scenario analysis was an average, it should be noted that each individual coupled climate change model yielded relatively similar simulation results, again likely owing to the short-term future projection.
As the most complex scenario to be tested in this study, the five coupled climate change scenarios were run again with a MODFLOW model that included the higher pumping rate of 4680 m 3 /d per well.The effects of the elevated pumping rates caused a more noticeable change in the GW-SW interaction pattern of the watershed for the 2010-2034 forecasts.The substantial increase in the volume of pumped water left less stored GW in the primary aquifers (i.e., the Paskapoo Fm.), as the overall discharge rates for each river cell are observed to decrease (Figure 9).This storage trend can also be observed in the MODFLOW model before coupling, with model-wide fluxes from storage spiking to correspond with pumping intervals (from un-pumped totals of under 1 m 3 /d up to total maximums of ~1,000,000 m 3 /d).While the total average discharge into the rivers was 1294 m 3 /d in the initial SWAT-MODFLOW model , this average was reduced to 1174 m 3 /d after simulating with the RCP8.5 GCMs and the elevated pumping schedule, with cell-by-cell flowrate differences ranging from under 5 m 3 /d to over 3000 m 3 /d (discharge decreasing in 320 of 405 river cells).A comparison of the change in aquifer discharge under the tested simulation scenarios, i.e., with pumping and without pumping, indicates that anthropogenic influence can have more immediate effects on the watershed than those of climate change, even at the lower rate of 468 m 3 /d.
Although water well pumping simulations were the only anthropogenic influence applied to this study, the change in the GW-SW exchange was measurable on a regional scale.However, the environmental footprint in the study area will only grow when considering additional sectors such as agriculture.The fact that the potential effects of these operations can now be quantified with respect to both GW and SW resources is a key benefit of using a coupled model.
The ability of SWAT-MODFLOW to simulate the hydro(geo)logy of a watershed is fully dependent on the accuracy of its component models.Due to much time being spent on the calibration of both the SW and GW models, the coupled result proves to be robust in terms of its representation of river discharges (Figure 7).Similar to the Sprague River model built by Bailey et al. [4], high spatial variation in GW-SW exchange rates can be seen in each SWAT-MODFLOW figure in Section 3. The ability to gain information on these patterns may prove to be crucial for the informed management of both SW and GW resources in this watershed, as well as any other in which a GW-SW model has been applied.Moreover, when additional stressors such as climate change and pumping simulations are included, it yields more information from a risk management perspective, which can aid operators and regulators to better understand the effects of anthropogenic influence within a watershed.The ability of SWAT-MODFLOW to simulate the hydro(geo)logy of a watershed is fully dependent on the accuracy of its component models.Due to much time being spent on the calibration of both the SW and GW models, the coupled result proves to be robust in terms of its representation of river discharges (Figure 7).Similar to the Sprague River model built by Bailey et al. [4], high spatial variation in GW-SW exchange rates can be seen in each SWAT-MODFLOW figure in Section 3. The ability to gain information on these patterns may prove to be crucial for the informed management of both SW and GW resources in this watershed, as well as any other in which a GW-SW model has been applied.Moreover, when additional stressors such as climate change and pumping simulations are included, it yields more information from a risk management perspective, which can aid operators and regulators to better understand the effects of anthropogenic influence within a watershed.

Limitations
If each component model is built and calibrated effectively, SWAT-MODFLOW acts as a powerful tool for simulating integrated GW-SW interaction.Although this has been established in the current study and those previously done [4,6,11,12], the current novelty of the model and

Limitations
If each component model is built and calibrated effectively, SWAT-MODFLOW acts as a powerful tool for simulating integrated GW-SW interaction.Although this has been established in the current study and those previously done [4,6,11,12], the current novelty of the model and dependence on its components give rise to some limitations.The performance of a SWAT-MODFLOW model is governed by the quality and refinement of its components.Therefore, in study areas that may be lacking in either GW-or SW-related inputs or observational data, the ability to produce a faithful model may be limited.In addition, this coupled model encounters a limitation shared with all models: the trade-off between scientific accuracy and computational burden while keeping the study goal in mind.One has the flexibility with SWAT-MODFLOW to discretize both of its component models to any desired refinement [4,12].However, the goals of the study (and lack of access to computational power/data) may limit the ability to produce models of such fine detail.In such cases, the limitations of SWAT and MODFLOW may, in fact, resurface.The coupled model used the best simulation results for SWAT, and the final results of manual calibration for MODFLOW, which were deemed to be satisfactory for the purposes of this study.As these components are subject to non-uniqueness, future studies using multiple calibration methods are an opportunity to address uncertainty in the coupled model.In addition, the SW component of SWAT-MODFLOW, SWAT, is primarily used for regional studies of hydrology due to its efficient HRU-based computational scheme [3,4].Because of this, the ability of SWAT-MODFLOW to simulate short-term site-specific processes may be limited.It is therefore important for an adequate model selection process to occur so that the goals of the study are met.
Study-specific limitations of our SWAT-MODFLOW model include the fact that the component models were calibrated over a relatively brief period with respect to the overall hydro-climate.As the calibrated parameters reflect short-term dynamics over 25 years, large-scale temporal variations may not have been taken into consideration.In light of this, further studies are required to assess the robustness of this model when subjected to a wider range of natural climate variability, including the addition of more climate change scenarios.To more completely assess the various demands within the water-food-energy nexus, water conflicts between sectors will also be a growing area of study.As such, conducting research where the water extraction rates from the petroleum, agriculture and domestic sectors (in addition to the needs of ecosystems) are all considered under various climate change scenarios will increase the ability of decision-makers to quantify the interrelated demands on water as a holistic GW-SW system

Conclusions
It is becoming increasingly important to understand the interrelated nature of water resources, and to consider them in a unified manner.This is especially true in areas like Alberta, where significant demands are placed on these resources from multiple sectors.Hydro(geo)logic models can help us to forecast and visualize complex scenarios, as well as to make more informed decisions about the management of these resources.Coupled or integrated hydro(geo)logic models are one of the critical emerging tools for this purpose, as they track the interactions between both the vital components in the hydrologic system.
This study tested the robustness of one such model, SWAT-MODFLOW, in an area of diverse climate and hydrogeology, and included both pumping and climate change to determine the relative effects on the study watershed.Although some fully integrated models may offer a more faithful representation of the site conceptual model [8], the relatively user-friendly and cost-effective code tested here performed successfully under highly variable hydro(geo)logic and weather conditions while incorporating a deeper, more complex bedrock system than any prior SWAT-MODFLOW model.The model also included additional factors, such as snow influence, that had not been included in previous SWAT-MODFLOW studies.The implementation of this SWAT-MODFLOW model proved that GW has a heavy influence on the hydrology of the Little Smoky River watershed, discharging significant volumes of water into the rivers and tributaries of the study area.This reinforces the importance of examining SW and GW as one system, as the demands placed on one have direct effects on the other.
The modelling workflows developed here can be used to apply new GW-SW interaction scenarios with a multitude of both natural and anthropogenic influences.Use of this model can help scientists track water resources more accurately, and can inform the decisions made by policymakers whose priority is to ensure that this vital resource can remain available to all who need it.1. Export the MODFLOW grid, and import it into SWAT as a shapefile.
2. Geo-reference the grid (was done manually in this study) to overlap with the watershed of the study area.3. Right click the MODFLOW grid shapefile, select "Joins and Relates," and then "Join".4. Select "Join data from another layer based on spatial location."5. Choose the point feature class as the joining layer.6. Choose the second option: Each polygon will be given all the attributes of the line that is closest to its boundary, and the result will be a joint layer.The distance filed (last one in the attribute table ) shows how close the line is to each cell.Therefore, a line falling inside a polygon is treated as being closest to the polygon (i.e., a distance of 0). 7. Select the cells you are interested in by filtering based on the "Distance" filed.8. Export the selected cells.3. Delete all columns except "HRU," "SUB," "MON," and "GW_RCHGmm."This will yield a spreadsheet with monthly GW recharge values for the entire simulation period.4. Delete the yearly average rows in each column.5. Create two columns for Time 1 (T1) and Time 2 (T2), with the initial MODFLOW day in T1 and the final MODFLOW day in T2 for each monthly stress period.Add rows until the ending day of simulation.Repeat for each sub-basin.6. Organize columns in the following format: Geo-reference the grid (was done manually in this study) to overlap with the watershed of the study area.
Select "Join data from another layer based on spatial location."5.
Choose the point feature class as the joining layer.6.
Choose the second option: Each polygon will be given all the attributes of the line that is closest to its boundary, and the result will be a joint layer.The distance filed (last one in the attribute table ) shows how close the line is to each cell.Therefore, a line falling inside a polygon is treated as being closest to the polygon (i.e., a distance of 0). 7.
Select the cells you are interested in by filtering based on the "Distance" filed.8.
Export the selected cells.
Workflow 2. Complete process of creating a .RCH file from the output of a SWAT model to import into the MODFLOW model of the project.

1.
Open the "output.hru"file for the calibrated SWAT model in Microsoft Excel.

2.
Perform the Text to Columns operation if necessary (to separate data by column).
Delete the yearly average rows in each column.

5.
Create two columns for Time 1 (T1) and Time 2 (T2), with the initial MODFLOW day in T1 and the final MODFLOW day in T2 for each monthly stress period.Add rows until the ending day of simulation.Repeat for each sub-basin.6.
Organize columns in the following format:  Table 5. List of SWAT-MODFLOW outputs and their contents (modified from Beets [46]).Table A5.List of SWAT-MODFLOW outputs and their contents (modified from Beets [46]).In the model built in this study, the same protocol was followed as in the SWAT-MODFLOW manual [34]: each identifier number was changed to be 5000 more than the original identifier number (16 would become 5016, etc.).After completely calibrating both models separately (explained further in Section 2.3), the SWAT-MODFLOW model was run.The model took roughly six hours to run completely-significantly less time than the initial hypothesis predicted.Each subsequent run after the initial was performed on a more powerful computer, which improved the run time of the coupled model.These runs took slightly less than four hours for a complete simulation, allowing more than one SWAT-MODFLOW run to occur in a day.

Appendix B.3 Hydraulic Head Calibration
To evaluate the overall representation of the observed data, various statistics were used, similar to the calibration process of the SWAT model.The principal simulation statistic used was the normalized root mean square error (RMS), which evaluates the overall closeness of a simulated dataset to the observed dataset [47].Generally given as a percentage, values closer to 0% indicate less deviation from the observed dataset.The model was also calibrated using the average correlation coefficient, a value that evaluates the linear relationship between two variables [48], which in the case of this study are observed head and simulated head.Ranging from −1 to 1, a value closer to 1 represents a better correlation.
The vast number of observation points and the Visual MODFLOW Graphical User Interface (GUI) allowed the pinpointing of locations with larger hydraulic head discrepancies, which facilitated parameter adjustment for future iterations.Many parameters and boundary conditions can be changed to achieve satisfactory calibration/validation results, including specific storage, specific yield and hydraulic head-related boundaries, but the key parameter in many MODFLOW models is the hydraulic conductivity of the layers.As these values need to be changed largely based on region, it is therefore important that the values of hydraulic conductivity in a given model be based on those found in the corresponding subsurface units to the extent possible.To that end, the calibrated values used in this study were kept within the ranges found in studies such as those by Hughes et al. [22] and Smerdon et al. [29].

Figure 1 .
Figure 1.The study area (Little Smoky River Watershed), shown within the province of Alberta, Canada.

Figure 1 .
Figure 1.The study area (Little Smoky River Watershed), shown within the province of Alberta, Canada.

Figure 2 .
Figure 2. Little Smoky River watershed; (a) displaying the delineated sub-basins, hydrometric stations, and simulated pumping well locations (based on water well trends in study area) included.Each well is proximal to or within the study area, with 16 of the 21 pumping wells located in the downstream portion of the modelled area; (b) displays the topographic classes and varying elevation.

Figure 2 .
Figure 2. Little Smoky River watershed; (a) displaying the delineated sub-basins, hydrometric stations, and simulated pumping well locations (based on water well trends in study area) included.Each well is proximal to or within the study area, with 16 of the 21 pumping wells located in the downstream portion of the modelled area; (b) displays the topographic classes and varying elevation.

Figure 3 .
Figure 3. South-north cross-section of the MODFLOW modelled area, displaying the geological units present.

Figure 3 .
Figure 3. South-north cross-section of the MODFLOW modelled area, displaying the geological units present.

Figure 4 .
Figure 4. Contour map of the hydraulic head for the Little Smoky River watershed, produced by Visual MODFLOW.The hydraulic head roughly follows topography, with higher heads found in the upstream regions of the area.Each small number on the map represents a calibration well within the Little Smoky River watershed.

Figure 4 .
Figure 4. Contour map of the hydraulic head for the Little Smoky River watershed, produced by Visual MODFLOW.The hydraulic head roughly follows topography, with higher heads found in the upstream regions of the area.Each small number on the map represents a calibration well within the Little Smoky River watershed.

Figure 5 .
Figure 5.Comparison of the monthly simulated (green band) vs observed (blue signal) streamflow for the Outlet Gauge in the Little Smoky River Watershed over the entire historical simulation period.(a) displays a simulation with the default basin parameters, and (b) shows the calibration and validation results.

Figure 5 .
Figure 5.Comparison of the monthly simulated (green band) vs. observed (blue signal) streamflow for the Outlet Gauge in the Little Smoky River Watershed over the entire historical simulation period.(a) displays a simulation with the default basin parameters, and (b) shows the calibration and validation results.
the groundwater recharge output values from the SWAT model are used directly in the MODFLOW model, each model will yield slightly different results based on their governing flow mechanisms (HRU-based flow vs. grid-based Finite-Difference flow).Water 2019, 11, x FOR PEER REVIEW 12 of 30 Although the groundwater recharge output values from the SWAT model are used directly in the MODFLOW model, each model will yield slightly different results based on their governing flow mechanisms (HRU-based flow vs grid-based Finite-Difference flow).

Figure 6 .
Figure 6.Calibration graph for MODFLOW performance, comparing simulated hydraulic head values to the observed hydraulic head values at each observation well.Similar graphs can be displayed at any time step throughout the simulation period.

Figure 7 .
Figure 7.Comparison of the monthly simulated vs. observed streamflow for the watershed outlet (Outlet Gauge in Figure 3).

Figure A1 .Workflow 1 .
Figure A1.(A) Little Smoky River catchment soil types, and (B) land use types.The soil and land use types roughly correspond to the elevations at which they are found.

Workflow 2 . 1 .
Complete process of creating a .RCH file from the output of a SWAT model to import into the MODFLOW model of the project.Open the "output.hru"file for the calibrated SWAT model in Microsoft Excel.2. Perform the Text to Columns operation if necessary (to separate data by column).

FigureWorkflow 1 .
Figure A1.(A) Little Smoky River catchment soil types, and (B) land use types.The soil and land use types roughly correspond to the elevations at which they are found.Workflow 1. Complete process of exporting river data from a SWAT model to create a .RIV boundary for a MODFLOW model.1.Export the MODFLOW grid, and import it into SWAT as a shapefile.2.Geo-reference the grid (was done manually in this study) to overlap with the watershed of the study area.3.Right click the MODFLOW grid shapefile, select "Joins and Relates," and then "Join".4.Select "Join data from another layer based on spatial location."5.Choose the point feature class as the joining layer.6.Choose the second option: Each polygon will be given all the attributes of the line that is closest to its boundary, and the result will be a joint layer.The distance filed (last one in the attribute table)shows how close the line is to each cell.Therefore, a line falling inside a polygon is treated as being closest to the polygon (i.e., a distance of 0).7.Select the cells you are interested in by filtering based on the "Distance" filed.8.Export the selected cells.

Figure A3. 3 -
Figure A3.3-D diagram of the MODFLOW model, with each color representing the hydraulic conductivity found at that point.Many of the model layers have more than one conductivity value.
Content swatmf_out_MF_gwsw GW-SW exchange (m 3 /day) for each MODFLOW River cell.Positive values are for river water entering the GW system, and negative values are for GW entering the river.swatmf_out_MF_recharge Recharge (m 3 /day) to the water table for each grid cell, for the day.The cell values are printed as a grid, at each daily time step of the simulation.swatmf_out_MF_riverstage MODFLOW river stage (m) for each river cell, listed in order of the river cell ID.swatmf_out_SWAT_channel Channel depth (m) of each sub-basin stream or river, listed in order of sub-basin # across the columns (i.e., each row in the file is one simulation day).swatmf_out_SWAT_gwsw Exchange (m 3 /day) for each sub-basin, listed in order of sub-basin # (1 to ...) for the day.Positive values are for GW entering the river, and negative values are for river water entering the aquifer.swatmf_out_SWAT_recharge Deep percolation (mm) calculated for each HRU.The values are listed in order of HRU # (1 to ...) for the day and are printed out for each day of the simulation.

Figure A3. 3 -
Figure A3.3-D diagram of the MODFLOW model, with each color representing the hydraulic conductivity found at that point.Many of the model layers have more than one conductivity value.
Content swatmf_out_MF_gwsw GW-SW exchange (m 3 /day) for each MODFLOW River cell.Positive values are for river water entering the GW system, and negative values are for GW entering the river.swatmf_out_MF_recharge Recharge (m 3 /day) to the water table for each grid cell, for the day.The cell values are printed as a grid, at each daily time step of the simulation.swatmf_out_MF_riverstage MODFLOW river stage (m) for each river cell, listed in order of the river cell ID.swatmf_out_SWAT_channel Channel depth (m) of each sub-basin stream or river, listed in order of sub-basin # across the columns (i.e., each row in the file is one simulation day).swatmf_out_SWAT_gwsw Exchange (m 3 /day) for each sub-basin, listed in order of sub-basin # (1 to ...) for the day.Positive values are for GW entering the river, and negative values are for river water entering the aquifer.swatmf_out_SWAT_recharge Deep percolation (mm) calculated for each HRU.The values are listed in order of HRU # (1 to ...) for the day and are printed out for each day of the simulation.

Figure A4 .
Figure A4.Comparison of simulated SWAT-MODFLOW results to observed values at the Little Smoky River watershed tributary gauge.
* Canadian Earth System Model Version 2, Community Climate System Model 4.0, Centre National de Recherches Météorologiques Climate Model Version 5, Commonwealth Scientific and Industrial Research Organization MK3, Model for Interdisciplinary Research on Climate Version 5.

Table 2 .
Final combined summary statistics for the calibrated and validated SWAT model.Initial refers to prior to calibration.

Table 2 .
Final combined summary statistics for the calibrated and validated SWAT model.Initial refers to prior to calibration.

Table A2 .
Input data required for the MODFLOW model, and their sources.

Table A6 .
Summary statistics for the MODFLOW component.