Regional recharge assessment in the crystalline bedrock aquifer of the Kenogami Uplands, Canada

Abstract

Assessing recharge in unconfined fractured rock aquifers is a challenging task due to the discontinuous nature of this medium. In this study, a methodology is proposed based on an analytical interpretation of the regional hydraulic head profiles of an unconfined fractured rock aquifer, to determine its recharge. This calculation is conducted for a regional groundwater system based on a one-dimensional Dupuit-Forchheimer model in steady-state conditions that yields the value of the ratio of the aquifer recharge to the aquifer hydraulic conductivity. The approach is tested on a crystalline bedrock aquifer (Kenogami Uplands) of the Canadian Shield located in the Province of Quebec in Canada. The results show an average hydraulic conductivity of the aquifer of 4.3 × 10⁻⁷ m/s and a recharge equivalent to approximately 0.4% (3.5 mm/year) of precipitation. These results are in agreement with other studies conducted on similar crystalline aquifers in the Canadian Shield.

Editor D. Koutsoyiannis

Citation Chesnaux, R., 2013. Regional recharge assessment in the crystalline bedrock aquifer of the Kenogami Uplands, Canada. Hydrological Sciences Journal, 58 (2), 421–436.

Résumé

L'évaluation de la recharge des aquifères libres de socle fracturé est une tâche ardue compte tenu de la nature discontinue de ces milieux. Nous présentons dans cette étude une méthodologie permettant d'estimer la recharge d'un aquifère de socle fracturé libre à partir d'une interprétation analytique des profils de charge hydraulique observés dans l'aquifère. Ce calcul a été appliqué à un système régional d'écoulement d'eau souterraine considéré comme unidimensionnel et permanent lorsque les conditions de Dupuit-Forchheimer sont respectées. Il permet l'estimation de la valeur du rapport de la recharge de l'aquifère et de sa conductivité hydraulique. Cette approche a été testée sur un aquifère de socle cristallin du bouclier canadien (Hautes terres de Kénogami) situé dans la province du Québec au Canada. Les résultats ont permis la détermination d'une conductivité hydraulique moyenne pour cet aquifère de 4,3 × 10⁻⁷ m/s et d'une recharge équivalant à environ 0,4% (3,5 mm/an) des précipitations totales. Ce résultat concorde avec ceux d'autres études qui ont été menées dans des conditions similaires de socle cristallin du bouclier canadien.

Key words:

• regional hydrogeology
• crystalline rock aquifer
• specific capacity tests
• analytical solutions
• groundwater recharge

Mots clefs:

• hydrogéologie régionale
• aquifères de socle cristallin fracturé
• essai de débit spécifique
• solutions analytiques
• recharge des eaux souterraines

INTRODUCTION

Groundwater recharge is defined as the volume, or quantity, of water entering an aquifer at the water table. An aquifer can be naturally recharged by either precipitation or surface water infiltration from rivers and lakes, usually both. Recharge may be equivalent to net infiltration, assuming that the hydraulic conditions remain uniform when water is moving through the vadose zone of the aquifer to the water table. Quantifying this recharge is a very important objective in many applications of hydrogeology. This topic is especially relevant to the study of groundwater sustainability, which aims to preserve the dynamic equilibrium of the groundwater resources of an aquifer over time. Knowing the recharge allows us to estimate the quantity of groundwater that may be withdrawn with minimum effect on the baseflow provided by the aquifer. Recharge is also a critical hydrogeological parameter for calculating groundwater travel times and evaluating the vulnerability of aquifers. Excellent overviews on groundwater recharge processes and models are provided by Sanford (2002) as well as De Vries and Simmers (2002).

Recharge, however, is recognized as one of the most difficult components to quantify in groundwater systems (Bredehoft 2009, Stephens 2009). Recharge has been extensively investigated during the last few decades, principally in granular aquifers, and multiple tools are now available to quantify recharge. Scanlon et al. (2002) present a good review of these tools and the different techniques for assessing recharge.

A regional approach for the assessment of groundwater recharge consists of conducting a water budget analysis within an aquifer, considering the inflow and outflow parameters of water within a watershed system such as precipitation, runoff, infiltration and evapotranspiration. Another approach, known as the water-table fluctuation method, consists of monitoring piezometers installed in the vadose zone of an unconfined aquifer and investigating how the groundwater levels react to infiltration and drainage events (Healy and Cook 2002). By simulating infiltration with commercial codes that discretize the Richards equation (Richards et al. 1956), recharge can be simulated if we know the unsaturated hydraulic properties (volumetric water content and hydraulic conductivity curves as a function of pore water pressure) of the porous medium (Lu et al. 2011).

The special case of recharge in fractured rock aquifers remains a poorly understood phenomenon. As described by Heppner et al. (2007), recharge processes in fractured rock aquifers present additional complexities when compared to granular aquifers. In the case of fractured rock aquifers, individual infiltration events can cause large-scale rises of the water table over short time periods. Factors such as the effects of fractures on flow, storage of fast preferential flow, fracture–matrix interaction, film flow, fracture fillings and their swell/shrink nature, as well as fracture roughness and aperture make the investigation of recharge more complicated (Salve et al. 2008). A reliable technique for quantifying recharge in a context of fractured rock aquifers is the tracer technique (Heilweil et al. 2006, Gleeson et al. 2009b, Praamsma et al. 2009, Pronk et al. 2009, Levison and Novakowski 2012). Milloy (2007) used a modified water table fluctuation method to calculate recharge in a crystalline bedrock environment of the Canadian Shield. Physical techniques have also been investigated, not directly for measuring recharge, but for characterizing the properties of a fractured rock aquifer under unsaturated conditions (Liu and Bodvarsson 2001, Flint et al. 2006, Katsura et al. 2006). These results have been critical in furthering the understanding of and ability to quantify the recharge through the study of infiltration processes occurring in the fractures and the matrix in the vadose zone of unconfined fractured rock aquifers.

