Abstract
Watersheds located within a mountain to coast physiographic setting have been described as having a highly interconnected surface water and groundwater environment. The quantification of groundwater—surface water interactions at the watershed scale requires upscaling. This study uses MIKE SHE, a coupled numerical model, to explore the seasonally and spatially dynamic nature of these interactions in the Cowichan Watershed on Vancouver Island, British Columbia, Canada. The calibrated model simulates a transition of the Cowichan River from mostly gaining within the valley, to losing stream near the coast where groundwater extraction is focused. Losing and gaining sections correlate with geological substrate. Recharge across the watershed accounts for 17% of precipitation. Climate change is projected to lessen snowpack accumulation in the high alpine and alter timing of snowmelt, resulting in higher spring and winter river discharge and lower summer flows.
1. Introduction
Watersheds located within a mountain-to-coast physiographic setting are unique in that they have been described as having a highly connected surface water and groundwater environment [1]. The high degree of coarse alluvial material, coupled with a steep topographic setting, creates conditions whereby the surface water and groundwater systems strongly interact. In regions where the climate is seasonally dry, the principal source of water within a stream is often from the discharge of groundwater [1–3]. Streams, however, may also recharge the aquifer, particularly during the freshet (e.g., [4]). These relationships are often poorly understood aspects of the hydrology within a mountainous watershed. Water balances, including estimates of recharge and discharge, are also highly variable within this type of setting, especially since the climate gradient (heavy precipitation in the mountains to relatively low precipitation near the ocean) is both seasonally and spatially variable. Also, there is a high degree of geological variability (shallow or exposed bedrock near the crest of the valley, and alluvium of variable thickness and composition within the valley). Management of water in such watersheds thus requires sound understanding of a range of hydrologic processes and particularly those factors that influence the interaction of groundwater and surface water at a range of spatial and temporal scales [1, 2, 5].
Coupled groundwater–surface water models are being increasingly used to examine a variety of environmental interactions, including interactions in small mountainous catchments [6]; solute transport [7, 8], flood wave modeling [9], catchment water resource management and understanding [10–17], the role of aquifer heterogeneity [18, 19], and wetland alterations and aquaculture [20, 21].
This study aims to contribute to the knowledge of surface water and groundwater interactions within a mountain-to-coast watershed, specifically through investigating how these interactions may be influenced by different stressors within the watershed. The study area is the Cowichan Watershed, situated on Vancouver Island in British Columbia, Canada (Figure 1). The steep topographical setting and geological conditions create conditions whereby the groundwater and surface water systems are dynamically coupled. The watershed is comprised of several catchments, covers an area of approximately 930 km2, attains a maximum elevation of approximately 1483 metres above mean sea level (masl) in the headwater region to the west, and terminates at sea elevation near its eastern extent. Cowichan Lake has a surface area of 62 km2 and stretches nearly 31 km from west to east. The Cowichan River flows from the headwaters at Cowichan Lake eastward for nearly 45 km to the estuary in Cowichan Bay near Duncan.
The watershed itself is a vast valley with a large accumulation of valley fill sediments, flanked by valley walls with thin veneers of soil. The climate is temperate with cool and wet fall and winter seasons, while the spring and summer months are warm and typically much drier. There is a strong precipitation gradient (decreasing to the east) due to a rain shadow effect. The lower coastal section of the watershed receives half the amount of precipitation (~1000 mm/year) than that received at Lake Cowichan (~2000 mm/year). It is estimated that the mountainous regions at the western boundary of this watershed can receive up to 4500 mm of precipitation annually [22]. Most precipitation occurs during the winter months, while very low amounts are measured in the summer months. At most, snow accounts for ~5–15% of the total precipitation.
The watershed provides fresh water to over 43,000 people, as well as agriculture and several forms of industry (fish farms, paper mills, etc.). More than 530 surface water licenses have been issued to divert water from streams and lakes in the watershed, and more than 1,300 wells have been drilled to pump water from the aquifers [23]. Water users within the watershed include agriculture, industry, and urban and rural water supply. There are several large water users within the watershed including, the pulp and paper industry (Catalyst Paper), fish hatcheries, and municipal water supply (Figure 2).
Recently, seasonal fluctuations, and changes in the timing of rainfall events have created challenges for managing water in the watershed. The variability in seasonal rainfall is extremely large; flooding conditions can occur in the winter, while drought conditions can prevail in the summer. Water demand puts added stress onto the hydrologic system, as peak demands for water occur during the summer low flow season. In 2012, the seriousness of the problem became obvious as returning salmon struggled to reach their spawning locations; this gained attention from the press. To address these issues and gain insight on the hydrologic conditions, a coupled groundwater-surface water model was developed using MIKE SHE [24]. The calibrated model is used to assess groundwater recharge and discharge, estimate the contributions of groundwater to the surface water system, identify key gaining portions of the Cowichan River, and evaluate the impact of localized pumping on the system. Lastly, the model is used to project how future climate may affect the dynamics of the hydrogeological system (over the next 40 and 70 years).
2. Materials and Methods
2.1. The MIKE SHE Modeling Interface
MIKE SHE is a deterministic and distributed modeling system that uses finite difference representations in mass and energy and measured empirical relationships to simulate aspects of the hydrologic cycle [24]. At its core is a framework of modules that are used to simulate the following processes: interception and evaporation, overland flow, unsaturated zone flow, saturated zone flows, and water quality. Rivers, lakes, and other channels are operated in the one-dimensional model, MIKE 11, which is coupled directly to the MIKE SHE model. The interception and evaporation module computes the actual evaporation (AET) from an area using user-defined potential evaporation (PET), using the Kristensen and Jensen [25] model. This model requires vegetation dependent parameters such as leaf area index (LAI), root characteristics, and an interception parameter. Unsaturated flow is calculated in 1D, vertically. A soil moisture retention curve along with the saturated hydraulic conductivity are defined for each soil class spatially, and the vertical layer. Richards’ equation is solved, and water flows from the unsaturated zone to the saturated zone, or vice versa. The overland flow component simulates runoff when the infiltration capacity of the soil is exceeded, when groundwater discharges to the surface, or when streams flood their banks. In this study, the flow solution utilized the diffusive wave approximation of the Saint-Venant equation, whereby topography, and the Manning’s coefficient control the direction and rate of runoff, respectively. The saturated zone flow component in MIKE SHE is 3D and is based on Darcy’s equation. Boundary conditions such as fixed head, zero flux, gradient, and specified flux are options which control the flow of groundwater within the model. Subsurface conditions are modelled as layers and lenses, with representative hydraulic properties assigned.
As mentioned, MIKE 11 controls the routing of water in rivers and lakes. The rivers module comprises four main components: the river network, river cross-sections, boundary conditions, and hydrodynamic parameters. MIKE 11 solves channel flow through the use of a 1D St. Venant equation based on the complete dynamic wave formulation [26]. MIKE SHE and MIKE 11 are coupled through the use of river links (-points). During a simulation, the amount of water entering or exiting the linking cells is calculated based on Darcy’s equation. Lateral inflows and outflow from overland flow as well as river-aquifer exchanges are completed for each computational time step [24].
2.2. Model Setup
The simulation period spanned January 1, 1998, to December 31, 2012 (the most recent date of available data). It was important to capture the year 2012, as the motivation for this study was the anomalously low river discharge of the Cowichan River during the salmon spawning season and the timing of late summer rains in that year. The initial groundwater levels were assigned at ground surface, as discussed below; therefore, the model had to spin up to achieve a dynamically stable state in which the deep groundwater levels no longer decline over time. It was found that starting the model in 1998 was adequate for both spin up and providing a suitable time frame (12 years) for analysis. All data output time steps were set to 24 hours with the exception of groundwater, which was set to 48 hours.
The model grid size was 200 m by 200 m. Topography was assigned using a 200 m digital elevation model (DEM). The model boundary conditions consist of a zero flux boundary to represent the topographical boundary of the watershed, and a specified head (sea level) within the alluvial layer where the model meets the ocean and Cowichan Bay. Underneath the alluvial layer, the bedrock layer is set to a zero flux boundary. These boundary conditions attempt to mimic groundwater discharge in a coastal environment, whereby deep groundwater is directed upward when it intersects the freshwater-saltwater interface. Thus, any discharge from the bedrock will be directed upward to the surficial sediments and subsequently out of the model. Overall, the assigned boundary conditions route whatever precipitation falls onto the model domain out of the model along three potential pathways: evaporation, surface water termination at the ocean, and groundwater discharge upward along the coast and directly into the ocean.
2.3. Meteorological Data
Precipitation was imported into MIKE SHE using a “station based” time varying format. Throughout the watershed, annual precipitation ranged from approximately 1000 to 5600 mm/yr and was used as a basis to divide the watershed into 10 zones based on an increase of 500 mm/yr per zone (Figure 3). Zones 1 and 2 were represented by the precipitation recorded at the Kelvin Creek climate station, while Zones 3 through 10 were represented by the Forestry Research Station (Figure 3). To model the increase in precipitation due to orographic and rain shadow effects, a proration was applied to each zone, based on the median amount of precipitation observed within that zone as illustrated on Figure 3.
Air temperature was defined using the dataset from the Forestry Research Station as there were negligible differences between this station and the Kelvin Creek Station. The spatial variation in temperature was modified from the station data according to a fixed temperature lapse rate (Figure 4(a)). The temperature recorded at the Jump Creek Snow Pillow Station was used as a calibration for the temperature lapse rate parameter.
(a)
(b)
(c)
Potential evapotranspiration (PET) was calculated using the Penman-Monteith method, carried out using a software package AWSET [27]. Daily climate data consisted of mean air temperature, humidity, solar radiation, and wind speed. Due to the variability of these climate parameters at the watershed scale, PET was estimated by recording the variability of daily mean temperature as it relates to altitude and location (Figure 4(a)), and then mapping the spatial variability of solar radiation (slope and aspect) within the region using the solar radiation analysis tool in ArcGIS [28] (Figure 4(b)). These two spatial datasets were then merged to create a PET zone map consisting of 42 zones (Figure 4(c)). The PET values were adjusted to an appropriate range based on values in an area with similar climate [29, 30]. All 42 permutations were added to MIKE SHE as time series.
Precipitation, PET, and temperature were set input parameters not subject to calibration, while parameters such as temperature lapse rate and snowmelt parameters were adjusted, through calibration, to match observed snowmelt data recorded at the alpine Jump Creek Snow Pillow Station. Table 1 shows the final snowmelt parameters.
2.4. Land Surface Data
The Cowichan region is classified into two main biogeoclimatic zones: Coastal Douglas-Fir and Coastal Western Hemlock (BC Biogeoclimatic Ecosystem Classification (BEC) system). Leaf area index (LAI) was estimated using satellite reflection imagery and calibrated to a published statistical relationship [31]. Observed data were grouped into five land use classes (Figure 5). For example, LAI values between 0 and 1 were grouped into Class 1 and assigned a LAI value of 0, which spatially represents open bodies of water, and urban areas. LAI values ranging from 5 to 10 were grouped into Class 3 and assigned a LAI value of 7.5 to represent a moderate LAI and distinguish previously logged areas. Classes 4 and 5 represent old growth and biologically dense areas, respectively.
These classes were also defined in terms of phenology to represent seasonal variations in vegetation density. Due to the fact that the Landsat satellite image used in the calculation of LAI was recorded in August of the year 2002, it likely represented a maximum LAI value. The BEC system suggests that most of the mature trees within the study area are coniferous (Western Hemlock, Redcedar, Douglas-Fir, lesser Arbutus, and Garry Oak); however, the understory is made up of ferns and shrubs. Therefore, some variation in LAI is expected [32]. Accordingly, Class 3 was assigned the highest seasonal variation as it represents a deforested system and, therefore, would largely comprise quickly inhabiting deciduous trees (Cottonwood, Aspen, Alder, and various shrubs). The seasonal variation for Classes 4 and 5 was set lower, based on the assumption that the forest cover is likely more mature and reflects the descriptions in the BEC. Within MIKE SHE, a consistent rooting depth of 400 mm was applied. Rooting depths are typically dependent on soil fertility and structure; however, Curt [33] suggests that majority of root mass is concentrated within the top 400 mm of soil, regardless of soil quality.
Along with the spatially and time varying LAI dataset, the vegetation module also contains ET parameters, which can be altered based on site specific data. The parameters include canopy interception (value of 0.05 mm), which needs be met before stem flow and ground infiltration can occur, and empirical coefficients labelled C1, C2, and C3, which relate to the Kristensen and Jensen equation used to calculate actual transpiration and soil evaporation. Coefficients C1 and C2 are plant dependent and influence the distribution between soil evaporation and transpiration. These parameters were set to 0.3 and 0.2 mm/d, respectively. The coefficient C3 is soil dependent and controls the release of water at certain matrix potentials and root densities; this parameter set to 20 mm/day. Finally, AROOT controls the fraction of ET extracted as a function of depth, as larger values have a greater range of ET and approach uniformity as the value nears 0. A value of 0.25 m−1 was used.
Overland flow is defined as the portion of runoff that occurs as sheet flow. If rainfall exceeds the infiltration capacity of the soil, water will move horizontally across the surface, being routed by surface topography at a rate that is calculated using the diffusive wave approximation. The resistance to flow overland is controlled by the “roughness” of the land surface, which can be inferred from land use/cover maps. Each land classification or surface then needs to be transformed directly into a number that assigns hydraulic “roughness.” Within MIKE SHE, the Manning’s coefficient (reciprocal of Manning’s ), which is equivalent to the Strickler roughness coefficient, controls the amount of friction and the velocity at which water can move horizontally. The value of is typically in the range of 0.01 (smooth channels) to 0.10 (thickly vegetated channels). There are several literature sources for estimating Manning’s coefficients over a variety of surfaces (e.g., [34–36]), although most values tend to be modified through the calibration process. The initial settings of Manning’s were obtained from Engman [35] as described below.
To represent the land surface within the Cowichan region, a present land use/cover dataset was used (Figure 6(a)) [37]. To represent mountainous streams, which are often very steep and rocky, all streams within the model were converted to points, assigned a high Manning’s coefficient (100 m1/3/s) and merged with the land use dataset. Accordingly, urban, forested, recently forested, agricultural, streams, and so forth, are represented in terms of a Manning’s within the MIKE SHE overland flow module (Figure 6(b)).
(a)
(b)
2.5. Unsaturated Zone Data
The surficial soils dataset was obtained from Liggett and Gilchrist [38], which is a simplified form of the soil type classifications of Jungen [39], based on the soil’s drainage ability (very poorly to well drained). Surficial geology maps [40, 41] illustrate a thin coverage, or veneer, of soil towards the valley side walls. To capture this morphology, the soil map was further defined to include additional underlying strata. Within ArcGIS, the watershed was divided into two main zones, alluvium and bedrock. The soil classification map was imposed on top of the geology layer and a “merge analysis” was performed, adding the underlying geological contact to the soil layer. Therefore, each of the soil classes contains either an “A” (alluvium) or “B” (bedrock) to signify the underlying material. The unsaturated zone geology was then defined vertically within the model. All soil classes, with the exception of “10A and 10B,” were assigned a base of 2 m below ground surface. Below 2 m depth, the underlying material (alluvium or bedrock) was assigned. Classes 10A and 10B were described as thin soil, and therefore, the unsaturated soil depth is 1 m. Within MIKE SHE, each unsaturated zone requires a “to and from” depth, which must be defined for the full range of the unsaturated zone. Therefore, the depths of the underlying units (alluvium and bedrock) were extended to a range deeper than the maximum simulated thickness of the unsaturated zone (up to 60 m). Figure 7 illustrates the spatial pattern of unsaturated soils, as well as an example of the vertical representation along a cross-section.
The soil class properties were initially defined based on the UNSODA unsaturated soil hydraulic database [42]. Texture class was assigned to each soil class (drainage ability) according to values of saturated conductivity () in the literature. Table 2 shows the van Genuchten soil parameters. As mentioned previously, to represent the deep unsaturated zone, additional unsaturated zone materials and properties were included for alluvium and bedrock.
2.6. Saturated Zone Data
The saturated zone is based on a conceptual model of an alluvial valley and bedrock, whereby a deeply incised valley has been infilled with alluvial sediments, while the mountainous upland areas are covered with a thin veneer of unconsolidated material. Therefore, the saturated zone consisted of two geological layers, “Alluvium” and “Bedrock.” The layers’ setup requires both geological layers to be present throughout the model domain and a measurement to the bottom of the unit specified. To represent the “thinning” of the alluvium material outside of the valley (up the mountain sides), the thickness was set to near zero, while the thickness of the alluvium within the valley ranged from 0.1 to 125 m. The underlying bedrock extended to 500 m so that the “active flow zone,” considered to be in the upper 200 m, was fully captured [45–47].
To represent each mapped aquifer in the valley, geological lenses were added to the model. Each lens represents a hydrogeological unit (HGU), which contains a top and bottom elevation and a spatial extent to represent the unit’s limits. The hydraulic properties for each of the designated aquifer HGUs were obtained from a summary of all available pumping and recovery test data [44]. Not all identified aquifers contained a well that had a pumping test completed and, conversely, pumping tests were completed within areas that did not contain a classified aquifer. Where pumping test data were not available, estimates from the literature were used according to material type [48]. The geomeans of hydraulic conductivity and specific storage provided the initial estimates of parameter values for the model. Throughout the calibration process, several of these values were altered based on model performance. Table 3 shows the final parameter values and Figure 8 shows the hydraulic conductivity distribution throughout the watershed.
2.7. MIKE 11 Stream Network and Hydrometric Flow Data
MIKE 11 models lake and river interactions using cross-sections and assigned elevations. The lake and river network was obtained from the BC Watershed Atlas [49]. Lake and rivers were represented in 1D as single line segments, with the extent of the feature defined by the width of the cross-section. Once MIKE 11 is coupled with MIKE SHE, bed topography and the extent of Cowichan Lake are specified in detail (3D).
Cowichan Lake has a major influence on the hydrology of this region; therefore, special care was taken to represent the lake as accurately as possible. The bathymetry of Cowichan Lake was defined by digitizing published bathymetry maps [50]. Cowichan Lake has a surface area of 62 km2, a shoreline distance of approximately 106.78 km, an estimated average volume of 2.5 billion m3, an average depth of approximately 50 m, and a maximum depth of approximately 160 m.
The Cowichan River network was converted to river nodes (-points) and river branches. Discharge and stage levels are calculated at “” and “” points, respectively. The discharge measurements at -points (positioned half way between each -point) are extrapolated from points in between input cross-sections. Stage measurements are calculated at all -points and are determined by the dynamics of flow within the cross-sections. The coupling of between MIKE 11 and MIKE SHE takes place at the -points.
Conductance values that control the flow of water between the river and the groundwater system were estimated solely from the subsurface geology hydraulic conductivity values as shown previously in Figure 8. These values were adjusted during model calibration.
Boundary conditions were assigned to the MIKE 11 river network. A closed boundary was assigned to the upstream end of the network, while an open boundary was assigned to the mouth of the Cowichan River at Cowichan Bay. This coastal open boundary consists of a water level condition, whereby the tide variations observed at Patricia Bay were used as input. The tide varies from 1.79 to 3.02 m over the simulation period.
The lake stage was initially specified at elevation of 160 masl, with a global bed resistance Manning number of 30. All other parameters (wind factor, computation scheme, and computation parameters) were set to MIKE 11 default values.
2.8. Groundwater and Surface Water Extraction
To model the influence that large water users have on the groundwater and surface water levels within the Cowichan, the estimated extraction rates were added to the model. Six large groundwater users and one large surface water user were included (see Figure 2). The majority of the pumping occurs near the City of Duncan (see Figure 1) clustered around the lower reaches of the Cowichan River. Most groundwater extraction values were provided as monthly totals and, therefore, were modified to be a constant pumping rate in cubic metres per second (m3/s). The groundwater extraction rates for the municipal wells peak during the summer season, nearly doubling relative to the other seasons. The hatcheries generally have an opposite withdrawal schedule, with extraction rates doubling in the winter season compared to the summer season. The identified small and medium water groundwater users were not included in the model as most represent single domestic wells. These small domestic users of groundwater also likely have a septic system on the property (which recycles a large portion of the groundwater back to the subsurface), and therefore, the amount of water lost to the system is thought to be minimal.
Only one large user of surface water was included in the model. Catalyst Paper has an intake on the lower reach of the Cowichan River near Duncan and withdraws water directly from the river. The water leaves the watershed. An annual withdrawal of approximately 50 to 60 million m3 is extracted annually from this location. To model this abstraction, a point-source inflow boundary condition, at the location of the intake, was defined in MIKE 11. The inflow boundary condition was set to a maximum withdraw of −2 m3/s for the entirety of the model simulation. This rate equates to the 63 million m3 of water extracted annually.
2.9. Model Calibration and Validation
The model was setup to run for a simulation period of 14 years (1998–2012). The calibration period for this model was 2002–2010 (8 years). The validation period was from 2010 to 2012. The degree of model fit or calibration was determined by correlation statistics including the mean error (ME), residual mean square error (RMSE), correlation (), and the Nash-Sutcliffe efficiency.
Model calibration first focused on the climatic conditions (snowmelt modeling). Snowmelt calibration consisted of adjusting model parameters including degree day coefficient, temperature lapse rate, and the max wet snow fraction. Each parameter affected the simulated timing (onset and release of snow) and the amount of snow accumulation. Mean daily temperatures measured at the Jump Creek Snow Pillow Station were used to calibrate the temperature lapse. The snow water equivalent (SWE) recorded at the climate station was used to calibrate the amount of water held in snow storage in the alpine regions. The second phase of calibration focused on the hydrometric characteristics. The calibration consisted of adjusting the physical conditions of the MIKE 11 stream system, Manning’s for overland flow and channel flow, and the streambed leakage coefficients. The next phase of calibration included comparing the measured stage and discharge from MIKE 11 to observed lake level and hydrometric data. The Water Survey of Canada (WSC) maintains three stations within watershed, measuring Cowichan Lake levels (08HA009), Cowichan River stage/discharge near the junction of Cowichan Lake to Cowichan River (08HA002), and the stage/discharge of the Cowichan River near Duncan (08HA011). The final phase of calibration focused on the groundwater flow within the region, while still evaluating the hydrometric characteristics. Calibration of groundwater levels within the saturated zone involved adjusting the horizontal and vertical hydraulic conductivity and specific storage. Calibration for groundwater level used the hourly data from the Ministry of Environment (MOE) observation well #204 within Aquifer 186, which has the longest period of record. Historical static groundwater levels measured at the end of drilling of domestic wells (both in alluvium and bedrock) were used to verify the calibration. All other modules (evaporation, soil hydraulic parameters) remained constant during the model simulations due to the limited availability of calibration data.
2.10. Climate Change Simulations
The general consensus from the results of climate modeling in British Columbia indicates that temperatures will generally rise, with the largest increases occurring in the summer [51]. Precipitation is projected to increase in the winter months and decrease in the summer months. These trends are expected to increase atmospheric evaporative demand, decrease snow accumulation, accelerate snowmelt, alter groundwater storage and recharge, alter timing and magnitude of steamflow, and result in a variety of ecological changes [51].
In order to assess how vulnerable the Cowichan Watershed is to the potential impacts of climate change, future climate change data were used to force the MIKE SHE model. Two MIKE SHE simulations were run (one representing the 2050s and one the 2080s). The projected climate change impacts were assessed using the BC Regional Analysis Tool [52]. Specifically, the climate projections from the “TreeGen ensemble” were used [53–55]. The TreeGen downscaling tool was applied to an ensemble of global climate models (GCMs) and SRES AR4 emissions scenarios, with the results compiled for the Province of BC. The results from Canadian Global Coupled Model 3 (CGCM3)-A2 (five model runs) and the Max-Planck Institute for Meteorology (MPI) ECHAM5-A2 (one model run) were used in this study. The A2 emissions scenario was selected because it represents a “worst case” scenario in terms of emissions, CO2 concentrations, and the resulting temperature increase [56].
Several datasets were extracted for the Cowichan area: absolute temperature change (max, min, mean, and medium) and percent change for precipitation and relative humidity for the time periods 2050s (2039–2069) and 2080s (2070–2099). Figure 9 illustrates the absolute change in mean monthly temperature and relative change (as a percent) in monthly precipitation averaged across the study area. Temperature is expected to increase between 1 and 3°C during the period 2050, and by as much as 2–5°C for the 2080s time period (Figure 9, left). The largest temperature differences are expected from July to August and from December to January. (Figure 9, right) indicates that by the 2050s an increase in precipitation of 10–20% is expected for the winter months and a reduction by up to 20% in the summer months. This trend continues throughout the 2080s, increasing by up to 30% in the winter months, and decreasing by 40% in the summer months.
The mean monthly climate shift factors (from the selected models in the ensemble) for each future time period were applied to historical data (1998–2012) from the Cowichan Lake Forestry Research Climate Station and the Kelvin Creek Climate Station. Specifically, the mean monthly climate shifts were applied directly (subtraction or addition to the mean daily temperatures or % increase/decrease to the precipitation rates) to the temporal climate datasets in MIKE SHE. The model was rerun for two 14 years’ period (representing a shift in the 1998–2012 climate data to each of the 2050s and 2080s climate). The baseline climate data and future climate data are discussed in the section on results of the climate change models (see Results and Discussion).
PET was also adjusted for the climate change simulations. The projected minimum, maximum, and mean temperatures, as well as the projected changes in relative humidity were used to calculate new PET values to reflect the projected climate. Again, the AWSET program was utilized to generate daily PET using the Penman-Monteith equation. The shifts to the temperature and humidity were added to the AWSET program by subtracting or adding the absolute temperature change to the min, max, and mean historical daily values, as well as the relative percent change to the historical relative humidity daily values. Modeled solar insolation and wind speed remained the same. By the 2050s, PET is expected to increase by 6.4 to 12.1% and by the 2080s by 11.9 to 21.2% for climate zone 22 (dominant zone in the watershed). The relative shifts in PET closely reflect the projected shifts in temperature. The same relative change in PET (% change) obtained for zone 22 was applied to the daily PET estimations for the other 41 PET permutations.
3. Results and Discussion
3.1. Model Calibration and Validation
The calibration fit statistics are given in Table 4. Figure 10 shows the fit for the two hydrometric stations, and Figure 11 shows the fit for groundwater level. In addition, measurements of groundwater levels made following drilling and reported in the BC WELLS database were compared to the simulated groundwater levels for wells completed in alluvium and bedrock. For the wells completed in alluvium, the was 0.97 and the root mean squared error (RMSE) was 10.78 (). For the bedrock wells, the was 0.70 and the RMSE was 21.5 () with observed groundwater levels being slightly overestimated, likely due to the fact that water levels measured in bedrock wells following drilling have not fully recovered and so would tend to be too low.
(a)
(b)
(c)
3.2. Water Balance
A water balance extraction was performed following calibration. Of interest to this study are the overall exchanges of water between different parts of the model (e.g., between the river and groundwater), the amount of recharge to the saturation zone, and the effect of pumping on the hydrologic system. The total input of water to the model occurs solely as precipitation (100% in input). Water is then partitioned (runoff or overland flow, infiltration or recharge to saturated zone, evaporation) and leaves the model through evaporation, boundary flow from the saturated zone into the ocean, river boundary flow to ocean, surface water extraction, or groundwater extraction, with some water in various stores at any one time (e.g., snow storage, canopy storage overland storage, subsurface storage, etc.).
Table 5 reports the total water balance for the Cowichan Region including error (mm/year). Values are reported for the water year (October 1 to September 30). Recharge is shown in the last column as a separate item. Recharge is computed from the exchange between the unsaturated zone and the saturated zone, and therefore, does not appear in the overall water balance for the watershed.
The water balance results must be examined carefully because there are numerous exchanges that take place. Therefore, the annual percentages do not add up. Overland flow to river (Cowichan River) and ET are the dominant fluxes of water within the Cowichan, constituting 55 and 43% of the annual budget water budget. ET is lost from the watershed; however, overland flow to river may, at other points in the watershed, contribute to groundwater (through the river to baseflow component) and perhaps return to the river downstream (baseflow to river). Thus, these terms are linked and likely elevate the overland flow to river component. The baseflow (groundwater) to river and river to baseflow (groundwater) represent exchange flows between the MIKE SHE and MIKE 11 models. These exchanges take place at the -points.
The water balance results suggest that the Cowichan River is approximately equal in the amount of water the river loses and gains along its length. This relationship is very consistent throughout each water year. The spatial representation of this relationship is explained in detail later. Small negative and positive values are reported for changes in overland flow and snow storage, while 3% of the average annual budget is accounted for by storage changes in the saturated zone. Over the long term, unless the saturated zone is being depleted, this should be zero. The amount of water pumped from the major groundwater users in the lower valley accounts for less than 1% of the total water balance. The average error associated with the convergence of processes in the model was approximately 1% over the calibration and validation periods of the model.
Based on the detailed saturated zone water balance (not shown), annual recharge (determined as the amount of water exchanged from the unsaturated to the saturated zone) is 438 mm/yr, or 17% of the annual precipitation (last column in Table 5). This amount is determined from a yearly average over the calibration and validation period (2002–2012). During this period, the amount of recharge to groundwater varies (253–630 mm/yr) accordingly with yearly variations in precipitation. Taking into account the total variation in precipitation, recharge to groundwater ranges from 14 to 21% of the total annual (WY) precipitation. Hydrogeological studies in close proximity [57–59] have estimated recharge rates to be from 23 to 45% of annual precipitation. However, these recharge rates reflect recharge to individual unconsolidated sand and gravel aquifers, rather than the net recharge across the entire watershed (including low conductivity bedrock).
At a monthly time scale (results not shown), the temporal variation in exchanges with surface water and groundwater mimic closely precipitation variations. Groundwater entering the Cowichan River dominantly occurs from December to May (6-7 mm/month) and is slightly lower from June to November. The exchange from surface water to groundwater follows a similar trend; as one might expect, a higher exchange occurs during the summer when the groundwater table is depressed. Recharge also varies significantly throughout the year. The highest recharge occurs in October and November (>100 mm/month), while a recharge deficit (P-ET) is indicated in the months of June, July and August, with peak deficits at −28 mm/month (loss of water from the saturated zone to the unsaturated zone). This deficit is not only evident in recharge, but also when comparing ET to precipitation over that same time period. May also experiences negative moisture conditions (ET being greater than incoming precipitation); however, recharge is still positive. These results likely reflect the effect of the melting snowpack in the alpine. As the snow melts, it infiltrates the unsaturated zone and eventually reaches the saturated zone.
The year 2012 was particularly bad in terms of sustained discharge within the Cowichan River. Discharge was extremely low, and there was very little precipitation in the later summer months. August and September of 2012 differed the greatest from the average conditions in the Cowichan Watershed. August 2012 had unseasonably low precipitation, resulting in a very large moisture deficit (−102 mm/month) compared to the average of –73 mm/month. This also resulted in greater than 100% reduction in recharge during that month. September was much the same; the moisture deficit in September was −57 mm/month compared to the +21 mm/month average. The climatic variations also caused a recharge deficit in September (−26 mm/month) as compared to the average groundwater recharge of 13 mm/month.
3.3. Recharge and Discharge Areas
Recharge is highly variable across the watershed (Figure 12), which reflects the range of parameters that influence recharge: the rate and annual amount of precipitation, the rate of evapotranspiration, topography and surface roughness coefficients, the hydraulic properties of the unsaturated soil, and, likely most importantly, the depth to the water table from ground surface (unsaturated zone thickness). Areas with a thin soil cover, high amounts of precipitation, and a permeable subsurface material with a groundwater table close to surface will have recharge that is orders of magnitude greater than areas will less precipitation, a low permeability substrate and a deep groundwater table. Figure 12 show a gradient of recharge from west to east. This gradient results primarily from the precipitation patterns within the valley, as yearly precipitation values in the west are several times larger than the east. There are several relatively small circular areas of highly focused recharge. These anomalous areas likely represent topographic depressions in the DEM, where water ponds and infiltrates throughout the simulation.
To assess the accuracy of the discharge features simulated by the model, the location of observed groundwater discharge features, such as springs and wetlands, were superimposed over the simulation results (Figure 12). The simulated linear seepage faces along the northern valley slopes correspond well to observed locations of springs. As well, observed wetland features tend to correspond with low topographic depressions within the lower valley where groundwater discharge occurs. Most discharge features throughout the watershed are situated in the valleys flanked by steep ridges. The discharge occurs as saturated zone to overland exchange.
3.4. Groundwater/Surface Water Interactions
Exchanges at the watershed scale are largely controlled by variations in subsurface lithology, including depth to bedrock and aquifer properties [60, 61]. For example, exchanges that occur in reaches of the Cowichan River where surface water overlies bedrock directly are controlled largely by the hydraulic conductivity of the bedrock, whereas, in other locations, the Cowichan River passes through zones of permeable alluvial deposits where the river channel is deeply incised into the alluvium. Valley width may also affect exchanges [60, 61].
To illustrate the influence of geology on exchanges, the Cowichan River itself (A-A′) was used as a cross-section (Figure 13(a)). This cross-section illustrates the material in contact with the river bed, the depth to bedrock or where bedrock is exposed in the river bed, and the thickness of the alluvial sediments. Also imposed on the figure are the relative positions (at -metres away from the river) of the surface water diversion and groundwater extraction wells to the nearest point of the river. All of the groundwater extraction wells are within unconfined sand and gravel Aquifer 186. Figure 13(b) shows the gaining and losing portions of the river, alongside the geology based on the annual exchanges simulated in 2008. For the majority of the up-river reaches, the Cowichan River is a gaining system (with the exception of a reach from 19500 to 21000 m, near Stoltz Pool). However, further down valley where the relief is lower, the river becomes predominantly losing. Large volumes of water are lost where the river crosses Aquifer 186. Right at the coast, the Cowichan River gains water, as would be expected in a coastal setting due to the presence of the saltwater-freshwater interface at depth, which directs fresh groundwater discharge upwards along a seepage face. As this was a freshwater model, the actual interface was simulated by placing zero flux boundaries in the bedrock and forcing discharge to exit the model domain through the alluvium.
(a)
(b)
Groundwater discharge into the river is highest during the spring season (when groundwater levels are highest) and lowest during the fall (when groundwater levels are low) (Figure 14). Losing conditions are the greatest (most negative) during the winter months when the stage of the river is high and the groundwater table may still be low (due to lag time), resulting in a higher hydraulic gradient. At 44000 m, the exchange conditions shift from predominantly losing to predominantly gaining, but the magnitude of the exchange varies seasonally. When groundwater levels are greater than the river stage (evident in March of 2008), the river is gaining, which illustrates how important groundwater levels are to conditions in the river.
As shown in Figure 15, the river is dominantly losing in the area where a number of wells are concentrated. To assess whether the pumping conditions within the lower reaches of the river are the cause of losing conditions, the model was rerun with the groundwater extraction rates set to zero. Figure 15 shows the results of the simulation with and without pumping for 2008. While the overall shapes of the curves are consistent, there are differences in the magnitudes of exchanges (highlighted within the ovals). With no pumping, the losing condition that is evident at 43000 and 44000 m during pumping becomes dominantly a gaining condition. Within the losing segments, the large negative peaks are lessened with no pumping, nearby, and at a fairly large distance (kms) from the wells. This result suggests that the pumping wells can lower the water table such that the effects are manifested at large distances.
3.5. Comparison of Simulated GW-SW Exchanges with Field Data
In-stream data throughout the Cowichan Region are limited due to data collection challenges including: the bedrock and gravel river substrate makes installing piezometers difficult; river discharge is high throughout much of the year rendering it unsafe to make in-stream measurements; and the perceived dangers of using of chemical tracers (e.g., solute and fluorescence tracers) on a Canadian Heritage River. However, some data were collected during the summer low flow season in 2013 at a few in-stream locations (S. Barroso, BC Ministry of Forests, Lands and Natural Resources Operations personal communication). The data include a series of in-stream mini-piezometer measurements of hydraulic head differences between the river stage and shallow groundwater levels within the river bed (using a pressure manometer board), as well as seepage rates between the shallow aquifer and the riverbed (using the same piezometer apparatus as a seepage meter). The seepage measurements (volumetric flow) and the modeled MIKE SHE exchange flow values were converted to a flux (m/s), by dividing the measured flow by the surface area. An additional source data came from snorkel surveys (fish count and habitat) that have been historically conducted within the Cowichan River (Mike MCCulloch, BC Ministry of Forests, Lands and Natural Resource Operations, personal communication). Fish count numbers, as well as descriptions of the habitat, were made. Indications of groundwater welling (gaining reaches of the river) often coincide with areas where fish counts are large and decrease in the temperature of the water.
Figure 16 shows the geographic positions of the gaining portions observed from snorkel surveys (blue markers), the locations of losing portions from seepage measurements (red markers), alongside the model results. The model is accurate overall in representing the gaining and losing conditions along the Cowichan River. Groundwater welling indicated by the snorkel surveys correlate well with the gaining conditions in the majority of the upper Cowichan River, although the gaining conditions from the model are not strong due to the low hydraulic conductivity of the sediments and bedrock within that portion of the river. Overall, the first 40,000 m (40 km) of the river is dominantly gaining (small magnitude), while the bottom 10,000 m (10 km) is losing (large magnitude).
3.6. Climate Change Simulation Results
The climate change results were analyzed over the last 10 years of the full 14-year simulation period to avoid the model spin-up time. The results represent a ten year time span during each of the 2050s and 2080s. Compared to the annual water balance values for the baseline model, the following trends are observed over time (baseline to 2050s to 2080s):(i)precipitation increases, with subsequent increases in runoff (overland flow) to the Cowichan River;(ii)evapotranspiration increases;(iii)all other aspects of the water balance remain fairly constant, including recharge, which is shown to increase only slightly; the estimated changes in recharge are within the uncertainty (error) range in the model.
Table 6 summarizes the changes to precipitation, ET, and recharge on a monthly basis for the baseline model and the 2080s (as amounts and percent changes). Key observations are as follows:(i)precipitation rates increase (relative to baseline) from September through to June, with the greatest increases in April, October, and November;(ii)ET rates increase (relative to baseline) throughout the entirety of the year, with the greatest increases from December to January;(iii)recharge rates increase (relative to baseline) for all months except June and August; the greatest increases (63%) occur in September.
The most noticeable effects of climate change within the Cowichan Watershed are related to snow. The continued increases in temperature consistently decrease the amount of snow accumulation (water storage) and alter the melt timing (earlier melt) as projected for other regions of BC and the Pacific Northwest [62–64]. Snow accumulation within the Cowichan is especially sensitive to climate change due to the dependancy of altitude for snow accumulation (currently simulated at above the 200 masl snow line). A warmer climate means that rain, as opposed to snow, will fall at progressively higher elevations during the winter months and elevations where snow accumulation is currently limited may have less winter snowpack and that snowpack will melt rapidly. Figure 17 illustrates the simulated spatial snowpack for the Cowichan region under the current climate condition, the 2050s, and the 2080s. A drastic decrease in snow accumulation is projected for the 2050s and 2080s. The snowpack becomes increasingly restricted to higher elevations, controlled largely by the temperature lapse rates, as temperatures within the valley are largely above 0°C. Both the spatial extent of the snowpack and the amount of accumulation within snowpack areas are greatly reduced. The timing of peak snowmelt also shifts from May-June to early January (results not shown).
As larger portions of winter precipitation fall as rain in future, the amount of water stored as snowpack decreases significantly, which greatly alters river flow dynamics throughout the year [51]. In general, in the Cowichan, the freshet will occur approximately 44 days earlier by the 2050s, and >100 days earlier by the 2080s. The simulated earlier freshet season results in increased peak flows during the winter months and lower flows during the summer and fall. Figure 18 shows the Cowichan River discharge (at the 08HA011 hydrometric station) near Duncan throughout the simulation for the baseline and climate change simulations. The higher resolution time series (bottom) shows that the peak flows in the winter increase by as much as 100 m3/s, while snowmelt-driven flows are no longer observed and summer flows are more than 50% less. These trends are fairly consistent for all model years. The hydrologic results are consistent with results of studies for other areas of BC [65, 66]. The simulation results suggest that the decreased summer flows may put additional stress on already sensitive aquatic habitat.
4. Conclusions
The MIKE SHE model was developed for the Cowichan Watershed with the intent to simulate the regional hydrology. Simplifications and assumptions were necessary to represent the unsaturated zone, the saturated zone, and the surface hydrology at a large scale. The following points summarize the key findings of the study:(i)the Cowichan River is dominantly gaining in the upper reaches except at a few isolated locations. At lower elevation, the river becomes dominantly losing;(ii)the aquifer hydraulic properties appear to be the main control on the magnitude of exchange that occurs, as most exchange occurs through the aquifers with the higher hydraulic conductivities;(iii)groundwater recharge over the extent of the watershed is spatially variable and ranged from approximately 253 to 630 mm/yr, with a mean of 438 mm/yr (17% of the annual precipitation);(iv)recharge varies significantly throughout the year. The highest recharge occurs in October and November (>100 mm/month), while a recharge deficit (P-ET) is indicated in the months of June, July, and August, largely reflecting precipitation patterns;(v)simulated groundwater discharge locations coincide with mapped springs and wetland areas;(vi)evapotranspiration ranges from 0.5 to 10 mm daily and is estimated at 1126 mm annually (44% of the annual precipitation);(vii)the water balance for year 2012 (extreme low flow conditions in the Cowichan River) shows significantly lower amounts of recharge and precipitation, with increased evapotranspiration, when compared to average conditions;(viii)groundwater pumping noticeably affects exchanges between the Cowichan River and the aquifer within the lower valley (near Duncan). Exchange conditions at this location change from gaining (no pumping included in the model) to losing (pumping included). Within the losing segments of the river, the large negative peaks in losses are lessened with no pumping;(ix)climate change is expected to influence the Cowichan Watershed in the following ways: precipitation and subsequent runoff increases; evapotranspiration increases; while all other aspects of the water balance remain fairly constant, including recharge, which is shown to increase only slightly;(x)climate change simulations show significant alteration to the accumulation of snow within alpine regions, as the snowpack in the 2080s simulation become increasingly limited to higher elevations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was funded primarily through a research grant to Diana Allen by the Cowichan Valley Regional District. The authors also acknowledge the feedback provided by Pat Lapcevic and Sylvia Barosso at the BC Ministry of Forests, Lands and Natural Resource Operations.