A modeling study of heterogeneity and surface water-groundwater interactions in the Thomas Brook catchment , Annapolis Valley ( Nova Scotia , Canada )

A modeling study of heterogeneity and surface water-groundwater interactions in the Thomas Brook catchment, Annapolis Valley (Nova Scotia, Canada) M. J. Gauthier, M. Camporese, C. Rivard, C. Paniconi, and M. Larocque INRS-Eau, Terre et Environnement, Québec, QC, Canada Geological Survey of Canada, Québec, QC, Canada Université du Québec à Montréal, Montréal, QC, Canada now at: Golder Inc., Montréal, QC, Canada now at: Università di Padova, Padova, Italy Received: 4 March 2009 – Accepted: 5 March 2009 – Published: 27 March 2009 Correspondence to: C. Rivard (crivard@nrcan.gc.ca) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
The hydrologic response of a small catchment was modelled using a 3-D coupled model.The main objective of this study was to investigate the level of complexity required to simulate hydraulic connections and interactions between surface water (springs, overland flow, and streamflow) and groundwater (within unconsolidated sediments and the bedrock aquifer).This work is relevant to integrated water resources management, and to any application where detailed knowledge of water budget components and their interactions is crucial.
Physically-based models that feature some form of coupling between the governing equations of surface flow (stream and/or overland) and subsurface flow (unsaturated zone and/or groundwater) make it possible to investigate surface-subsurface interactions and to quantify several parameters (such as overland flow, return flow, saturation, and recharge) that are difficult if not impossible to measure in the field.Simplified models, in some cases analytically solved, can be applied when appropriate assumptions can be made that constrain or reduce the geometry of the study domain, the degree of heterogeneity within the system, the predominant directions of flow, the boundary conditions of interest, and the number and type of parameters (e.g.Pohll et al., 1996;Barlow and Moench, 1998;Singh and Bhallamudi, 1998;Anderson, 2005;Hantush, 2005).For more general problems however, such as applications to natural catchments with complex geometries and distributed parameters, numerical models based on the fully three-dimensional Richards equation for variably saturated subsurface flow and on one-or two-dimensional approximations to the Saint-Venant equations (e.g.kinematic or diffusion wave) are more appropriate (e.g.VanderKwaak and Loague, 2001;Morita and Yen, 2002;Panday and Huyakorn, 2004;Kollet and Maxwell, 2006;Kampf and Burges, 2007;Jones et al., 2008).Applying such models to real catchments can be data-and compute-intensive and quite challenging, thus a study investigating complexity in hydrogeological representation and response is relevant from both hydrological and numerical perspectives.For the present study the CATHY (CATchment HYdrology) model was applied; this model was first presented in Bixio et al. (2000) and Putti and Paniconi (2004), and fully described in Camporese et al. (2009a).
The impact of heterogeneity on surface and subsurface processes, and to what degree, and how, this complexity is to be represented in a simulation model is a much-studied problem in hydrology (de Marsily et al., 2005).In the example of subsurface hydraulic conductivity (K), heterogeneity can be integrated using a homogeneous representation based on some hypothesized equivalent or effective K value, or it can be described in detail based on the observed spatial variability of the porous medium.The former approach raises issues of whether an equivalent homogeneous medium actually exists, under what conditions it may do so, and what techniques (e.g.averaging schemes) are appropriate for establishing an effective parameter value (e.g.Philip, 1980;Binley et al., 1989;Bachu and Cuthiell, 1990;Kim and Stricker, 1996).The latter approach relies on having sufficient data to characterize the heterogeneity.This data that can be obtained through well-established techniques such as pumping tests and still-evolving methods such as borehole geophysics, both of which were applied in the Annapolis Valley study (Morin et al., 2006;Rivard et al., 2009).
In practice, it is difficult to know how much data is "sufficient", and obtaining even a few sparsely distributed mea-surements is very costly.Sophisticated upscaling, homogenization, renormalization, and other techniques are useful in reconstructing geological complexity from available data, and numerous theoretical sensitivity studies on simplified or idealized systems have been conducted to investigate the importance of heterogeneity on flow and transport processes (e.g., Sykes et al., 1985;Kabala and Milly, 1990;Sánchez-Vila et al., 1995;Renard et al., 2000;Cushman et al., 2002).For the real catchment studied in this work, a more straightforward approach was pursued, starting with a very simple representation of the Thomas Brook catchment and then gradually increasing the level of geological complexity in the model, based on field measurements, maps, and other databases.The coupled, distributed model allowed detailed representation of parameters and close examination of processes and responses, providing, in particular, insights on surface-subsurface interactions.From a first, homogeneous representation of the catchment, on which a rudimentary calibration of the model was performed, eight additional scenarios were simulated, introducing successively more realistic geology, hydraulic conductivity, and, in a final scenario, also snow cover.These nine scenarios were used to assess the impacts of different levels of detail in the representation of heterogeneity on processes such as return flow, recharge, groundwater levels, surface saturation, and streamflow.
The Thomas Brook catchment (8 km 2 ) was selected for this study because it was reasonably well-characterized during a regional hydrogeological study of the Annapolis Valley, Nova Scotia (Rivard et al., 2007(Rivard et al., , 2009) ) and other studies and because it is considered representative of the Valley in its geology, topography, and land use.