To this day, however, there remains a lack of research focused on the relationship between the fractures and the matrix under unsaturated flow conditions. As stated recently by Stephens (2009), “ … there is room for improving our understanding of recharge processes and for developing more accurate, simpler, and less expensive methods for quantifying recharge directly”. This phrase clearly summarizes the philosophy of the approach proposed here for quantifying recharge at a regional scale. Gleeson et al. (2009b) comment that recharge in fractured rock aquifers can be extremely localized, while Levison and Novakowski (2012) state that infiltration and consequently recharge in fractured rock aquifers can be extremely rapid. These findings also concur to justify the need for employing a regional approach to assess recharge, considering implicitly that tracer techniques may be appropriate for estimating local recharge, but not for averaging recharge in fractured rock aquifers on a larger scale.

In this study, an alternative method is proposed to assess mean recharge on a regional scale in an unconfined fractured rock aquifer from a one-dimensional Dupuit-Forchheimer model. The method is based on: (a) calculating the hydraulic conductivity of the aquifer; (b) calibrating the hydraulic head profiles (measured by piezometers) by theoretical equations along the principal direction of groundwater flow; and (c) calculating recharge from the parameters of the obtained equation and the hydraulic conductivity determined in the first step. This approach is applied to a sample region located in the Canadian Shield within the province of Quebec, Canada.

DESCRIPTION OF THE STUDY AREA AND ITS HYDROGEOLOGICAL FEATURES

The province of Quebec in Canada supports a growing population and industry that is dependent on groundwater. A number of the aquifers currently being investigated in Quebec as part of a groundwater mapping programme are constituted of fractured bedrock (e.g. the Saint-Lawrence Lowlands, the Appalachians, and the Precambrian Shield). This study will focus on a representative fractured rock aquifer located in the Saguenay Region where the Precambrian crystalline rocks are part of the Grenville Province within the Canadian Shield. The crystalline rocks of the Canadian Shield are mainly composed of three families of rock: anorthosite/granite/syenite-migmatite-gneiss. This bedrock aquifer is known to be fracture-controlled and is characterized by a low primary porosity and low hydraulic conductivity (Milloy 2007, Gleeson et al. 2009a, 2009b, Praasma et al. 2009, Sykes et al. 2009). The aquifer of interest, the Kenogami Uplands, is a fractured unconfined rock aquifer composed of anorthosite (Fig. 1), which is a Phanerozoic-aged structure in meridional Quebec (south of the 52nd parallel). The Kenogami Upland consists of a 600 km² rock plateau surrounded by the Saguenay River to the north and the Kenogami Lake to the south. Its eastern and western borders are marked by a sharp topographical transition with the plain. The surface of the plateau undulates at elevations ranging from 150 to 200 m a.s.l., whereas the surrounding lower plain varies between 50 and 150 m a.s.l. The Kenogami Uplands form their own watershed due to the plateau configuration of this region. Many outcrops are distributed throughout the area and much of the Kenogami Uplands is covered by a thin and minimal soil over the Precambrian crystalline anorthosite (Fig. 2). There is no overburden, explaining the unconfined conditions observed in the crystalline bedrock. The land cover consists mostly of boreal forest as well as small bushes; there is no significant agricultural or industrial activity, which ensures the preservation of the water quality of the aquifer. There are 326 private wells located in the unconfined fractured rock aquifer of the Kenogami Uplands. These give a depth to water table of between 0 and 50 m below ground surface depending on their position on the Uplands.

Fig. 2 Outcrop view of the fractured crystalline bedrock (anorthosite) of the Saguenay-Lac-St-Jean region located in Quebec, Canada.

Fig. 1 Study location: crystalline rock aquifer of the Kenogami Uplands located in the Saguenay-Lac-St-Jean region in the province of Quebec, Canada.

Climate is generally humid with short, hot and humid summers, cold and snowy winters, and rainy springs and autumns. Average monthly temperatures range from –16°C in January to 18°C in July. The Uplands capture a total mean annual precipitation of 930 mm/year (standard deviation 21 mm) (Fig. 3) recorded at station 7060400 (Fig. 1) at Bagotville (Environment Canada 2011). Note that the period spanned by the recorded data used to characterize the climate at this station is 1970–2011. The Kenogami Uplands constitute a zone of recharge for groundwater. However, from October to March, as most precipitation is in the form of snowfall, recharge is minimal. Most of the recharge occurs in the spring during snowmelt, as six months of accumulated snowfall suddenly become available for infiltration.

Fig. 3 Mean monthly precipitation, Station 7060400 in Bagotville (Environment Canada 2011).

There is a network of permanent lakes underlain by anorthosite. These lakes have a total surface area of 76 km² covering almost 13% of the Kenogami Uplands (see the hydrography network on Fig. 1). The Uplands are the recharge area of an unconfined, fractured anorthosite aquifer; as such the large number of lakes probably indicates that the permeability of the aquifer is low. There is an east/west groundwater divide almost through the centre of the Uplands (Fig. 1), from which groundwater discharges to the Saguenay River to the north, and to Kenogami Lake to the south. This large water divide and the two surface water boundaries cause the regional groundwater flow to be relatively bidirectional. The Kenogami Uplands therefore constitute an appropriate field for observing groundwater head profiles in a predominant groundwater flow direction within an unconfined fractured rock aquifer. It should also be noted that the groundwater discharges to the Saguenay River and to Kenogami Lake in the form of springs. Field observations by the author show evidence of these springs. A few small creeks flow over exposed bedrock at the topographic transition between the Uplands and the discharging zones (Fig. 1). The springs, forming the creeks, are fed by groundwater flowing out of the bedrock through water bearing fractures. The year-round low temperature (between 5° and 8°C) and high conductance (between 1000 and 1500 μS/cm) of the water measured in these creeks confirm a groundwater contribution.

