Next Article in Journal
Five-Year Experimental Study on Effectiveness and Sustainability of a Dry Drainage System for Controlling Soil Salinity
Next Article in Special Issue
Water Footprint and Water Pinch Analysis in Ethanol Industrial Production for Water Management
Previous Article in Journal
Effect of Mimic Vegetation with Different Stiffness on Regular Wave Propagation and Turbulence
Previous Article in Special Issue
Assessing the Impacts of Population Growth and Climate Change on Performance of Water Use Systems and Water Allocation in Kano River Basin, Nigeria
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, AB, T6G 2E3, Canada
2
Alberta Geological Survey, Alberta Energy Regulator, Edmonton, AB, T6B 2X3, Canada
*
Author to whom correspondence should be addressed.
Water 2019, 11(1), 110; https://doi.org/10.3390/w11010110
Submission received: 11 December 2018 / Revised: 24 December 2018 / Accepted: 3 January 2019 / Published: 10 January 2019

Abstract

:
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.

1. 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.
Collaborative measures must be taken to manage limited water resources to meet the increasing water demands of every sector. Eco-hydro-geologic modelling can be useful in quantifying the resources and needs of specific areas and understanding the natural and anthropogenic factors affecting their distribution in time and space [3,4,5]. Depending on the scope of the research studies, models can be applied to regions of significantly varying size and complexity. In order to make accurate forecasts, computer models that quantify water resources should be holistic, focusing on groundwater (GW) and surface water (SW) together instead of treating them individually. In the last two decades, the number of models that considered GW and SW together has increased [6,7,8,9]. Developing models that are consistently accurate is a key step toward more intelligent water management. To achieve this objective, considerations must be made during the selection process to ensure that the model is able to meet the demands of the study, including data availability, ease of use, desired accuracy of simulations and budget. Of those available, four coupled GW–SW models have been used in multiple published basin-scale studies. These include HydroGeoSphere [8], MIKE SHE [10], GSFLOW [7] and SWAT–MODFLOW [4,6,11,12]. While the first two of these models are capable of incredibly complex simulations, they are both fully distributed and grid-based, therefore requiring significant amounts of high-resolution data and computing capacity. The latter two are semi-distributed, and incorporate the most important physical processes to efficiently simulate large-scale basins over longer timespans.
SWAT–MODFLOW is a coupled hydro(geo)logical model that employs both the Soil and Water Assessment Tool—SWAT [13] and MODFLOW [14] (SW and GW models, respectively) to yield an integrated output. Coupling the component models covers the limitations inherent in each component model to yield a solution that is more faithful to real-world hydro(geo)logy. SWAT covers processes associated with SW hydrology, such as precipitation, temperature, river flow, surface runoff, soil water, actual evapotranspiration and GW recharge, while MODFLOW is responsible for the GW processes, including saturated flow and GW discharge into rivers.
Although multiple SWAT–MODFLOW codes have been developed over the past twenty years for study areas of highly varying scale [4,6,12,15,16], each study was carried out to overcome the shortcomings of each of its component models (SWAT and MODFLOW), and to evaluate SW and GW resources as one system. While SWAT simulates physical processes associated with soil, plants, water in the root zone and weather, MODFLOW is responsible for providing outputs for three-dimensional GW flow, as well as groundwater pumping simulations. Each SWAT–MODFLOW code follows the same general logic, in that MODFLOW is run as a subroutine of SWAT [4,11,12], effectively replacing the GW module in SWAT to yield outputs at finer resolution than SWAT can achieve by itself for a given time step. Due to its ability to efficiently simulate multiple physical processes in a large region with no associated cost, this study uses the publicly available SWAT–MODFLOW tool developed by Bailey et al. [4].
Models developed in previous studies covered smaller study areas with less complex hydrogeology, and were mainly intended for understanding local water changes [4,6,11,12]. Prior to this study, the Sprague River Watershed model [4] was the largest at roughly 4100 km2. Other models range in size to fit various goals, with some as small as 198 km2 [11], and most having study areas in the high hundreds to low thousands of square kilometers [12]. In addition, the successfully implemented SWAT–MODFLOW models thus far have had hydrology and hydrogeology of low complexity, with each of the MODFLOW components focusing on shallow, local GW flow patterns [11,12].
This study uses a portion of western Alberta, Canada as a case study to simulate a hydro(geo)logically variable watershed (Little Smoky River). Since the rise in the economic potential of hydraulic stimulation in combination with horizontal drilling to develop unconventional resources, the region has become a key area for Alberta’s oil and gas development [17]. For many operators in the area, the Duvernay Formation is targeted for its unconventional resources, with wells using as much as 30,000 m3 of water to effectively fracture the reservoir [18], although this value can increase further when targeting formations elsewhere in western Canada [19]. Historically, much of the fresh water required for operations was sourced from the rivers and streams in the region [20], with occasional extraction of GW from freshwater (or near-fresh) aquifers. With increased unconventional resource development comes the potential to use GW to support hydraulic fracturing operations [21,22]. However, a lack of understanding of the total GW and SW availability, reliability and interactions, and limited available projections under future climate change scenarios, may restrict future increases in production of energy resources in the region. Therefore, the integration of both SW and GW supply-demand processes in the region is needed to systematically study the future opportunities and conflicts for industry. Due to similarly increasing water use trends in other semi-arid basins with strongly defined seasons, such as North Dakota [23], the findings and workflows of this study can be applied to evaluate GW–SW interactions in many similar watersheds worldwide.
By virtue of being in western Alberta, the case study area is subjected to highly variable, yet predictable, climate patterns. Like most Canadian watersheds, the Little Smoky River experiences higher discharge during the spring and summer months, with next to no flow during the winter months. This, in turn, affects the hydrology of the basins in the study area. With the headwaters of the Little Smoky River located at elevations of more than 1200 m above sea level, snowfall is an influence that must be considered in this area for the most realistic representation of the climate. In addition, the underlying Western Canadian Sedimentary Basin (WCSB) is comprised of geological formations that dip upwards from west to east, creating a wedge-shape geological setting in the study area. Many of the formations are both thicker and deeper along the Rocky Mountains but become shallower and are exposed or only covered by a thin layer of sediments to the northeast. Combined with these dipping layers, variability of rock properties within a given geological formation also varies spatially, resulting in a complex hydrogeological environment with widely varying hydraulic conductivity and storage properties [22,24]. These complex but realistic hydrologic and hydrogeologic conditions and the importance of water to support unconventional resource development and sustain rivers, makes it an ideal case study to test the SWAT–MODFLOW model.
The primary objective of this study was to apply SWAT–MODFLOW to simulate a large region to study the regional hydrogeological conditions and the potential impacts of climate change and water withdrawal on GW–SW interactions. Our specific objectives were: (1) to calibrate, validate, and provide uncertainty analyses for both model components; (2) to assess the effects of future climate change on the GW–SW exchange pattern; and (3) to examine the impacts of groundwater use scenarios on GW–SW interactions for the interval of 2010–2034. This was done to achieve the goal of further understanding of the use of SWAT–MODFLOW as a planning tool to predict the effect of natural and anthropogenic influences on the watershed.