Description of the study area
The Annapolis Valley in Nova Scotia is a major agricultural region of Canada.The Valley is 100 km long and is located between the North and South Mountains, along the Bay of Fundy.The Valley includes five watersheds: Annapolis, Cornwallis, Canard, Habitant, and Pereau (Fig. 1).Concomitant with population growth and changes in land use in the last decade, water has become a crucial resource for the continued development of the region.Limited surface water supplies and surface water contamination have made groundwater the primary source of drinking water.In this rural region, groundwater contamination by nitrates is also an issue (KCEDA, 2000;Rivard et al., 2009).
The Thomas Brook catchment is located in the Cornwallis watershed, in Kings County near the town of Berwick.It is located in a rural region with apple, corn, strawberry, and potato as the major crops.Agriculture and Agri-Food Canada (AAFC, 2002) has been studying the surface water of this catchment using six gauging stations.The outlet station collects flow data throughout the year, while the other five stations are only operational during the growing season.Elevations on the Thomas Brook catchment range from 30 to 220 m.Surface water drainage is mainly north to south, beginning at the top of the North Mountains and flowing south into the Cornwallis River, which ultimately discharges into the Bay of Fundy.Flow within the regional bedrock aquifer is also topography-driven to a large degree, with groundwater gradients directed from the North Mountains to the centre of the Valley (Fig. 2).Because of the steep slope of the North Mountain cuesta, a discharge zone with numerous springs is located at the foot of this slope.These springs are commonly used as a water supply source by the Valley residents.
Average annual precipitation and potential evapotranspiration are, respectively, 1211 mm/y (from Environment Canada, http://climate.weatheroffice.ec.gc.ca for the Kentville CDA station) and 689 mm/y (based on the Penman-Monteith equation; see Gauthier, 2009), with rainfall heaviest in spring and autumn.Groundwater recharge has been estimated to be approximately 315 mm/year based on hydrograph separation applied to 2005 data (Gauthier, 2009).

Geological context
The Thomas Brook catchment comprises the three main bedrock formations of the Annapolis Valley and is covered by surficial deposits resulting from glacial events that are typical of the region.

Bedrock geology
Formations of the Triassic Fundy Group (from the Mesozoic) overlay the Paleozoic rocks that form the basement of the Annapolis Valley.They comprise, in ascending order: the Wolfville, Blomidon, and North Mountain formations.The bedrock geology map is presented in Fig. 3. Figure 4 shows a reconstruction of the geology of the catchment for a simplified vertical cross section along a north-south transect.The Wolfville Formation, from the Late Triassic age, is composed of reddish thickly-bedded medium to coarse-grained sandstone, cross-stratified, with subordinate conglomerate, typically in lenticular beds, and it is sometimes interbedded with siltstone (Hamblin, 2004).The strata generally dip 5-10 • NW and thicknesses up to 833 m are reported (Hamblin, 2004).This formation represents the best aquifer unit of the Valley.From fieldwork and existing data analysis, its average hydraulic conductivity was estimated to be 5×10 −5 m/s (Gauthier, 2009).
The Blomidon Formation, also assigned a Late Triassic age, overlies the Wolfville Formation and underlies the North Mountain formation.It contains the same rock types as the Wolfville Formation, although in different proportions, with fine-grained beds reported to be predominant.Thicknesses up to 363 m are documented and the strata dip 5 • NW (Hamblin, 2004).This formation represents the second best aquifer unit of the Valley and its average K is 1×10 −5 m/s (Gauthier, 2009).
The North Mountain Formation is the youngest unit of the Valley.It corresponds to a series of tholeiitic basalt flows capping the North Mountains and overlaying the Blomidon Formation (Trescott, 1968).It is characterized by a cuesta with the steep slope facing the Valley (Trescott, 1968).Columnar jointing, vesicular flow tops, and abundant vertical fracturing define these basalts, providing potentially significant transmissivity (Rivard et al., 2009).However, this formation generally represents a poor aquifer due to lack of horizontal fracture connectivity.Thicknesses up to 427 m are recorded for this formation and the strata generally dip 3-5 • NW (Hamblin, 2004).Its average Kwithin the Thomas Brook catchment was found to be 1×10 −6 m/s (Gauthier, 2009).

Surficial deposits
The Quaternary geology of the Valley presents a high degree of diversity associated with the major depositional systems of glacial settings.There are four glacial phases known in chronological order: the Early Wisconsinan-Caledonian phase, the Late Wisconsinan glacial maximum Escuminac phase, the Scotian phase, and the Chignecto phase.For more information on the glacial history of this region see Stea (2004) and Rivard et al. (2007Rivard et al. ( , 2009)).The Quaternary sediments are mainly composed of tills, glaciofluvial sands, and glaciomarine and/or glaciolacustrine clays.Figure 5 presents the surficial units of the Thomas Brook catchment.At the surface, a sandy till (Tb and Tv) is present on two-thirds of the basin and glaciolacustrine deposits (Lb.v), which consist of silty to sandy mud, are observed in the lower part of the catchment.These units are considered semipermeable.More permeable glaciofluvial sands can be found underneath the glaciolacustrine deposits at the southern end of the basin.
These surficial deposits are generally thin on the Thomas Brook catchment (less than 10 m) and do not represent good aquifers.Their hydrogeological role is nonetheless important as they have an impact on bedrock aquifer recharge and can offer some protection against diffuse source contamination (Trépanier, 2008).Hydraulic conductivities for these sediments estimated with a Guelph permeameter range from 1×10 −7 m/s to 1×10 −5 m/s (Gauthier, 2009).

Description of the coupled model
The distributed, physically-based CATHY model (Camporese et al., 2009a) integrates land surface and subsurface flow processes.The three-dimensional Richards equation, solved by the finite element method, represents variably saturated flow in porous media, whereas a path-based onedimensional diffusion wave equation is used for hillslope (rivulet) and stream channel flow, with a different parameterization for these two elements of surface runoff.The distinction between grid cells belonging to the hillslope and stream network systems can be made according to three different threshold-based options, based on criteria such as upstream drainage area, local terrain slope, and land surface curvature.A pre-processor for the model analyzes digital terrain data to identify the drainage network, and sets up the surface discretization from which the three-dimensional subsurface grid is projected.Cell drainage directions can be identified by the simple D8 scheme (the method used for this study, whereby one of eight inflow/outflow directions is taken on each cell) or by more recent nondispersive and dispersive methods (Orlandini et al., 2003).The exchange fluxes, or coupling term, between the surface and subsurface are computed via a boundary condition switching procedure as the balance between atmospheric forcing (rainfall and potential evaporation) and the amount of water that can actually infiltrate or exfiltrate the soil.Switching between specified head and specified flux boundary conditions occurs at surface saturation (zero pressure head) in the case of rainfall and at the "air-dry" pressure head value in the case of evaporation.Nested time stepping in the numerical resolution scheme allows multiple steps to be taken for the surface routing component within each time step for the subsurface flow equation, which is in turn linearized using Newton-based iterative schemes.Two different sequential data assimilation schemes, nudging and the ensemble Kalman filter (Paniconi et al., 2003;Camporese et al., 2009b), allow model predictions to be updated with spatio-temporal observation data of surface and subsurface state variables.The data assimilation feature was not used for this study.
Outputs from the CATHY model include surface ponding depths, overland fluxes, subsurface pressure head and moisture content values, and groundwater velocities.Numerous other variables can be derived from these main outputs, and in this way the model allows spatio-temporal quantification of aquifer recharge, catchment saturation, streamflow, and other processes that arise from the interactions between surface water and groundwater.Outlet streamflow is computed at the outlet cell of a catchment and includes both surface and subsurface contributions.These two hydrograph components are not readily distinguishable, however, since subsurface flow may become overland flow before it reaches the channel, due to runoff generation from infiltration excess or saturation excess mechanisms, and overland runoff may in turn reinfiltrate the soil further downslope.The CATHY model, which handles these various runoff and infiltration processes, could also accommodate seepage face boundary conditions, for instance along the vertical face of a streambank segment, but this type of boundary condition was not used for the Thomas Brook application, since the brook is quite shallow and has negligible banks.Seepage of course still occurs whenever the water table intersects the land surface (saturation excess runoff generation).
Overland flow computed by the model at a simulation time t is the sum over all surface nodes (hillslope and channel) of the surface fluxes generated at that instant.Return flow is the part of overland flow that comes directly from the subsurface (groundwater that returns to the surface).Specifically, it occurs when the subsurface module in CATHY computes an outgoing flux normal to the land surface (i.e.exfiltration) at a surface node that is ponded or saturated (note that an exfiltration flux computed when the surface is below saturation is instead a contribution to evaporative demand).
The recharge calculation in CATHY considers the vertical component of Darcy velocities across the water table when these are negative (i.e.downward).At each simulation time step and for each surface node, the water table position is determined by examining, from the bottom of the catchment to the top, each node in the vertical profile.The water table will lie between the first two nodes that are found for which the pressure head transitions from a positive value to a negative value.Locating this switch from the bottom up ensures that the water table corresponding to the bedrock aquifer is found, and not some localized (perched) water table that forms at the edge of an infiltration front.Negative vertical velocities for each water table node are converted to length units (multiplying by the time step size), summed over all nodes that are recharging, and then summed over all time steps to produce total recharge (e.g. annual in the case of a 1-year simulation).

Model implementation
Topographic data used for the Thomas Brook catchment was obtained from Geomatic Canada's Canadian Digital Elevation Data based on the National Topographical Data Base at a 1:50 000 scale and with a resolution of 20 m (http: //www.ctis.nrcan.gc.ca/).These data were processed in a geographic information system to coarsen the digital elevation model (DEM) to a resolution of 60 m from which the surface and subsurface discretizations for the model implementation were derived.
A main focus of this study was on the representation and effects of geological heterogeneity.To do so, features of the CATHY model allowing heterogeneity to be represented by layer (vertically) and by zone (laterally) were used.The defined zones represent lateral heterogeneity, as deduced from the bedrock and surficial geology data and from the inclined nature of the bedrock formations.Each layer can have its own pattern of lateral heterogeneity, thus in order to represent some of the more complex geological scenarios it was necessary to define as many zones as there are surface cells.This process and the subsequent attribution of material properties (for instance K) to each zone were performed using the gOcad and ArcGIS software packages.The seven principal units present in the Thomas Brook catchment are illustrated in the conceptual model of Fig. 4. In geochronological ascending order these units are the: Wolfville Formation; Blomidon Formation; North Mountain Formation; till; glaciofluvial deposits (sands); glaciolacustrine deposits (mud); and colluviums.
A flat base was used for the bottom of the study area.At the lowest topographic point, corresponding to the outlet of the catchment, a total thickness (surficial deposits plus bedrock aquifer) of 50 m was assigned.With a total topographic relief of 190 m for the catchment, the flat base configuration resulted in a maximum thickness (at the northern end of the catchment) of 240 m.A total of 17 layers were used for the vertical discretization, with one or more layers for each major geological formation.Each layer, except the bottommost one, was aligned parallel to the surface and assigned a uniform thickness.The thinnest layers (0.1 m) were those closest to the surface, needed to accurately capture the interactions between surface water and groundwater, including In passing from the DEM-based discretization of the catchment surface (2234 cells in total) to the finite element discretization of the subsurface, each cell was divided into two triangles.The 4468 triangles (connected by 2390 cell corners or nodes) were then projected vertically to the 17 layers, with each pair of subtended triangles giving rise to 3 tetrahedra.The subsurface domain was thus discretized into 227 868 tetrahedral elements (4468×3×17) and 43 020 grid points or nodes (2390×(17+1)).A grid of this size required calculation times of several hours for 1-year simulations run on a laptop computer with a 1.90 GHz processor.
Boundary conditions assigned to the discretized domain representing the catchment were no-flow across its lateral boundaries and its base and atmospheric forcing, subject to boundary condition switching, over its surface.Figure 6 presents the atmospheric fluxes used as input to the model, as well as the measured streamflow at the outlet for year 2005 (the simulated results shown in this figure are discussed later).The input fluxes are the difference between daily precipitation and daily potential evaporation data.They are usually positive (i.e.water input, eventually producing recharge) from mid-October to mid-March and predominantly negative the rest of the year.Only one year, 2005, was retained because it was the only one containing a nearly complete streamflow time series on record.Streamflow varies from 0.06 to 0.42 m 3 s −1 .Figure 6 shows that a long recession period, when baseflow is the main contributor to streamflow, occurred between early June and early October.Initial conditions for the model were obtained by simulating the drainage of the catchment from full saturation (see Sect. 5.1).

Description of simulation scenarios
Nine scenarios were developed to investigate the influence of heterogeneity and other factors on the model's response for the Thomas Brook catchment.The response variables examined included streamflow (outlet discharge), soil-aquifer interactions (recharge to the bedrock aquifer, groundwater levels), and surface-subsurface interactions (saturation dynamics at the land surface, return flow).The scenarios, of increasing geological complexity for the first seven, introducing regional-scale parameter values in the eighth, and with the addition of snow cover for the last one, are summarized below and illustrated in Fig. 7.Note in Fig. 7 that the topography for all scenarios, though not depicted in these schematic figures, is taken into consideration and is as depicted in Fig. 4. Hydrogeological parameters were assigned based on fieldwork results and existing databases.
Scenario 1 corresponds to a homogeneous medium.A mean bedrock hydraulic conductivity of 1×10 −5 m/s, a porosity (n) of 15%, and a specific storage (S s ) of 5×10 −3 m −1 were used.In scenario 2, a homogeneous surficial sediment layer was added with a mean Kof 1×10 −6 m/s and a porosity of 20%.The specific storage was maintained at 5×10 −3 m −1 for both the bedrock and the surficial deposits.For scenario 3, the bedrock aquifer was vertically subdivided into three parts representing the Wolfville, Blomidon, and North Mountain formations.The calibrated K values were based on fieldwork performed on the catchment.Scenario 4 incorporates individual values of n for the bedrock formations, based on Freeze and Cherry (1979) for the North Mountain (basalt) formation and on estimated values (Rivard et al., 2009) for the other two formations.In scenario 5, a more realistic representation (i.e.variable thicknesses) of the surficial sediment units were introduced, and in scenario 6 a sandy unit was added beneath the glaciolacustrine deposits, in the southern portion of the catchment near the outlet.For scenarios 1 to 6, there is no inclination of the bedrock formations.
From scenario 7 onwards the reconstruction obtained using the gOcad software (Fig. 4) was represented.This is the most realistic geological model of the system, incorporating relative thicknesses, shapes, and inclinations of the bedrock formations and additional details in the unconsolidated sediments, including permeable colluviums at the foot of the North Mountains.Based on pump-test results from fieldwork, the specific storage values for the North Mountain, Blomidon, and Wolfville formations were reduced relative to the S s values for the surficial sediment units.For scenario 8, K estimates obtained from the more extensive database of the Annapolis Valley study were used instead of the few available measurements taken on the Thomas Brook catchment itself, i.e. smaller K for the North Mountain basalts, glaciolacustrine deposits, and colluvium layer (Rivard et al., 2009).The K values for all other formations, as well as the porosity and specific storage values, were unchanged.
Hydrol.Earth Syst.Sci., 13,[1583][1584][1585][1586][1587][1588][1589][1590][1591][1592][1593][1594][1595][1596]2009 www.hydrol-earth-syst-sci.net/13/1583/2009/The final, ninth scenario used the same hydraulic property values as scenario 8, and added snow accumulation and snowmelt processes, a new feature of the CATHY model introduced specifically for this study in consideration of the potential importance of snow cover on the temporal patterns of infiltration, recharge, and streamflow.A simple algorithm was used to treat snow dynamics, whereby precipitation occurring on days with a recorded average land surface temperature below 0 • C was accumulated as snow, and allowed to potentially infiltrate whenever the daily temperature rose above 0 • C (in scenarios 1-8 all precipitation was treated as rainfall).For the 2005 dataset used in this study, this treatment of temperature and precipitation data produced six snowmelt events over the winter.This is considered to be reasonably representative of the mild winters that are a characteristic feature of the Annapolis Valley micro-climate.The hydrogeological parameter values used for scenarios 8 and 9 are summarized in Table 1.

Model calibration and initial conditions
Model calibration was first performed for scenario 1, adjusting only unsaturated zone and surface routing parameters based on streamflow (see Table 2), the other parameters  The same values were as a whole, and characterize the "at-a-station" and b =0.5 [/] used for both overland "downstream" relationships of Leopold and Maddock (1953); and channel flow regimes.k s (A s , 1) and W (A s , 1) denote the Gauckler-Strickler k s (A s , 1)=0.1 [L 1/3 •T −1 ] coefficient and the water-surface width at a site draining W (A s , 1)=60 [L 0.22 •T 0.26 ] area A, for a flow discharge equal to unity.* dimensionless were kept fixed for all subsequent scenarios.This rudimentary calibration thus provided a common basis for all nine scenarios, where complexity was introduced in an incremental manner by varying the geological representation, the hydraulic conductivity, porosity, and specific storage parameter values, and the processing of precipitation data.
For each scenario, initial conditions were established by simulating drainage of the catchment starting from fully saturated conditions and with zero atmospheric forcing.When the streamflow generated from this pure drainage simulation roughly matched the Thomas Brook baseflow, the corresponding pressure head values for each grid node were read in as initial conditions for the 2005 simulation of that scenario.Baseflow was taken to be about 0.1 m 3 s −1 (see Fig. 6), and the match was obtained after 1-2 months of drainage, depending on the scenario.Regeneration of initial conditions was necessary because of the changes in domain configuration and parameter values introduced with each new scenario.The procedure used is analogous to capturing the steady state water balance of the catchment (represented by baseflow in this case), and is a quite standard method for obtaining adequate starting conditions for long-term transient simulations (e.g.Stephenson and Freeze, 1974).

Model response
The various scenarios show different performance according to the observed (streamflow) or estimated (recharge) variables.In general, model performance based on streamflow improved from scenario 1 to scenario 8, but for recharge, results are variable (see Table 3).Scenarios 4 to 6 produced similar values for hydraulic components, with the best results for mean streamflows compared to observations; nonetheless, recharge was much too high.Scenario 7, with the most realistic representation of the system thus far, nevertheless yielded water levels that were too low in the North Moun-  tain formation and an annual recharge that was too high, at 675 mm (Table 3), likely due to the available local data not being very representative.This was corrected in scenario 8 when regional values were incorporated.When considering recharge, model performance also improved significantly from scenario 8 to scenario 9, after snow accumulation was taken into account (see below).
With the adjusted hydraulic conductivities, scenario 9 produced 349 mm of recharge, close to the estimated value for 2005 (315 mm/y).The match in groundwater levels at the end of the year 2005 for this scenario is shown in Fig. 8.The wells located in the North Mountain formation are those in the high range of the graph, while the lower elevation points correspond to wells in the Blomidon and Wolfville formations in the central and southern parts of the catchment.There was a good match over the entire study area, with perhaps a slight tendency for the model to overestimate Hydrol.Earth Syst.Sci., 13,[1583][1584][1585][1586][1587][1588][1589][1590][1591][1592][1593][1594][1595][1596]2009 www.hydrol-earth-syst-sci.net/13/1583/2009/ groundwater levels in the Blomidon Formation.A coefficient of determination (R 2 ) of 0.99, a mean error of 0.33 m, a mean absolute error of 3.7 m, and a root mean square error of 5.5 m were obtained.It should be noted that the measured data presented in Fig. 8 correspond to measurements taken at different times in private wells with varying (and often unknown) depths.The calibrated model was considered satisfactory.Figure 9 presents measured and simulated groundwater levels in two monitoring wells (their location is shown in Fig. 2).Both wells are located in the Blomidon Formation, but well #1 is very close to the boundary with the Wolfville Formation. Wll #1 is used for domestic purposes, which explains occasional drops in water levels.The simulated groundwater levels for these two points compared reasonably well with the observation data, with discrepancies over time ranging from 0 to about 2 m, well within the accuracy of the DEM, i.e. 5 m.Observed and simulated groundwater levels responded in a similar way to the rainfall pulses and to the extended drier period in the summer.However, the model's response to this atmospheric forcing was much more pronounced.This is an indication that despite the geological heterogeneity introduced in scenarios 8 and 9, additional local heterogeneities in the aquifer not represented in the model probably act to dissipate the fluctuations in atmospheric forcing.These heterogeneities could be related for instance to the presence of an overlying less permeable strata or to the effect of spatially variable fracturing.
Figure 10 presents the simulated outlet discharge for scenario 9 compared to the measured streamflow for 2005.In general, the model succeeded in reproducing the individual events throughout the year as well as the overall response, with an average annual simulated flow of 0.18 m 3 s −1 , compared to 0.15 m 3 s −1 for the observed hydrograph.However, the highest flow rates are overestimated and summer base-  flows are underestimated.Mismatches could perhaps be reduced with further refinement of the model parameterization, including surface routing parameters for which data were not available for this study.These results are a significant improvement over the streamflow results from scenario 1, shown in Fig. 6 (see also Table 3, discussed below).

Effects of heterogeneity
Streamflow, overland flow, return flow, and annual recharge from the nine scenarios are summarized in Table 3.In this table, the average value (over the 1-year simulation) of the instantaneous overland fluxes is larger than the average streamflow value because these fluxes have not yet been propagated.
Overland propagation brings with it losses to evaporation and re-infiltration, as well as time lags due to the distribution of travel times generated by the surface routing equation.The ratios between return flow and overland flow reported in Table 3 are quite high (as much as 90% for scenario 1).This reflects the influence of topography, especially the steep drop along the cuesta that generates a lot of return flow that is quickly propagated.Indeed this is an area of the Thomas Brook catchment where numerous natural springs are found.As geological complexity increased, the ratio of return flow to overland flow tended to decrease.This was particularly true as lower permeability surface layers were introduced in scenario 2. These had the dual effect of reducing exfiltration of water already in the subsurface (thus reducing return flow) and reducing infiltration of rain and ponded water (thus increasing purely surface-generated overland flow).The lowest ratio between return flow and overland flow (less than 50%) occurs for scenarios 8 and 9, for which the hydraulic conductivity of three zones was further reduced.
Introducing a surficial sediment layer in scenario 2 with a lower hydraulic conductivity than the underlying bedrock unit reduced recharge by 21% compared to scenario 1. Similarly from scenario 7 to scenario 8, the annual recharge decreased significantly (45%) due to the reduced conductivities   introduced for the basalts, the till, and the glaciolacustrine sediments.Another (smaller) decrease occurred for scenario 9, because snow accumulation favours overland flow and evaporation compared to scenario 8 where all rain and snow can potentially infiltrate as it falls.The effect of snow accumulation can be seen in Fig. 11, where the monthly recharge for scenarios 8 and 9 is plotted for the simulation year 2005.
Recharge was lower for scenario 9 compared to scenario 8 during the winter months, higher in the spring when the snow melts entirely (but not sufficiently high to compensate for the lower recharge during the winter, as some of the snow becomes runoff), and identical to scenario 8 during the summer months.Recharge was high in January for both scenarios because the model was adjusting to the initial conditions: the starting water table was quite low as a result of the pure drainage simulation used to generate the initial conditions.These high recharge values for January are not representative of mean values for this region, as observed data over 30 years suggest that the main recharge periods are in the spring and fall.
Varying the porosity (for example between scenarios 3 and 4) had, in general, only a minimal effect on the response variables shown in Table 3, suggesting a much greater sensitivity of the model to hydraulic conductivity, expected.Increasing heterogeneity had a significant impact on numerical performance.CATHY adapts the step size within a given range, increasing it when there is rapid convergence in the iterative scheme used to linearize Richards equation and reducing it when convergence is slow.Should convergence fail for a given time step, the simulation steps back and attempts the integration anew with a smaller step size.A sevenfold increase in the number of subsurface time steps was observed between scenario 1 and scenario 8, implying that much smaller time steps were needed as heterogeneity increased.Likewise, the number of convergence failures doubled over this range of scenarios, from 128 for scenario 1 to 246 for scenario 8.

Catchment behavior for different response variables
In     The pattern of simulated groundwater levels for 14 June is shown in Fig. 15 for a longitudinal (north-south) transect of the catchment.The water table is shallow almost everywhere along this transect, but shows a significant dip around the cuesta.For comparison, observed water levels varied between 0 and 30 m below the surface, with a median of 6.3 m and with very shallow water tables downstream.Even wells in close proximity showed large differences in water level measurements, especially in the North Mountains, likely due to the fractured nature of the aquifer.These results point to some of the difficulties in dealing with heterogeneity both in the field and in modeling, with measurements of even seemingly straightforward variables such as well water levels showing large spatial variability over geologically complex terrain, and with the model probably requiring a much finer local grid discretization to accurately capture the effects of www.hydrol-earth-syst-sci.net/13/1583/2009/ Hydrol.Earth Syst.Sci., 13,[1583][1584][1585][1586][1587][1588][1589][1590][1591][1592][1593][1594][1595][1596]2009 the steep topography and the large variations in conductivity that characterize the area around the cuesta, where the North Mountains transition into the Blomidon formation.

Conclusions
A field and modeling study of a small catchment (8 km 2 ) in Nova Scotia, eastern Canada, was conducted to investigate groundwater-surface water interactions and other hydrological processes, with a particular focus on the effects of hetero-The coupled, distributed model used was suited to generate outputs of the spatial and temporal distributions of a number of response variables such as infiltration, recharge, groundwater levels, flow velocities, surface saturation, overland ponding, and channel discharges.These variables cannot be realistically monitored or measured in such detail at the catchment scale, so simulations play an important role in analyzing the factors that affect hydrological processes and interactions.
The simulation of seven scenarios of increasing geological complexity, an eighth scenario using regional-scale conductivity values, and a ninth scenario incorporating snow accumulation and melting, showed the importance of integrating certain characteristics to adequately represent the various hydrological processes in the catchment.For instance, they showed that the description of subsurface heterogeneity must include, as a minimum, hydraulic conductivity values for the three bedrock formations, because of their strong influence on water table position, and for the various surficial units given the dominant role of the surface cover in rainfall-runoff-recharge partitioning.Other parameters had a relatively minor impact on the catchment results on an annual basis (for example geological formation dip, porosity, and snow accumulation and melting).The simulated heads, aquifer recharge, and streamflow at the outlet for scenario 9 were comparable to observed or previously estimated values.Natural springs at the foot of the North Mountain cuesta, an important feature of the Thomas Brook catchment, were also adequately reproduced.Although the fully coupled model used here is not designed for management applications, results from this study of groundwater-surface water interactions provide useful information for water management and aquifer protection in this rural area.For instance, the insights provided on the spatial and temporal distribution of recharge can be used to limit potential contamination.
As in any modeling study of real aquifers and catchments, there are numerous potential improvements that can be suggested with regards to the representation of processes, parameters, and boundary conditions used in this study.Even with the most complex model, simplifications are always necessary, owing to model structure limitations, lack of data, etc.For example, a model structure that considers not just classical porous media flow but also preferential flow through fractures could be an important consideration, especially for the North Mountain basalts.Various theoretical approaches for doing so are available, but require additional data on fracture geometry, density, and connectivity.As another example, further simulations over multi-year periods could be used to verify and refine the parameter values used in this single-year study.The model could then also be used to study the long-term response of stream discharge to hydrologic input scenarios arising from climate changes and other humaninduced stresses.
This study confirms that a distributed, physically-based coupled model can be used to simulate groundwater and surface water flow on a small and moderately characterized watershed.A significant level of heterogeneity was required to achieve adequate simulation results, with representation of surficial deposits and realistic estimates of bedrock hydraulic conductivities being particularly important.This work adds to the relatively limited body of literature reporting field applications of fully coupled models.

Figure 6 : 27 Fig. 5 .
Figure 6: Daily net atmospheric forcing (rainfall -potential evaporation) and measured and simulated (scenario 1) outlet streamflow for 2005 for the Thomas Brook catchment

Figure 6 :Fig. 6 .
Figure 6: Daily net atmospheric forcing (rainfall -potential evaporation) and measured and simulated (scenario 1) outlet streamflow for 2005 for the Thomas Brook catchment

Figure 7 :
Figure 7: Schematic of seven scenarios of increasing geological complexity: 1) homogenous conditions; 2) addition sediments; 3) definition of the principal bedrock formations; 4) incorporation of heterogeneity in porosity; 5) variab deposits; 6) addition of a sandy unit in the area around the catchment outlet; 7) most realistic representation corre reconstruction model.