Most of the residents of the Kenogami Uplands have their own individual wells and rely on groundwater supply. Data on the well characteristics (operating pumping rates, well design, static water level, specific capacity) were obtained from drillers' logs for 326 wells (Fig. 4). Most of the interpretations in our study are based on the observations made from individual wells. Indeed, the individual wells constitute the only available information on the fractured rock aquifer in the Kenogami Uplands. Note that a spatial geodatabase for the Saguenay-Lac-St-Jean region (Fig. 1) has been implemented within ArcGIS for the hydrogeological characterization of this region. This database provides the user with relevant information on aquifers and groundwater properties and is based on information from boreholes, among other sources of hydrogeological information (Chesnaux et al. 2011); it was used for this study. Note that the geodatabase offers a screening process allowing the removal of data suspected to be inaccurate or abnormal. This allows a good degree of confidence regarding the quality of the hydrogeological information on the different settings extracted for the purpose of this study.

Fig. 4 Location of the individual wells installed in the anorthosite bedrock in the Kenogami Uplands and of cross-section A–A′.

METHODOLOGY

The methodology used to calculate regional recharge in the unconfined fractured rock aquifer of the Kenogami Uplands is based on an analytical interpretation of hydraulic head profiles of the groundwater flow system in the aquifer. The static groundwater levels (water table position) were obtained from 326 wells; the hydraulic properties of the fractured rock were obtained from specific capacity tests. The methodology is explained in the following sections.

Analytical solution for calculating regional recharge based on a unidirectional regional profile of the water table

Dupuit-Forchheimer model of groundwater flow and recharge calculation

In the Dupuit-Forchheimer model (Dupuit 1863, Forchheimer 1886), the saturated groundwater flow occurs within an unconfined aquifer underlain by an impervious boundary. The aquifer is bounded by two fixed-head boundaries. The flow in the aquifer is considered horizontal and one-dimensional according to the Dupuit approximation, which ignores unsaturated flow and is essentially a 1D approach to solve vertical 2D plane problems in an unconfined aquifer (Strack 1989). The slope of the tangent to the water table is assumed to be small and thus the vertical component of the gradient vector is assumed to be negligible. The flow velocity is assumed constant. A uniform, annual average rate of recharge, W [L T⁻¹] is applied over the entire upper boundary of the aquifer, establishing a regional groundwater flow system. Flow is assumed to be steady state, which implies that the downstream discharge equals the groundwater recharge and that the water table position does not change over time. The conceptual flow system can be described using the Dupuit-Forchheimer differential equation, written as:

(1)

where K is the hydraulic conductivity of the aquifer [L T⁻¹] and h(x) the hydraulic head [L] along the x-axis. The general solution to Forchheimer's equation for the phreatic surface elevation (h) vs distance (x) is given by Bear (1972):

(2)

where L is the length of the aquifer [L], and h ₀ and h_L are the fixed upstream (x = 0) and downstream (x = L) heads, respectively (it is assumed here that h ₀ > h_L ).

After Forchheimer (1901), the phreatic surface h(x) in each case describes a semi-ellipse with a maximum elevation h _max at x _WD (groundwater divide position) either somewhere between the two boundaries or beyond the upgradient boundary. The gradient dh/dx is zero at h _max.

Equation (2)(2) can be used to determine the position of x _WD which can be expressed as (Bear 1972):

(3)

When the Dupuit-Forchheimer assumptions are followed, it is possible to calculate the ratio W/K from equation (2)(2) by plotting h ² vs x. Given a predominant 1D flow along the x-axis, if static water levels are monitored in observation wells distributed along the groundwater flow path x-axis, the saturated thicknesses h(x) of the unconfined aquifer can be calculated at different positions. The coefficient W/K can then be extracted from the h ² vs x curve after application of a quadratic regression. Assuming steady-state conditions of the groundwater flow system, W/K is constant; knowing the average homogenous and isotropic (K component in the x-axis direction) hydraulic conductivity K of the aquifer, it is then possible to calculate the average steady-state recharge of the aquifer. This method for assessing the recharge of an unconfined aquifer is therefore limited by a certain number of assumptions. However, in particular cases of regional groundwater flow systems where most of these assumptions and the Dupuit-Forchheimer assumptions are satisfied, the approach can ensure satisfactory estimates for aquifer recharge with a minimum of effort compared with other methods discussed earlier.

In the case of fractured rock aquifers, a critical aspect concerns the definition of the lower limit representing the datum of the aquifer from which hydraulic heads (or saturated thicknesses in the case of unconfined aquifers) are calculated. Indeed, the base of the aquifer, arbitrarily horizontal, will be the datum of the hydraulic heads; its position impacts the resolution of equation (2)(2). The datum for hydraulic heads must be calculated by defining the limit between the active zone of the aquifer and the non-active zone of the aquifer (aquiclude). The first zone is characterized with a higher hydraulic conductivity where most of the groundwater flow takes place, whereas the groundwater flow in the second zone is negligible compared to groundwater flow in the active zone. The separation between the two zones can be difficult to define because of the complexity of groundwater flow. One of the challenges underlined by Maréchal et al. (2010) lies in defining the effective thickness of the hydraulically active part of the unconfined fractured rock aquifer.