2. Materials and Methods

2.1. Study Area

2.1.1. Surface Hydrology

The Little Smoky watershed of Alberta (Figure 1) covers 11,494 km2 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.

2.1.2. 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.

2.2. Model Development and Input Data

2.2.1. 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.

2.2.2. 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 km2, 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).
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.

2.2.3. 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.

2.3. Calibration, Validation and Uncertainty Assessment

2.3.1. 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 R2 and bR2. The R2, or coefficient of determination, relates how well the simulated data trend replicates the observed trend. The primary factor observed in this study, however, is bR2, which takes both trend and closeness to observed results into account, and is calculated by multiplying R2 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 bR2 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 bR2, 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.

2.3.2. 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 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.

2.4. 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 m3/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 (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 m3/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 (1983–2007) 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.

3. Results and Discussion

3.1. Sensitivity, Calibration, Validation and Uncertainty Analysis

3.1.1. 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 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 Section 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.

3.1.2. 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 R2, 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 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).

3.2. 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 m3/s in most cases (Figure 5). However, river flux values of above 0 m3/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 R2 statistic at the outlet and tributary gauges for this run of SWAT–MODFLOW were 0.54 and 0.484 respectively.
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 m3/d convert to m3/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 m3/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 m3/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.

3.3. 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 m3/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 m3/d to −1261 m3/d. This discharge rate decreased in 329 of the 405 river cells, and the maximum observed decrease in discharge was 938 m3/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 m3/d, the discharge decreases in terms of the overall discharge rates (~150 m3/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 scenario 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.04 m3/d, while that of the averaged climate change scenarios was −1294.14 m3/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 CO2 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 m3/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 m3/d up to total maximums of ~1,000,000 m3/d). While the total average discharge into the rivers was 1294 m3/d in the initial SWAT–MODFLOW model (1983–2007), this average was reduced to 1174 m3/d after simulating with the RCP8.5 GCMs and the elevated pumping schedule, with cell-by-cell flowrate differences ranging from under 5 m3/d to over 3000 m3/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 m3/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.

3.4. 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

4. 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.

Author Contributions

D.S.A. and M.F. conceived and designed the experiments and supervised the modelling and findings of this work; D.C. performed the experiments and analyzed the data; B.S. assisted in developing the hydrogeological framework; D.C. wrote the initial draft of the paper, and all other co-authors contributed in the revision process of the paper.

Funding