Fig. 7 .
Fig.7.Schematic of seven scenarios of increasing geological complexity: 1) homogenous conditions; 2) addition of a layer for the surficial sediments; 3) definition of the principal bedrock formations; 4) incorporation of heterogeneity in porosity; 5) variable thickness of the surficial deposits; 6) addition of a sandy unit in the area around the catchment outlet; 7) most realistic representation corresponding to the geological reconstruction model.

Figure 8 :Figure 9 :
Figure 8: Simulated (scenario 9) and measured groundwater levels for diff Thomas Brook catchment.The solid line represents a perfect match.

Fig. 8 .
Fig. 8. Simulated (scenario 9) and measured groundwater levels for different wells in the Thomas Brook catchment.The solid line represents a perfect match.

Figure 8 :
Figure 8: Simulated (scenario 9) and measured groundwater levels for different wells in the Thomas Brook catchment.The solid line represents a perfect match.

Fig. 12 .
Fig. 12. Distribution of simulated recharge (in m/s) for the Thomas Brook catchment on 29 July for scenario 9.The top chart shows the atmospheric input for year 2005, with the red pointer on 29 July.

Figure 13 :
Figure 13: Distribution of simulated surface saturation state for the Thomas Brook catchment during the drier summer period (August 28) for scenario 9.The top chart shows the atmospheric input for year 2005, with the red pointer on August 28.