Dupuit-Forchheimer model applied to the Kenogami Uplands fractured rock aquifer

The conceptual Dupuit-Forchheimer model presented in Fig. 5 is applied to the Kenogami Uplands along the cross-section A–A′ shown in Fig. 4. This cross-section of the aquifer is appropriate, since it intercepts the water divide in the Kenogami Uplands sub-parallel to the two fixed-head boundaries, the Saguenay River and the Kenogami Lake. Groundwater flows along this cross-section towards the Saguenay River, to the north and towards Kenogami Lake, south of the water divide (Fig. 1), if the water table is considered to be recharge-controlled (Haitjema and Mitchell-Bruker 2005). It can also be anticipated that the length of cross-section A–A′ is sufficient (15.6 km) when compared to the thickness of the aquifer (active zone thickness) to allow for dominant horizontal flow in the direction of the cross-section to satisfy the Dupuit-Forchheimer assumptions. Both ideas will be discussed later when assessing the thickness of the active zone of the fractured rock aquifer.

Fig. 5 Conceptual Dupuit-Forchheimer regional flow model applied to the Kenogami Uplands along cross-section A–A′ shown in Fig. 4.

To summarize, the particular hydrogeological context observed in the Kenogami Uplands presents an opportunity to conduct the analytical quantification of recharge in a fractured rock aquifer of the Canadian Shield. Figure 4 shows that there is a good number of relatively evenly distributed wells, which allows interpolation of the hydraulic head profile along the cross-section before calculating the constant W/K. Also, the many lakes intercepted by the cross-section will be used to restrain the water table position, since the unconfined aquifer is in direct hydraulic contact with surface waters. Although groundwater levels fluctuate over time, it will be assumed that on a regional scale these fluctuations are negligible.

Characterization of regional hydraulic conductivity from specific capacity tests

Specific capacity tests have been conducted in 326 wells that intercept fractured anorthosite. The specific capacity test consists of pumping the well a short time after the well's completion to determine the potential yield of the well. The specific capacity of a well observed at a certain pumping flow rate defines the ratio of the observed yield and the corresponding final quasi-stabilized drawdown.

Data obtained from these tests have been extracted from the Quebec Ministry of Environment drillers' database. The driller's report usually contains the coordinates of the well location, well diameter, depth to rock, final depth to static water level, pumping rate (yield), test duration and dynamic water level at the end of the test.

Assessing hydraulic conductivities of aquifers is challenging, especially in the case of discontinuous medium like fractured rock aquifers and due to the scale effect (Guéguen et al. 1996). Studies of scale effects on hydraulic properties of fractured rock aquifers (Nastev et al. 2004, Lavigne et al. 2005) have concluded that the specific capacity test, because of the important amount of data, gives a sufficient representation of the average hydraulic properties of the aquifer.

Transmissivity and hydraulic conductivity values can be obtained from the specific capacity test by using the Cooper-Jacob approximation (Cooper and Jacob 1946) of the Theis (1935) equation during the transient phase of the pumping test:

(4)

with T the transmissivity of the aquifer [L² T⁻¹], Q the pumping rate [L³ T⁻¹], t the time [T], s(t) the drawdown [L] measured at time t, S the storage coefficient of the aquifer [-] and r_w the radius of the well [L]. Note that the ratio Q/s defines the specific capacity of the well when quasi steady-state conditions are reached when pumping at a Q flow rate.

The Cooper-Jacob approximation expressed in equation (4)(4) theoretically only applies to confined aquifers. However, in the case of unconfined aquifers, the drawdown in equation (4)(4) can be corrected. This is the case in our study, which considers that the crystalline bedrock of the Kenogami Uplands is unconfined. The corrected drawdown s_c must be defined according to Jacob (1947):

(5)

where h_i [L] is the initial saturated thickness of the unconfined aquifer. Note that knowing h_i supposes that we also know the static water level and the thickness of the aquifer represented by the active zone of the groundwater flow in the crystalline bedrock, as mentioned previously (Maréchal et al. 2010). Using s_c given in equation (5)(5), equation (4)(4) is commonly applied to unconfined aquifers. Such use is reasonable and does not introduce significant error when evaluating transmissivity (Razack and Huntley 1991, El-Naqa 1994). Equation (4)(4) is applicable when well losses and skin effects are negligible. Total drawdown in a well can be calculated using a step-drawdown test (Rorabaugh 1953). However, such tests were not performed on the wells of the Kenogami Uplands. Well losses and skin effects are important to consider in unconsolidated sediments where there is a screen and a gravel pack in the well that may produce turbulent flow (Jalludin and Razack 1991). Well losses and skin effects, if not taken into account, lead to underestimation of the transmissivity of the aquifer (Razack and Huntley 1991). However, in the present study, the 326 wells of the Kenogami Uplands are located in fractured rock, which implies that the wells consist of an open borehole. In such cases, well losses are more limited than in the case of unconsolidated sediments, and thus it is reasonable to neglect well losses in these environments. Further discussion on this assumption is provided in the sensitivity analysis section of this article.

Another source of uncertainty when using equation (4)(4) to obtain the transmissivity is the wellbore storage effect, which affects drawdown behaviour at the early stage of pumping. Indeed, during this phase, the Cooper-Jacob assumptions are not respected and the calculation of the transmissivity would not be accurate (overestimation of T due to storage of water in the borehole) when using a drawdown value at too early a stage of pumping. The pumping time required to neglect the influence of the wellbore effects depends on the pumping rate, but usually for specific capacity tests conducted in private wells, the first minutes are critical. Consequently, it was verified that the 326 wells extracted from the geodatabase (Chesnaux et al. 2011) do have a pumping duration of at least 30 min before including them in the study (the pumping tests considered in this study lasted between 30 min and 2 h).

It must be noted, (a) that the transmissivity is contained in equation (4)(4) in both an arithmetic and logarithmic expression, and (b) that the equation contains two variables: K and S. Bradbury and Rothschild (1985) proposed a computerized technique to solve T in equation (4)(4) after establishing a value for S. They noticed that the result is not very sensitive to the S value, which can be arbitrarily chosen, as it is included in the Neperian logarithmic term (see equation (4)(4)). The value of S depends on the type of aquifer that is considered: for unconfined aquifers, S corresponds to the specific yield (or effective porosity) of the aquifer, whereas for confined aquifers, S corresponds to the elasticity properties of the aquifer. In the particular case of the Kenogami Uplands unconfined fractured rock aquifer, a reasonable value of S was established to be 0.01, since the specific yield tends to equal fracture porosity (Healy and Cook 2002, Caine and Tomusiak 2003). The porosity of crystalline bedrock in the Canadian Shield has been reported to be consistently of the order of 1% (Gleeson 2009, Sykes et al. 2009). Other studies have determined the same order of porosity (between 0.01% and 4%) in fresh crystalline rock (Belghoul 2007).

Equation (4)(4) can then be solved from the drillers' report data on the specific capacity test and S = 0.01; an iterative technique was set up in an Excel spreadsheet (instead of the original FORTRAN program defined by Bradbury and Rothschild 1985) and the value of T converges and becomes constant after only a few iterations. In this way, a transmissivity value was determined for each of the 326 wells located in the fractured bedrock of the Kenogami Uplands (density about 1 well per 2 km²). Hydraulic conductivity was then determined by dividing the calculated transmissivity by the saturated thickness b of the aquifer (K = T/b).

RESULTS

The results of the hydraulic conductivity are presented first, followed by the calculation of aquifer recharge based on the analytical interpolation of the piezometric observations in the 326 wells.

Hydraulic conductivity estimation

Results on hydraulic conductivity estimations

The 326 values of T calculated using the iterative method applied to equation (4)(4) range from 2.8 × 10⁻⁷ to 1.90 × 10⁻³ m²/s with a geometric mean T value of 2.3 × 10⁻⁵ m²/s with S = 0.01. Note that the initial saturated thickness of the unconfined aquifer h_i was calculated relative to a groundwater flow active zone considered to be of 100 m. The justification for this value is given in the next section. However, as mentioned previously, there is no significant impact on the calculation of T when considering the corrected drawdowns s_c (equation (5)(5)) instead of s. This is explained by the fact that the observed drawdowns s are relatively high compared with the initial saturated thickness, considering an active zone of 100 m for the aquifer.

The 326 values of hydraulic conductivity calculated by dividing transmissivity by the saturated length in the wells range from 3.5 × 10⁻⁹ to 1.0 × 10⁻⁴ m/s. The overall geometric mean K value is 4.3 × 10⁻⁷ m/s with a standard deviation of log(K) of 0.86, typical for fractured rock aquifers (Fig. 6). A regional map of the hydraulic conductivities is given in Fig. 7 to visualize the distribution of the hydraulic conductivity values.

Fig. 6 Log(K) data histogram with hydraulic conductivity calculated from specific capacity tests in the 326 wells of the Kenogami Uplands installed in the fractured bedrock (S = 0.01). Note that the K values are in m/s.

Fig. 7 Map of the lineament locations and the different ranges of hydraulic conductivities obtained from specific capacity tests in 326 wells located within the Kenogami Uplands and installed in the fractured bedrock.

Table 1 presents the results of a sensitivity analysis of the value of S on K and shows that the variation of K with S is negligible.

Table 1 Sensitivity analysis results of specific yield S on hydraulic conductivity K

S (-)	K min (m/s)	K max (m/s)	Std logK	K geometric mean (m/s)
0.01	3.5 × 10^−9	1.0 × 10^−4	0.86	4.3 × 10^−7
0.001	6.1 × 10^−9	1.2 × 10^−4	0.83	5.8 × 10^−7
0.0001	8.5 × 10^−9	1.5 × 10^−4	0.82	7.2 × 10^−7
0.00001	1.1 × 10^−8	1.7 × 10^−4	0.81	8.5 × 10^−7

The lineaments observed from the interpretation of aerial photographs have been included in Fig. 7. No correlation between proximity to lineaments and the presence of higher hydraulic conductivity values could be detected for the Kenogami Uplands. These observations are supported by recent findings of Gleeson and Novakowski (2009) in the Canadian Shield bedrock indicating that lineaments may act as a barrier to groundwater flow and, therefore, do not necessarily represent high conductivity areas in fractured rock.

Active groundwater flow zone of the crystalline rock aquifer

The hydraulically active zone of fractured rock aquifers is known to be located in both the shallow weathered layer and the underlying densely fractured zones of the crystalline basement. In order to identify the active zone, the hydraulic conductivity estimated in the wells was computed with the open borehole length of each well. Figure 8 shows that the hydraulic conductivity of the fractured rock decreases as the aquifer thickness increases (open borehole length), when the thickness is between 0 and approx. 100 m. In this range, hydraulic conductivity values vary from 10⁻⁴ to 10⁻⁷ m/s. However, at thickness greater than 80–100 m, the hydraulic conductivity appears to reach a plateau, represented by a low value of about 10⁻⁸ m/s. Several different authors (e.g. Lavigne et al. 2005, Maréchal et al. 2010) have documented a decrease of hydraulic conductivity with depth. These authors report that the shallower part of bedrock is more susceptible to chemical and physical weathering increasing the number of connected fractures and hence increasing the K value; at a thickness greater than about 100 m, hydraulic conductivity appears to plateau at a low value, lower than about 10⁻⁹ m/s, making the crystalline bedrock below this depth quasi-impermeable vis-à-vis the shallower bedrock above. According to these previous results in crystalline bedrock, which confirm the field observations at the study site (Fig. 8), it is assumed in this study that the first 100 m of crystalline bedrock represents the active zone of the fractured rock aquifer where regional groundwater flow mainly occurs. Below this depth, groundwater flow can be neglected in favour of the shallower flow. This statement is important when defining the hydraulic heads measured in the observation wells. In unconfined aquifers, the hydraulic heads are represented by the saturated thicknesses of the aquifer, which implies that the hydraulic head datum must be located at the base of the aquifer. Considering that below a depth of 100 m the groundwater flow in bedrock becomes negligible suggests that below this level there is an impermeable medium. Consequently, the datum for calculating hydraulic heads (or saturated thicknesses) in the aquifer is not constant and will be fixed at 100 m below the mean top of bedrock elevation in the Kenogami Uplands.

Fig. 8 Variation of hydraulic conductivity as a function of the length of the open borehole.

It has to be noted that Sykes et al. (2009) drew a similar conclusion when they noted that the hydraulic heads measured in the shallower crystalline bedrock of the Canadian Shield (<100 m) respond rapidly to the seasonal changes in groundwater recharge. They also mention that this rapid response corresponds to low storativity, high hydraulic conductivity, and good connectivity in the shallow, fractured rock networks. Some isotopic groundwater results in the Canadian Shield bedrock (Gleeson et al. 2009a) also suggest that young groundwater preferentially recharged by winter snow melt is found to a depth of at least 50 m. Finally, in another fractured rock aquifer in the St Lawrence Lowlands in Canada, and based on observations of K as a function of the open borehole length in 326 wells, Nastev et al. (2004) reported that an open borehole depth of 100 m may be considered as an average aquifer thickness.

Calculating regional recharge and conducting verifications

Interpolated hydraulic heads in the Kenogami Uplands aquifer

Regional piezometry has been defined using the 326 wells of the Kenogami Uplands, with the water levels of several lakes located within these Uplands. Piezometric levels are interpolated considering a quadratic relationship of the squared hydraulic heads and the spatial distance according to equation (2)(2). Figure 9(a) shows the interpolated water table position map (or piezometric map) of the Kenogami Uplands aquifer. Interpolation was realized using the ArcGIS software hosting the groundwater geodatabase of the Saguenay-Lac-St-Jean region (Chesnaux et al. 2011). The regional groundwater flow occurs as expected from each side of the east/west groundwater divide, discharging to the north in the Saguenay River boundary and to the south in the Kenogami Lake boundary. The mean regional hydraulic gradient i in both principal directions is equal to 9.25 × 10⁻³. A reasonable calibration between the hydraulic heads observed in the wells and the interpolated hydraulic heads is expressed by the root mean square (RMS), and yields a value of 5.3 m. The interpolated water table elevation profile observed along cross-section A–A′ is shown in Fig. 9(b). Note that the position of the groundwater divide is obtained at x = 8800 m. The RMS value of 5.3 m represents an acceptable uncertainty given that the water table elevations above the basement of the active zone ranged between 90 and 180 m along cross-section A–A′ presented in Fig. 9(b). Indeed, this value represents about 5% of the maximum difference between elevations.

Fig. 9 (a) Interpolated saturated thickness in the Kenogami unconfined fractured rock aquifer. (b) Interpolated water table elevation profile and DEM along cross-section A–A′. Note that the water table elevations in both pictures are expressed in metres above the basement of the active zone of the aquifer.

Calculation of saturated thickness profile and recharge

The conceptual model presented in Fig. 5 was applied along cross-section A–A′. The water table elevation profile obtained through interpolation and presented in Fig. 9(b) is transformed into a saturated thickness profile along A–A′ so as to respect the Dupuit-Forchheimer flow model conditions. According to Fig. 5, the Saguenay River discharge zone is located at x = 0, the water divide is located at x = x_WD and the Kenogami Lake discharge zone at x = L along the A–A′ x-axis. By plotting the squared saturated thickness (h ²) as a function of x, we obtain a second-degree polynomial relationship, since a quadratic interpolation was applied to generate the hydraulic heads in the piezometric map (Fig. 9(a)), satisfying the Dupuit-Forchheimer theory (equation (2)(2)). The saturated thicknesses are calculated assuming an active groundwater flow zone in the first 100 m of the crystalline bedrock (Fig. 9(b)). Note that a choice was made to consider a thickness of 100 m below the mean rock elevation. This choice fixes the datum at an elevation of –13 m a.s.l. A good fit is obtained after applying the quadratic regression to h ² vs x and the resulting equation is written as: h ² = –2.61 × 10⁻⁴ x ² + 4.78x + 9920 with h and x expressed in metres (Fig. 10). Therefore, the constant W/K in equation (2)(2) equals 2.61 × 10⁻⁴. The length of cross-section A–A′ is in the order of several thousands of metres (15 600 m). The mean K value of approximately 4.3 × 10⁻⁷ m/s can be considered (Fig. 6). In this case, the regional mean recharge in the Kenogami Uplands aquifer is calculated to be equal to 2.6 × 10⁻⁴ × 4.3 × 10⁻⁷= 1.1 × 10⁻¹⁰m/s, or 3.5 mm/year. This value corresponds to 0.38% of precipitation, the mean annual precipitation being 930 mm. Such a low rate of recharge is in agreement with previous studies reporting very limited recharge in crystalline bedrock corresponding to approximately 1% of precipitation (Rodhe and Bockgard 2006, Milloy 2007, Croteau et al. 2010).

Fig. 10 Interpolated squared saturated thickness (h ²) along cross-section A–A′ in the Kenogami Uplands unconfined fractured rock aquifer in the case of 100 m thickness of the active zone.

Sensitivity analysis and further characterizations

Figure 11 presents the interpolated saturated profiles obtained for different thicknesses of the active zone of the aquifer. It illustrates the regression obtained for the assessed 80-m thickness of the active zone in the case of the Kenogami Uplands, and it also presents other regressions for the active zone with a greater thickness. This allows us to perform a sensitivity analysis of the calculated recharge's dependence on the thickness of the active zone under consideration. Indeed, the choice of active zone thickness affects the saturated groundwater thickness in the aquifer and as such, the results obtained when applying the Dupuit-Forchheimer model. Table 2 compiles the results of the regression for active zone thicknesses ranging from 80 m to 200 m when considering the calculated geometric mean hydraulic conductivity, K = 4.3 ×10⁻⁷ m/s, (Fig. 6) and the mean annual precipitation of 930 mm.

Fig. 11 Interpolated squared saturated thickness (h ²) along cross-section A–A′. Different cases of the active zone thickness are presented in order to perform a sensitivity analysis.

Table 2 Sensitivity analysis of the considered active zone thickness on the calculation of recharge (mean precipitation is 930 mm/year)

Active zone thickness (m)	W/K	K (m/s)	Calculated W (mm/year)	% Recharge
80	2.3 × 10^−4	4.3 × 10^−7	3.1	0.34
100	2.6 × 0^−4	4.3 × 10^−7	3.5	0.38
120	3.0 × 10^−4	4.3 × 10^−7	4.1	0.44
140	3.3 × 10^−4	4.3 × 10^−7	4.5	0.48
160	3.7 × 10^−4	4.3 × 10^−7	5.0	0.54
180	4.0 × 10^−4	4.3 × 10^−7	5.4	0.58
200	4.4 × 10^−4	4.3 × 10^−7	6.0	0.64
Note: Data in italics represent the results for the actual estimated thickness of the active zone (100 m).

Table 2 shows that the recharge rate is actually not very sensitive to the active zone thickness, recharge remaining lower than 1% of precipitation in every case. The active zone thickness has been estimated at 100 m based on the observation of various authors in the Canadian Shield. In reality, we can modify this value and re-estimate the recharge rate if we decide to increase the active zone thickness. However, according to the results of the sensitivity analysis, the recharge rate would not change dramatically.

It must be noted that the smaller the thickness of the active zone, the higher the ratio of the aquifer length (L) to this thickness and the greater the validity of the Dupuit-Forchheimer model. In this condition, it becomes reasonable to ignore vertical flow in preference to horizontal flow. A rule of thumb for the Dupuit-Forchheimer approximation states that it is acceptable for L to be five times larger than the mean aquifer thickness (Kirkham 1967). In our case, we verified that the regional groundwater flow is predominantly horizontal according to this rule. This confirms that using a Dupuit-Forchheimer model is an acceptable assumption for the analytical determination of regional recharge.

The other parameter that has an influence on the calculation of recharge is the hydraulic conductivity. Indeed, when extracting the coefficient W/K which is a constant dependent on the mean thickness of the active flow zone as investigated in the previous section, W is calculated by multiplying this constant by K. Therefore, W is directly proportional to K which implies that the sensitivity on the calculation of W is the same as the sensitivity on the determination of K. In this study, K was calculated from the specific capacity tests realized on the 326 wells located in the Kenogami Uplands and a standard deviation of logK 0.86 was determined with a mean value of 4.3 × 10⁻⁷ m/s. This means that the average uncertainty on the determination of K is about of one order of magnitude and that consequently it is also of one order of magnitude on the calculation of the recharge W. This is a high sensitivity because it implies that recharge would be between 0.038% and 3.8% of precipitation, but these are still low rates of recharge. It is therefore important to relate the large uncertainty on K with the sensitivity of W to K. The mean annual recharge in this study is, however, still considered to be 0.38% according to the geometric mean value of K obtained from the 326 specific capacity tests conducted in the 326 wells. Recall too that the calculation of K was also sensitive to the value of the storage coefficient S (corresponding to the effective porosity of the aquifer chosen to be 1%) but not significantly. Finally, the calculation of K was sensitive to the well losses and wellbore storage which were neglected according to the previous discussion and justification. It should be noted that well losses and wellbore storage effects both tend to underestimate T and K, and consequently also to underestimate W.

In order to further characterize the groundwater flow system of the Kenogami Uplands, we can evaluate: (a) whether or not the aquifer is recharge-controlled (compared to topography-controlled), and (b) whether the calculated position of the groundwater divide x_WD along cross-section A–A′ calculated with equation (3)(3) agrees with the observed natural position of the surface water divide.

Haitjema and Mitchell-Bruker (2005) propose a condition that allows us to determine if the water of an unconfined aquifer is controlled by topography or recharge. They established a criterion based on the assumptions that the Dupuit-Forchheimer model is valid and that the water table is not perched. This criterion is expressed as follows, where the water table is respectively topography-controlled and recharge-controlled:

(6)