Funding for this work has been received from a Natural Sciences and Engineering Research Council of Canada (NSERC) Collaborative Research and Development grant (CRDPJ 469308-14) in collaboration with the Encana Corporation, Alberta Scholarships, and the Campus Alberta Innovation Program (CAIP) Chair Award (Grant #RES0030781).

Acknowledgments

We would like to thank the Alberta Geological Survey of the Alberta Energy Regulator for collaborating on the regional geological model.

Conflicts of Interest

The authors declare no conflicts of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, nor in the decision to publish the results.

Appendix A. Figures, Tables and Workflows

Table A1. List of major developments in SWAT–MODFLOW use, with alternative models listed below.
Table A1. List of major developments in SWAT–MODFLOW use, with alternative models listed below.
SWAT-MODFLOW Models
YearAuthorTitleKey Developments
1999Sophocleous et al.Integrated numerical modeling for basin-wide water management: The case of the Rattlesnake Creek basin in south-central KansasFirst fully coupled SWAT-MODFLOW model to be applied to a watershed
2008Kim et al.Development and application of the integrated SWAT-MODFLOW modelFirst published SWAT-MODFLOW model to include a detailed water well pumping schedule
2015Guzman et al.A model integration framework for linking SWAT and MODFLOWFirst SWAT-MODFLOW model to have simultaneously built component models (using SWATmf-app)
2016Bailey et al.Assessing regional-scale temporal patterns of groundwater-surface water interactions using a coupled SWAT-MODFLOW modelFirst SWAT-MODFLOW code to be endorsed by SWAT, and publicly available on the SWAT website
Other Available GW-SW Models
YearAuthor and Model NameTitleKey Developments
2008Sudicky et al.: HydroGeoSphereSimulating complex flow and transport dynamics in an integrated surface-subsurface modeling frameworkMost powerful GW-SW tool currently available
2008Markstrom et al.: GSFLOWGSFLOW – Coupled Ground-Water and Surface-Water Flow Model Based on the Integration of the Precipitation-Runoff Modeling System (PRMS) and the Modular Ground-Water Flow Model (MODFLOW-2005).GW-SW model developed by USGS. Direct competitor with SWAT-MODFLOW
Table A2. Input data required for the MODFLOW model, and their sources.
Table A2. Input data required for the MODFLOW model, and their sources.
Data Sources for SWAT Model
TypeDatasetResolution/CoverageTime SpanTime StepRegionReference
ClimateMeteorological StationsN/A1983–2007DailyLocalFaramarzi et al., 2015
Streamflow2 in Study Area1986–2007Daily (monthly calibration)Localhttp://www.ec.gc.ca/rhc-wsc/
Digital MapsDigital Elevation Map30 m × 30 mN/AN/AProvincial(Jarvis et al., 2008)
Land Use Map13 Feature Classes1983–2007N/ALocalhttp://www.geobase.ca/geobase/en/data/landcover/csc2000v/description.html
Soil Map20 Feature Classes1983–2007N/ALocalhttp://sis.agr.gc.ca/cansis/nsdb/slc/index.html
Climate Change DataPrecipitation300 arc sec2010–2034DailyDownscaled to LocalPacific Climate Impacts Consortium (PCIC) Cannon et al., 2015
Temperature300 arc sec2010–2034DailyDownscaled to LocalPacific Climate Impacts Consortium (PCIC) Cannon et al., 2015
Data Sources for MODFLOW Model
TypeDatasetResolution/CoverageTime SpanTime StepRegionReference
ElevationBedrock Formations—South1 km × 1 kmN/AN/ALocalAlberta Geological Survey
Bedrock Formations—NorthPicked Fm. tops. 20 wellsN/AN/ALocalAccuMap
Bedrock PropertiesHydraulic Conductivity, Specific Storage1 km × 1 kmN/AN/ALocalAlberta Geological Survey
Boundary ConditionsConstant HeadN/A1983–2007ConstantProvincial(Brinsky, 2014), SWAT Output
Rivers405 Grid Cells1983–2007ConstantLocalSWAT Output
RechargeSubbasin-scale1986–2007MonthlyLocalSWAT Output
Observation Wells/GaugesStreamflow2 in Study Area1986–2007Daily (monthly calibration)Localhttp://www.ec.gc.ca/rhc-wsc/
Water Well Records377 in Study Area2002–2007SporadicLocalAlberta Environment and Parks
Pumping WellsPumping Rate from Encana Well. Self-placed in model21 in Study AreaDays 1101–9131 of MODFLOW SimulationSemi-yearlyLocalEncana Corporation
Table A3. Weather gauges near the Little Smoky River watershed, with example average humidity and wind speed data for 2018 (modified from Alberta Agriculture and Forestry, 2018).
Table A3. Weather gauges near the Little Smoky River watershed, with example average humidity and wind speed data for 2018 (modified from Alberta Agriculture and Forestry, 2018).
Station NameDateAvg. Humidity (%)Avg. Wind Speed (km/h)
Pass Creek AutoJanuary 201881.21614.408
Pass Creek AutoFebruary 201874.66315.504
Pass Creek AutoMarch 201875.86213.368
Pass Creek AutoApril 201860.83113.703
Pass Creek AutoMay 201855.09813.886
Fox Creek AutoJanuary 201891.7425.351
Fox Creek AutoFebruary 201885.697.517
Fox Creek AutoMarch 201883.4036.156
Fox Creek AutoApril 201868.4817.047
Fox Creek AutoMay 201860.4185.788
Valleyview AGDMJanuary 201880.78813.637
Valleyview AGDMFebruary 201873.68914.826
Valleyview AGDMMarch 201873.70113.008
Valleyview AGDMApril 201864.14514.121
Valleyview AGDMMay 201851.20214.368
Tony AutoJanuary 201879.8969.475
Tony AutoFebruary 201873.37910.158
Tony AutoMarch 201873.3769.599
Tony AutoApril 201858.56310.172
Tony AutoMay 201853.6449.312
Berland Hills AutoJanuary 201877.7355.273
Berland Hills AutoFebruary 201866.3397.704
Berland Hills AutoMarch 201870.7437.092
Berland Hills AutoApril 201859.7748.02
Berland Hills AutoMay 201859.2237.958
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.
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.
Water 11 00110 g0a1
Workflow 1. Complete process of exporting river data from a SWAT model to create a .RIV boundary for a MODFLOW model.
  • Export the MODFLOW grid, and import it into SWAT as a shapefile.
  • Geo-reference the grid (was done manually in this study) to overlap with the watershed of the study area.
  • Right click the MODFLOW grid shapefile, select “Joins and Relates,” and then “Join”.
  • Select “Join data from another layer based on spatial location.”
  • Choose the point feature class as the joining layer.
  • 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).
  • Select the cells you are interested in by filtering based on the “Distance” filed.
  • 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.
  • Open the “output.hru” file for the calibrated SWAT model in Microsoft Excel.
  • Perform the Text to Columns operation if necessary (to separate data by column).
  • 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.
  • Delete the yearly average rows in each column.
  • 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.
  • Organize columns in the following format:
    Water 11 00110 i002
  • Insert the spreadsheet into the middle of the MODFLOW VMP file.