Fig. 13 .
Fig. 13.Distribution of simulated surface saturation state for the Thomas Brook catchment during the drier summer period (28 August) for scenario 9.The top chart shows the atmospheric input for year 2005, with the red pointer on 28 August.
Figs. 12 to 15, examples of the response variables derived from the model simulation for scenario 9 are presented.The recharge distribution for 29 July is shown in Fig.12.The recharge rates were highest at the foot of the North Mountain, mainly due to the presence of permeable colluviums; they are also high in the center-west due to a local topographic high.Figures 13 and 14 illustrate the differences in surface saturation response during dry and wet periods.In late summer (28 August; Fig.13), after an extended period of low rainfall, only the areas in or adjacent to the channel network were saturated, caused by the saturation excess (or Dunne) mechanism, i.e. the water table reaching the surface.During the wetter period around 11 November (Fig.14), a much greater portion of the catchment surface was saturated, and infiltration excess runoff (or Horton saturation) was also evident, although the Dunne mechanism for runoff generation was still dominant.In both figures, saturated zones at the base of the North Mountain cuesta (just below the recharge area; see water table profile in Fig.15below) are apparent, corresponding to the natural springs located in this region.It should be noted that what the model is labelling as Horton saturation is a lumping together of the classical infiltration excess mechanism for producing saturation or ponding at the surface as well as additional causes that may arise from overland or shallow subsurface flow contributions from neighbouring cells that converge to and saturate a given surface cell.

Figure 14 :
Figure 14: Distribution of simulated surface saturation state for the Thomas Brook catchment during the wetter autumn period (November 11) for scenario 9.The top chart shows the atmospheric input for year 2005, with the red pointer on November 11.

Fig. 14 .
Fig. 14.Distribution of simulated surface saturation state for the Thomas Brook catchment during the wetter autumn period (11 November) for scenario 9.The top chart shows the atmospheric input for year 2005, with the red pointer on 11 November.

Figure 15 :
Figure 15: Distribution of simulated water table depth along a transect of the Thomas Brook catchment on June 14 for scenario 9.The model grid layers are shown as dotted lines.

Fig. 15 .
Fig. 15.Distribution of simulated water table depth along a transect of the Thomas Brook catchment on 14 June for scenario 9.The model grid layers are shown as dotted lines.

Table 1 .
Hydrogeological properties of the various units for scenarios 8 and 9.
e. homogeneous representation of the catchment geology) an adequate agreement between simulation and observation based only on matching the general trend of streamflow at the catchment outlet for calendar year 2005 (see Fig.6).Once obtained, the parameter values in Table2

Table 2 .
Model parameters common to all scenarios.

Table 3 .
Summary of simulated streamflow, recharge, and two components of surface flow for the nine scenarios.Errors* relative to the observed streamflow and estimated recharge are given in parentheses.×100/Average measured value.** The annual mean has been calculated using the total volume over the 2005 calendar year divided by 365 days.