with η a unit-less constant, H [L] the average aquifer thickness and d [L] the maximum terrain rise defined as the largest distance between the average surface elevation and the terrain surface as measured at the highest terrain elevation. According to Haitjema and Mitchell-Bruker (2005), for one-dimensional flow between two parallel surface water boundaries, η = 8. This is the case in the presented model of the Kenogami Uplands aquifer.

Along cross-section A–A′ (Fig. 9(b)), H = 140 m, d = 57 m, and equation (6)(6) yields a value of 0.86 which, being less than 1, indicates that the water table of the Kenogami Uplands aquifer is recharge-controlled. In recharge-controlled water table systems, water tables are generally deeper than in topography-controlled water table systems and a shift is observed between groundwater divides and surface water divides.

By using equation (3)(3), the position of the groundwater divide can be calculated from the interpolated results of the saturated thickness along cross-section A–A′. The calculated groundwater divide position (Fig. 10) is at a distance x_WD  = 8800 m from the Saguenay River (or 6280 m away from Kenogami Lake along A–A′) and corresponds to a shift of about 600 m when compared with a surface water divide located at x = 8230 m. This observation confirms that the Kenogami Uplands water table is recharge-controlled and implies that considering the Dupuit-Forchheimer assumptions is correct and acceptable when establishing our regional simplified flow model for calculating mean annual recharge.

DISCUSSION AND CONCLUSION

A methodology for calculating regional recharge has been established for a regional aquifer satisfying Dupuit-Forchheimer conditions. This methodology may be applied to either granular aquifers or fractured rock aquifers, but is of particular relevance in the context of a fractured rock aquifer when the estimation of recharge is difficult due to the discontinuous nature of the media.

However, it is important to remember that recharge can greatly fluctuate over time, depending on patterns of precipitation and snowmelt. For example, in our study region, recharge is high in the spring due to melting of a significant snowpack, whereas recharge is very low during the winter due to limited infiltration. The recharge that was calculated using the described analytical approach assumes steady-state conditions and does not take into account seasonal fluctuations of recharge. This approach therefore provides a mean annual recharge value based on mean annual conditions yielding a mean water table profile of the regional aquifer. Therefore, in order to assess the seasonal fluctuations of recharge, other methods should be applied.

Furthermore, as underlined by Gleeson et al. (2009b), recharge in fractured rock aquifers can be extremely localized and rapid (Levison and Novakowski 2012), meaning that the regional approach would not be appropriate for a local investigation of recharge due to the spatial variability of the fracture networks and the variability of the fracture apertures. Local studies on recharge would involve a good knowledge of the structural characteristics, location and distribution of the main discontinuities (e.g. faults and fractures) controlling recharge at a local scale.

Even though this regional approach considers steady-state and homogenous conditions and, as a consequence, may be somewhat limited in terms of local accuracy, it remains a useful tool as it provides a mean estimation of regional recharge with minimum effort and investigation. Indeed, this method requires the investigator to simply gather hydrogeological information from regional wells. As illustrated in this study, if well data have already been compiled in a database, the required information is easily accessible and one needs only to estimate the hydraulic properties of the aquifer from hydraulic tests and to record the static groundwater levels observed in the wells. This study illustrates that, despite a limited budget and limited access for field studies, if a hydrogeological database is accessible, practitioners are easily able, under certain conditions, to assess regional recharge for multiple uses such as conducting a water resources budget analysis or estimating the vulnerability of an aquifer. The results obtained regarding the recharge and hydraulic properties of crystalline bedrock aquifers can later be used to construct a numerical model. The determination of these parameters offers interesting perspectives for groundwater flow modelling in crystalline rock aquifers.

We must also emphasize that this analytical methodology is only relevant under specified conditions and assumptions. The key parameters for calculating the regional recharge are the hydraulic conductivity, the groundwater static level measurements as well as the definition of the datum of the aquifer based on the definition of the active zone of the fractured rock aquifer. Regional recharge is sensitive to these parameters and they must be accurately assessed, as illustrated in this study, in the sensitivity analysis on recharge.

To conclude, the Kenogami Uplands aquifer is appropriate for a regional analytical estimation of recharge as it satisfies Dupuit-Forchheimer conditions. This may not be the case for all aquifers, and practitioners wishing to use the proposed methodology should be aware that appropriately similar conditions to those outlined in our example should be confirmed before applying this approach to other aquifers. It should be recalled that in order to meet the Dupuit-Forchheimer conditions, it is necessary to consider as a basic hypothesis that there is piezometric continuum of groundwater within the aquifer, even though in reality, we know that the structure of fractured rock aquifers is discontinuous. A continuous hydraulic connection between the fractures of the aquifer is implied in order to execute the calculations. This study illustrates how analytical studies may still be useful when other tools are not available. Regional estimations of recharge may be more accurate than local estimates, which require precise measurements as well as a good knowledge of the spatial variability of the hydraulic conductivity in order to obtain accurate values. Indeed, the variability of the local recharge in fractured rock aquifers is associated in part with the variability of the hydraulic conductivity which is characteristic of such types of aquifers. Modelling studies are an appropriate tool to perform recharge calculations in a fractured rock aquifer, but only when good knowledge of fractured networks and hydraulic conductivity fields is available. This is rarely the case.

Acknowledgements

The author would like to thank the “Ministère du Développement Durable et des Parcs (MDDEP)” of Quebec (Canada) and the different regional municipalities of the Saguenay-Lac-St-Jean region for their financial support of this study. The author would also like to thank Mélanie Lambert and Julien Walter for drafting of the figures and editing assistance. Professor Eduardo Kruse and another anonymous reviewer are thanked for their valuable comments which improved the quality of the manuscript. Grateful thanks are expressed to Ms Josée Kaufmann for editorial collaboration.