Table A4. List of definitions for all parameters used in the SWAT-CUP calibration process (definitions based on Abbaspour [35]). “sno” parameters are only applied to sub-basin 24. All other parameters were repeated with unique values for sub-basin 24. “v” denotes exact value changes, whereas “r” denotes relative changes (by a percentage) to the parameter value of interest. The “–” in some cells indicates that no value was used.
Table A4. List of definitions for all parameters used in the SWAT-CUP calibration process (definitions based on Abbaspour [35]). “sno” parameters are only applied to sub-basin 24. All other parameters were repeated with unique values for sub-basin 24. “v” denotes exact value changes, whereas “r” denotes relative changes (by a percentage) to the parameter value of interest. The “–” in some cells indicates that no value was used.
ParameterDefinitionInitialFinalFinal (Sub 24)
v__PLAPS.subPrecipitation change with increased elevation (mm H2O/km).--4.85--
v_TLAPS.subTemperature change with increased elevation (degrees C/km).--−13.65--
v__TMPINC().subMax/min temperature adjustment (degrees C)--−0.01575−0.25925
r__CN2.mgtSoil curve number (function of soil permeability)0.3−0.20125−0.44985
r__ALPHA_BF.gwBaseflow alpha factor. Gauges groundwater response to recharge changes (1/days)0.10.14050.312
v__REVAPMN.gwThreshold water depth in shallow aquifer required for deep percolation to occur (mm H2O)487.5238.75852.575012
v__GW_DELAY.gwEstimate for delay time312.50.158225 (r)1.7115
r__GW_REVAP.gwTransfer rate between shallow aquifer and root zone, ranging from 0.02 (restricted) to 0.20 (unrestricted).0.15050.0224880.023858
v__GWQMN.gwThreshold water depth in shallow aquifer required for return flow to occur (mm H2O)18750.195 (r)343.75
v__RCHRG_DP.gwFraction of deep water percolation from root zone (0.0 to 1.0)0.4750.124325 (r)0.147938
v__CH_N2.rteManning’s roughness coefficient (“n”) for a channel (0.025 to 0.150)0.183750.005 (r)0.09875
v__CH_K2.rteEffective hydraulic conductivity of the channel walls (mm/h)12.490250.04750.00325
v__ESCO.hruCompensation factor for soil evaporation ranging from 0.01 (more evaporative demand from deeper soils) to 1 (less deep evaporative extraction)0.3250.7911250.785125
v__EPCO.hruCompensation factor for plant uptake ranging from 0.01 (less water uptake from deeper soil) to 1 (more water uptake from deeper soil)0.7250.8270.75125
r__OV_N.hruManning’s roughness coefficient (“n”) for overland flow (0.010 to 0.480)−0.38−0.0750.2275
v__HRU_SLP.hruAverage slope steepness (y/x, or m/m).--0.326250.0625
v_SLSUBBSN.hruAverage slope length (m), measured at the point that flow starts to concentrate. Generally no more than 90 m.--66.44999731.4
v__SOL_AWC().solAvailable water capacity of soil layer (mm H2O/mm soil)0.180.2330.333425
r__SOL_K().solSoil saturated hydraulic conductivity (mm/h)−0.22−0.16465−0.137675
r__SOL_BD().solSoil moist bulk density (g/cm3 or Mg/m3). Generally between 1.1 and 1.9 g/cm3−0.3−0.08575−0.00035
r__SOL_ALB().solSoil albedo (moist). Ratio of sunlight reflected to sunlight that hits the soil.0.06−0.1250.095
v__SUB_SFTMP().snoAverage temperature at which snowfall occurs (degrees C)−19--−0.2425
v__SUB_SMTMP().snoThreshold temperature for snow melt to occur (degrees C)3--0.9315
r__SUB_SMFMX().snoSnow melt factor on June 21 (mm H2O/degrees C-day)−0.14--−0.135
r__SUB_SMFMN().snoSnow melt factor on December 21 (mm H2O/degrees C-day)−0.14--0.055
v__SUB_TIMP().snoLag factor of the snow pack temperature (influenced by previous day’s snow pack temperature)0.425--0.135322
Figure A2. Calibration and validation results for the Little Smoky River watershed tributary gauge.
Figure A2. Calibration and validation results for the Little Smoky River watershed tributary gauge.
Water 11 00110 g0a2
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.
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.
Water 11 00110 g0a3
Table A5. 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]).
NameContent
swatmf_out_MF_gwswGW–SW exchange (m3/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_rechargeRecharge (m3/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_riverstageMODFLOW river stage (m) for each river cell, listed in order of the river cell ID.
swatmf_out_SWAT_channelChannel 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_gwswExchange (m3/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_rechargeDeep 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. Comparison of simulated SWAT–MODFLOW results to observed values at the Little Smoky River watershed tributary gauge.
Figure A4. Comparison of simulated SWAT–MODFLOW results to observed values at the Little Smoky River watershed tributary gauge.
Water 11 00110 g0a4

Appendix B. Supplemental Information

Appendix B.1. MODFLOW Model Construction

Upon completion of the creation process, the river cells can then be imported into MODFLOW. The row and column of the river cells, as well as its stage, depth, thickness of the riverbed, and conductance, are all data required for a successful import. The RCH data was copied from the SWAT model “output.hru” file, which contains all output values for the entire model. The data brought the recharge schedule of MODFLOW exactly in line with that of the SWAT model. Similar to the process of creating a .RIV file for the MODFLOW model, multiple steps were required before importing the data, which are fully described in the Appendix (Workflow 2). Although the recharge data is ignored when the SWAT and MODFLOW models are coupled [4], it was included in the MODFLOW model as a test of its ability to handle the variable fluxes that would occur in SWAT–MODFLOW.
To allow for a full run of the MODFLOW model (and, by extension, the SWAT–MODFLOW model), some controls under the “Run” menu of Visual MODFLOW had to be changed from their default values. In particular, the “Maximum Outer Iterations” value, found under “Solver,” had to be increased to allow for more than 100 iterations per time step. Although the most outer iterations required for any single time step was 112 at the 2652nd simulation day (with the vast majority of time steps far under this value), the “Maximum Outer Iterations” value was changed to 1000 to account for possible increases to the iteration number when the model was coupled with SWAT.

Appendix B.2. SWAT–MODFLOW Execution

When all required files are compiled into one folder, the SWAT–MODFLOW program can be run. This study uses Bailey’s “SWAT_MODFLOW6.exe” executable, which is an updated version of the original code that was made available after the publication of Bailey et al. [4]. Although version 2 of this code is publicly available on the SWAT website (https://swat.tamu.edu/software/SWAT–MODFLOW/, including a tutorial), the author of the code will tailor the executable to fit the needs of the project, such as showing headers for each category (e.g., river cell location and flux rate, as in this study). Depending heavily on the complexity and simulation times of both component models, this SWAT–MODFLOW code may take anywhere from minutes to many hours to run to completion.
Many of the MODFLOW input files must be modified before attempting to run SWAT–MODFLOW. Each file has an “identifier number,” which helps the code to access it when required. However, some of these identifier numbers overlap with those of the SWAT files, therefore necessitating a change. 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].
Table A6. Summary statistics for the MODFLOW component.
Table A6. Summary statistics for the MODFLOW component.
MODFLOW Summary Statistics
Number of Data Points379
Max. Residual (m)−155.379
Min. Residual (m)0.134
Residual Mean (m)−21.349
Normalized RMS (%)6.506
Correlation Coefficient0.956

References

  1. St. Jacques, J.-M.; Sauchyn, D.J.; Zhao, Y. Northern Rocky Mountain streamflow records: Global warming trends, human impacts or natural variability? Geophys. Res. Lett. 2010, 37, L06407. [Google Scholar] [CrossRef]
  2. Waughray, D. Water Security—The Water-Food-Energy-Climate Nexus. World Economic Forum Water Initiative; Island Press: Washington, DC, USA, 2011. Available online: http://www3.weforum.org/docs/WEF_WI_WaterSecurity_WaterFoodEnergyClimateNexus_2011.pdf (accessed on 24 April 2018).
  3. Faramarzi, M.; Srinivasan, R.; Iravani, M.; Bladon, K.D.; Abbaspour, K.C.; Zehnder, A.J.B.; Goss, G.G. Setting up a hydrological model of Alberta: Data discrimination analyses prior to calibration. Environ. Model. Softw. 2015, 74, 48–65. [Google Scholar] [CrossRef]
  4. Bailey, R.T.; Wible, T.C.; Arabi, M.; Records, R.M.; Ditty, J. Assessing regional-scale temporal patterns of groundwater-surface water interactions using a coupled SWAT-MODFLOW model. Hydrol. Process. 2016, 30, 4420–4433. [Google Scholar] [CrossRef]
  5. Masud, M.B.; McAllister, T.; Cordeiro, M.R.C.; Faramarzi, M. Modeling future water footprint of barley production in Alberta, Canada: Implications for water use and yields to 2064. Sci. Total Environ. 2018, 616–617, 208–222. [Google Scholar] [CrossRef] [PubMed]
  6. Sophocleous, M.A.; Koelliker, J.K.; Govindaraju, R.S.; Birdie, T.; Ramireddygari, S.R.; Perkins, S.P. Integrated numerical modeling for basin-wide water management: The case of the Rattlesnake Creek basin in south-central Kansas. J. Hydrol. 1999, 214, 179–196. [Google Scholar] [CrossRef]
  7. Markstrom, S.L.; Niswonger, R.G.; Regan, R.S.; Prudic, D.E.; Barlow, P.M. GSFLOW—Coupled Ground-Water and Surface-Water Flow Model Based on the Integration of the Precipitation-Runoff Modeling System (PRMS) and the Modular Ground-Water Flow Model (MODFLOW-2005). In U.S. Geological Survey Techniques and Methods 6-D1; Geological Survey: Reston, VA, USA, 2008; 254p. [Google Scholar]
  8. Sudicky, E.A.; Jones, J.P.; Park, Y.J.; Brookfield, A.E.; Colautti, D. Simulating complex flow and transport dynamics in an integrated surface-subsurface modeling framework. Geosci. J. 2008, 12, 107–122. [Google Scholar] [CrossRef]
  9. Kollet, S.; Sulis, M.; Maxwell, R.M.; Paniconi, C.; Putti, M.; Bertoldi, G.; Coon, E.T.; Cordano, E.; Endrizzi, S.; Kikinzon, E.; et al. The integrated hydrologic model intercomparison project: A second set of benchmark results to diagnose integrated hydrology and feedbacks. Water Resour. Res. 2016, 53, 867–890. [Google Scholar] [CrossRef]
  10. Abbott, M.B.; Bathurst, J.C.; Cunge, J.A.; O’Connell, P.E.; Rasmussen, J. An introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”, 1: History and philosophy of a physically-based, distributed modelling system. J. Hydrol. 1986, 87, 45–59. [Google Scholar] [CrossRef]
  11. Kim, N.W.; Chung, I.M.; Won, Y.S.; Arnold, J.G. Development and application of the integrated SWAT-MODFLOW model. J. Hydrol. 2008, 356, 1–16. [Google Scholar] [CrossRef]
  12. Guzman, J.A.; Moriasi, D.N.; Gowda, P.H.; Steiner, J.L.; Starks, P.J.; Arnold, J.G.; Srinivasan, R. A model integration framework for linking SWAT and MODFLOW. Environ. Model. Softw. 2015, 73, 103–116. [Google Scholar] [CrossRef]
  13. Arnold, J.G.; Srinivasan, R.; Muttiah, R.S.; Williams, J.R. Large area hydrologic modeling and assessment Part 1: Model development. J. Am. Water Resour. Assoc. 1998, 34, 73–89. [Google Scholar] [CrossRef]
  14. McDonald, M.G.; Harbaugh, A.W. A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model; USGS Open-File Report 83-875; US Geological Survey: Reston, VA, USA, 1983.
  15. Kim, N.W.; Chung, I.M.; Won, Y.S. The development of fully coupled SWAT-MODFLOW model (I) model development. J. Korea Water Resourc. Assoc. 2004, 37, 499–507. [Google Scholar] [CrossRef]
  16. Kim, N.W.; Chung, I.M.; Won, Y.S. The development of fully coupled SWAT-MODFLOW model (II) evaluation of model. J. Korea Water Resourc. Assoc. 2004, 37, 509–515. [Google Scholar] [CrossRef]
  17. Chevron Canada. Our History in Fox Creek. 2017. Available online: http://www.chevron.ca/our-businesses/kaybob-duvernay-appraisal-program/our-history-in-fox-creek (accessed on 20 August 2017).
  18. Petrel Robertson Consulting Ltd. Integrated Assessment of Water Resources for Unconventional Oil and Gas Plays, West-Central Alberta: Aquifers in Shallow Bedrock and Surficial Sediments, Final (Year 2) Report; Petroleum Technology Alliance of Canada: Calgary, AB, Canada, 2014; 32p. [Google Scholar]
  19. Johnson, E.G.; Johnson, L.A. Hydraulic Fracture Water Usage in Northeast British Columbia: Locations, Volumes and Trends. In Geoscience Reports 2012; British Columbia Ministry of Energy and Mines: Victoria, BC, Canada, 2012; pp. 41–63. [Google Scholar]
  20. Canadian Society for Unconventional Gas. Understanding Hydraulic Fracturing. Information about Canada’s Emerging Energy Resources. 2016. Available online: https://www.scribd.com/document/203055543/CSUG-Hydraulic Frac-Brochure (accessed on 18 October 2017).
  21. Alessi, D.S.; Zolfaghari, A.; Kletke, S.; Gehman, J.; Allen, D.M.; Goss, G.G. Comparative analysis on hydraulic fracturing wastewater practices in unconventional shale development: Water sourcing, treatment and disposal practices. Can. Water Resour. J. 2017, 42, 105–121. [Google Scholar] [CrossRef]
  22. Hughes, A.T.; Smerdon, B.D.; Alessi, D.S. Hydraulic properties of the Paskapoo Formation in west-central Alberta. Can. J. Earth Sci. 2017, 54, 883–892. [Google Scholar] [CrossRef] [Green Version]
  23. Scanlon, B.R.; Reedy, R.C.; Male, F.; Hove, M. Managing the Increasing Water Footprint of Hydraulic Fracturing in the Bakken Play, United States. Environ. Sci. Technol. 2016, 50, 10273–10281. [Google Scholar] [CrossRef] [PubMed]
  24. Grasby, S.E.; Chen, Z.; Hamblin, A.P.; Wozniak, P.R.; Sweet, A.R. Regional characterization of the Paskapoo bedrock aquifer system, southern Alberta. Can. J. Earth Sci. 2008, 45, 1502–1516. [Google Scholar] [CrossRef]
  25. Alberta Water Portal Society. The history of climate in Alberta and effects of climate change on Alberta’s watersheds. In History and Effects of Climate Change on Alberta’s Watersheds; Calgary, AB, Canada, 2017; Available online: https://albertawater.com/history-and-effects-of-climate-change-on-alberta-s-watersheds (accessed on 16 March 2018).
  26. Wright, G.N. (Ed.) The Western Canada Sedimentary Basin, a Series of Geological Sections Illustrating Basin Stratigraphy and Structure; Canadian Society of Petroleum Geologists and the Geological Association of Canada: Calgary, AB, Canada, 1984. [Google Scholar]
  27. Bachu, S.; Michael, K. Hydrogeology and Stress Regime of the Upper Cretaceous-Tertiary Coal-Bearing Strata in Alberta; EUB/AGS Earth Sciences Report 2002-04; Alberta Energy and Utilities Board: Edmonton, AB, Canada, 2002. [Google Scholar]
  28. Smerdon, B.D.; Atkinson, L.A.; Hartman, G.M.D.; Playter, T.L.; Andriashek, L.D. Field Evidence of Nested Groundwater Flow along the Little Smoky River, west-central Alberta; AER/AGS Open File Report 2016-02; Alberta Energy Regulator: Edmonton, AB, Canada, 2016; Available online: https://ags.aer.ca/publications/OFR_2016_02.html (accessed on 22 September 2018).
  29. Dawson, F.M. Uppermost Cretaceous and Tertiary Strata; in Geological Atlas of the Western Canada Sedimentary Basin; Mossop, G.D., Shetsen, I., Eds.; Canadian Society of Petroleum Geologists and Alberta Research Council: Calgary, AB, Canada, 1994. Available online: http://www.ags.gov.ab.ca/publications/wcsb_atlas/atlas.html (accessed on 24 November 2017).
  30. Lyster, S.; Andriashek, L.D. Geostatistical Rendering of the Architecture of Hydrostratigraphic Units within the Paskapoo Formation, Central Alberta; ERCB/AGS Bulletin 66; Energy Resources Conservation Board: Edmonton, AB, Canada, 2012. Available online: http://www.ags.gov.ab.ca/publications/abstracts/BUL_066.html (accessed on 16 November 2017).
  31. Jarvis, A.; Reuter, H.I.; Nelson, A.; Guevara, E. Hole-Filled SRTM for the Globe Version 4, the CGIAR-CSI SRTM 90m Database. 2008. Available online: http://srtm.csi.cgiar.org (accessed on 26 July 2018).
  32. Cordeiro, M.R.C.; Lelyk, G.; Krobel, R.; Legesse, G.; Faramarzi, M.; Masud, M.B.; McAllister, T. Deriving a dataset for agriculturally relevant soils from the Soil Lanscapes of Canada (SLC) database for use in Soil and Water Assessment Tool (SWAT) simulations. Earth Syst. Sci. Data 2018, 10, 1673–1686. [Google Scholar] [CrossRef]
  33. Niswonger, R.G.; Panday, S.; Ibaraki, M. MODFLOW-NWT, A Newton formulation for MODFLOW-2005. US Geol. Surv. Tech. Methods 2011, 6, 44. [Google Scholar] [CrossRef]
  34. Park, S.; Bailey, R.T. SWAT-MODFLOW Tutorial—Documentation for Preparing Model Simulations; Department of Civil and Environmental Engineering, Colorado State University: Fort Collins, CO, USA, 2017; 56p. [Google Scholar]
  35. Abbaspour, K. User Manual for SWAT-CUP: SWAT Calibration and Uncertainty Analysis Programs; Swiss Fed. Inst. of Aquatic Science and Technology: Dübendorf, Switzerland, 2015; Available online: http://www.eawag.ch/organization/abteilungen/siam/software/swat/index_EN/ (accessed on 5 March 2017).
  36. Krause, P.; Boyle, D.P.; Base, F. Comparison of different efficiency criteria for hydrological model assessment. Adv. Geosci. 2005, 5, 89–97. [Google Scholar] [CrossRef] [Green Version]
  37. Faramarzi, M.; Abbaspour, K.C.; Schulin, R.; Yang, H. Modelling Blue and Green Water Resources Availability in Iran. Hydrol. Process. 2009, 23, 486–501. [Google Scholar] [CrossRef]
  38. Arnold, J.G.; Kiniry, J.R.; Srinivasan, R.; Williams, J.R.; Haney, E.B.; Neitsch, S.L. Soil & Water Assessment Tool—Input/Output Documentation; Version 2012; Texas Water Resources Institute: College Station, TX, USA, 2012; 650p. [Google Scholar]
  39. Brinsky, J. Basal Belly River Hydraulic Head and Show Map; Encana Corporation: Calgary, AB, Canada, 2014. [Google Scholar]
  40. IPCC. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the International Panel on Climate Change; Core Writing Team, Pachauri, R.K., Meyer, L.A., Eds.; IPCC: Geneva, Switzerland, 2014; 151p. [Google Scholar]
  41. Chen, J.; Brissette, F.P.; Leconte, R. Uncertainty of downscaling method in quantifying the impact of climate change on hydrology. J. Hydrol. 2011, 401, 190–202. [Google Scholar] [CrossRef]
  42. Van Vuuren, D.P.; Edmonds, J.; Kainuma, M.; Riahi, K.; Thomson, A.; Hibbard, K.; Hurtt, G.C.; Kram, T.; Frey, V.; Lamarque, J.-F.; et al. The representative concentration pathways: An overview. Clim. Chang. 2011, 109, 5–31. [Google Scholar] [CrossRef]
  43. Smerdon, B.D. A synopsis of climate change effects on groundwater recharge. J. Hydrol. 2017, 555, 125–128. [Google Scholar] [CrossRef]
  44. Pacific Climate Impacts Consortium. Statistically Downscaled Climate Change Scenarios. 2014. Available online: https://www.pacificclimate.org/data/statistically-downscaled-climate-scenarios (accessed on 25 June 2018).
  45. Faramarzi, M.; Abbaspour, K.; Adamowicz, W.L.; Lu, W.; Fennell, J.; Zehnder, A.J.B.; Goss, G. Uncertainty based assessment of dynamic freshwater scarcity in semi-arid watersheds of Alberta, Canada. J. Hydrol. Reg. Stud. 2017, 9, 48–68. [Google Scholar] [CrossRef]
  46. Beets, L. SWAT-MODFLOW Documentation; Word Document; Wageningen University & Research: Wageningen, Netherlands, 2016. [Google Scholar]
  47. Hyndman, R.J.; Koehler, A.B. Another look at measures of forecast accuracy. Int. J. Forecast. 2006, 22, 679–688. [Google Scholar] [CrossRef] [Green Version]
  48. Pearson, K. Notes on regression and inheritance in the case of two parents. Proc. R. Soc. Lond. 1895, 58, 240–242. [Google Scholar] [CrossRef]
Figure 1. The study area (Little Smoky River Watershed), shown within the province of Alberta, Canada.
Figure 1. The study area (Little Smoky River Watershed), shown within the province of Alberta, Canada.
Water 11 00110 g001
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. 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.
Water 11 00110 g002
Figure 3. South-north cross-section of the MODFLOW modelled area, displaying the geological units present.
Figure 3. South-north cross-section of the MODFLOW modelled area, displaying the geological units present.
Water 11 00110 g003
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. 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.
Water 11 00110 g004
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. 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.
Water 11 00110 g005
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 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.
Water 11 00110 g006
Figure 7. Comparison of the monthly simulated vs. observed streamflow for the watershed outlet (Outlet Gauge in Figure 3).
Figure 7. Comparison of the monthly simulated vs. observed streamflow for the watershed outlet (Outlet Gauge in Figure 3).
Water 11 00110 g007
Figure 8. Long-term (1986–2007) average monthly groundwater–surface water (GW–SW) exchange differences from the average annual (a) for each MODFLOW river cell, in the months of January (b), June (c), and October (d). Dev (m3/d) = monthly (m3/d) − yearly (m3/d).
Figure 8. Long-term (1986–2007) average monthly groundwater–surface water (GW–SW) exchange differences from the average annual (a) for each MODFLOW river cell, in the months of January (b), June (c), and October (d). Dev (m3/d) = monthly (m3/d) − yearly (m3/d).
Water 11 00110 g008
Figure 9. Long-term (1986–2007) average annual GW–SW exchange data (a), and differences from the average annual data for each MODFLOW river cell, with a pumping rate of 468 m3/d (b), the averaged RCP8.5 (the heaviest greenhouse gas emission) scenario for 2010–2034 (c), and the averaged RCP8.5 scenario with a pumping rate of 4680 m3/d (d). Dev (m3/d) = monthly (m3/d) − yearly (m3/d).
Figure 9. Long-term (1986–2007) average annual GW–SW exchange data (a), and differences from the average annual data for each MODFLOW river cell, with a pumping rate of 468 m3/d (b), the averaged RCP8.5 (the heaviest greenhouse gas emission) scenario for 2010–2034 (c), and the averaged RCP8.5 scenario with a pumping rate of 4680 m3/d (d). Dev (m3/d) = monthly (m3/d) − yearly (m3/d).
Water 11 00110 g009
Table 1. General information of selected General Circulation Models (GCMs) and scenarios available from the Pacific Climate Impacts Consortium (PCIC) [44].
Table 1. General information of selected General Circulation Models (GCMs) and scenarios available from the Pacific Climate Impacts Consortium (PCIC) [44].
AcronymCountrySourceSpatial Resolution (Long × Lat, Degrees)
CanESM2 *CanadaCanadian Centre for Climate Modelling and Analysis2.81 × 2.79
CCSM4 *United StatesNational Center for Atmospheric Research1.25 × 0.9
CNRM-CM5 *FranceCentre National de Recherches Meteorologiques/Centre Europeen de Recherche et Formation Avancees en Calcul Scientifique1.41 × 1.40
CSIRO-MK3 *AustraliaCommonwealth Scientific and Industrial Research Organization in Collaboration with the Queensland Climate Change Centre of Excellence1.875 × 1.86
MIROC5 *JapanMeteorological Research Institute1.41 × 1.39
* 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.
Hydrometric StationDrainage Area (km2)p-Factorr-FactorR2Initial R2bR2Initial bR2NSE
Outlet Gauge 11,1400.440.650.670.290.62820.05210.53
Tributary Gauge10410.520.850.580.460.56850.3774--
Average (watershed)--0.480.750.620.370.590.21--

Share and Cite

MDPI and ACS Style

Chunn, D.; Faramarzi, M.; Smerdon, B.; Alessi, D.S. 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. Water 2019, 11, 110. https://doi.org/10.3390/w11010110

AMA Style

Chunn D, Faramarzi M, Smerdon B, Alessi DS. 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. Water. 2019; 11(1):110. https://doi.org/10.3390/w11010110

Chicago/Turabian Style

Chunn, David, Monireh Faramarzi, Brian Smerdon, and Daniel S. Alessi. 2019. "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" Water 11, no. 1: 110. https://doi.org/10.3390/w11010110

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop