GIS-Based Hydrogeological-Parameter Modeling

INTRODUCTION

Hydraulic conductivity is one of the most important parameters of surficial deposits for hydrogeological study. At a watershed scale, the hydraulic conductivity changes vertically with decreasing downward from the earth surface and horizontally from lithological unit to lithological unit. The vertical variation of conductivity can be described by hydraulic conductivity decay factor that describes the exponential decay of hydraulic conductivity with the depth of the considered-point geologic profile. The horizontal variation can be described by hydraulic conductivity deviation that shows the difference between the hydraulic conductivity of individual deposit and the mean hydraulic conductivity in the watershed. A topographicsoil index model, used to estimate these two hydrogeological parameters for different lithology units using the Digital Elevation Model (DEM) and the water-well data from the study area, provides overall assessment of relative hydraulic conductivity of lithological units in a watershed without conducting expensive pumping tests. A regression model has been constructed to relate the depth of a water table obtained from each water well to the conductivity, decay factor, watershed area and slope derived from DEM using GIS technology. Results show that older glacial deposits such as Newmarket Till show the relatively lower conductivity whereas the newer deposits such as Oak Ridges Moraine (ORM) and Halton Till, late ponds and recent deposits, show relatively higher hydraulic conductivity. The former has been recognized as a block layer and the latter as excellent recharge layers for groundwater system. The results in this study not only lead to a general methodology for conducting similar tasks for the estimation of GIS-based hydrological parameters derived from the observed data, but also offer evidence to support the geological hypothesis proposed by geologists involved in the evaluation and assessment of groundwater resources in the Greater Toronto area (GTA).

GEOLOGICAL SETTING

The 30 000 km² study area of this research covers the Oak Ridges Moraine and its adjacent regions, including Lake Ontario in the south, Lake Simcoe in the north, Trent River in the east and Niagara Escarpment in the west. The corresponding elevation ranges from 234 meters (Lake Ontario in south) to 511 meters (Oak Ridges water divide in north). Southern Ontario, part of the Grenville Orogen, occurred 1 200 to 1 000 million years ago and ranged from Texas to Labrador (Johnson et al., 1992). Two structural layers are characteristic of this study area: a Precambrian crystalline basement layer that experienced long-term complex metamorphism and deformation and a Paleozoic layer that has remained virtually intact. The Oak Ridges Moraine area is located on the northeastern end of the Algonquin Arch, a small part of St. Lawrence lowland whose basement is complex structures such as folds and thrust with NE-SW orientations. Field works show that the dominant structures in the Paleozoic rocks are NE-SW fractures (Sanford, 1995). In addition, the seismic data that some of the vertical movement may take place during deposition of the overlying Quaternary and recent sediments beneath the floor of the lake. In southwestern Ontario, several thousand boreholes are available to extract sedimentary and deformational history during Paleozoic. Some of the fractures have been rejuvenated during or after Quaternary. Bedrock in the study area is an Ordovician or Silurian sedimentary rock underlain by Precambrian rocks and overlain by Quaternary glacier deposits. Precambrian rocks are not exposed in the Oak Ridges Moraine area. The bedrock is gener-ally composed of dense shale or limestone with low permeability. Average hydraulic conductivity is 10^-4 cm/s for regional aquifers, and 10^-7 cm/s for regional aquitards (Brennand et al., 1995).

The Quaternary stratigraphy in Southern Ontario was dominantly produced by the advance and retreat of the Laurentide Ice Sheet (Barnett, 1992). The period in which glaciation is active (both advance and retreat) is referred to stade. The period between two stades is referred to interstade when glacial activity is weak. Pre-late Wisconsinan deposits are scarcely exposed, while late Wisconsinan deposits widely cover the whole Ontario. The geological model is presented for a segment of the glacial deposits of the Oak Ridges Moraine area of southern Ontario (Sharpe et al., 1997). The model contains four units as well as incised channels dissecting these strata. The lowest element is the Paleozoic bedrock. The lower drift package overlying unconformably bedrock is composed of a number of units, including the Scarborough Formation and Don Beds. Overlying the lower drift is the Newmarket Till, a regional till sheet and aquitard. Through the till sheet are eroded the south-southwest-trending channels composed of silt, sand and gravel. The Oak Ridges Moraine truncates and infills the surface of the channels. Beneath the Oak Ridges Moraine, channels with sediment fill 150 meters thick may reach the bedrock. The channels eroded into the Newmarket Till especially beneath the Oak Ridges Moraine have been considered as excellent water aquifer for ground water assessment in the area. The main objective of this study is to provide statistical hydrogeological evidence to show the difference in terms of hydraulic conductivity. South of the Oak Ridges Moraine, a fine grained diamicton-lacustrine package forms the Halton drift.

TOPOGRAPHIC-SOIL INDEX MODEL

Topographic-soil index (I_ts) is a local drainage index derived for quasi-steady state conditions (Spivapalan et al., 1987). It is a synthetic index to represent the physical properties of topography and hydraulic properties of soil (in surficial deposits instead of in soils) at each point (pixel). The variation of hydraulic properties of soils is represented by the difference between local hydraulic conductivity (K_i) at each point and an average value (K_e) of hydraulic conductivity for a whole catchment. We determin the relationship between dynamic change of local water table relative to an average water table in the whole catchment and the variance of local hydraulic conductivity. The topographic feature, another parameter that influences variation of water table, is characterized by these two indices surface slope (β) and the area (a) drained through the local unit contour length at the point, as is represented below (Spivapalan et al., 1987)

where I_ts is topographic-soil index (free of a unit); a is area drained through the unit contour length at the pixel (m²); K_i is local soil hydraulic conductivity at the pixel (m/s); K_e is catchment average value of the saturated hydraulic conductivity (m/s); β is local slope angle (degree).

Spivapalan et al. (1987) proposed an equation for the local water table depth z_i, as shown below

where z_i is local water table depth at pixel; z is catchment average water table depth; λ is catchment average value of the topographic variable ln (λ/tan β); f is parameter that describes the exponential decay of hydraulic conductivity with the depth (z) of considered point at profile as in the following formula (Spivapalan et al., 1987)

where K_z is soil hydraulic conductivity with the depth of considered point at profile, and K₀ is soil hydraulic conductivity at surface. K₀ is greater than K_i in most cases due to loose compaction and more fractures.

Local water depth data from numerous wells in the area are available for our present study. The depth of water level in a well may be lower than water table (z_i) at the same location because the well efficiency may be less than 100%. The difference is called well loss. Theoretically, if the well diameter is very small and free of water pumping, the well efficiency can be supposed to be 100%. In this case, water depth in a well is equal to the water table depth. In other words, the well loss is zero. The well water depth is treated as water table depth in this study for the following two reasons. First, the majority of wells are measurement wells with a small diameter but without a water pumping when water depth was measured. Thus, the well efficiency is very high. Second, the difference of water depths from the mean value is used in statistics of the study (i.e., z_i in the equation 2); a constant well loss is canceled if it does exist. Therefore, the well depth data used in the equation 2 is reliable and reasonable. The average water table depth (z) in each catchment can be calculated with ARC/INFO. Therefore, the decay factor (f) and the deviation of hydraulic conductivity (d) of a lithology unit can be estimated using the following equation

where

Well depth z_i is a known variable and z is an average value of the well depths in a watershed. Both α and β can be derived from DEM with ARC/INFO software and c can be calculated. Therefore, the unknown variables f and d can be calculated using two or more well data for a specific lithologic unit. There are usually tens of wells in each unit, which can be used to get more precise estimations by regression. Given n wells from the same lithology unit, f and d values can be then estimated from the following equation

As in matrix format, the equation (5) can be rewritten as

where

The solution to the equation(7) can be obtained as shown below

Each watershed is covered by all or some of the following seven lithological unit bedrocks (very limited), lower beds, Newmarket Till, Oak Ridges complex, Halton Till late ponds and recent deposits. The area covered by each of the lithologic units varies from watershed to watershed. Four typical watersheds were selected to run the program and to calculate the deviation of conductivity (d) and decay factor (f) for each lithology unit, respectively.

The deviation of hydraulic conductivity (d) is the logarithmic deviation of conductivity away from its average, i.e., d= ln K_e-ln K_i, which is related but not equal to the hydraulic conductivity of lithology. Therefore, the value d can be positive or negative. The results obtained from the four selected watersheds consistently show strong correlation between lithology types and d values. The d value of the Halton Till and ORM are close to zero, which indicates that the hydraulic conductivity (K_i) of the Halton Till and ORM is close to the mean hydraulic conductivity of the whole watersheds K_i ≈ K_e. The d value of Newmarket Till is negative, which means that its hydraulic conductivity is lower than the mean hydraulic conductivity of the watersheds K_i < K_e. The d values of late ponds and recent deposits are positive, implying that their hydraulic conductivities are higher than the mean hydraulic conductivity of the watersheds, K_i>K_e. In short, the younger the glacial deposits in the study area, the greater the hydraulic conductivity. The result is reasonable because younger deposits tend to have the looser texture than older ones and also more important because it supports a popular hydrogeological hypothesis in the area (Sharpe et al., 1997) that the Newmarket Till serves as block layer of ground water system, whereas ORM and late deposits serve as good recharge layers.

The decay factors of the lithologic units vary from watershed to watershed but the estimated values of decay factors are close to zero, indicating that the vertical change of hydraulic conductivity on a large scale (basin scale) is less significant than that of its horizontal variation.

CONCLUSION

A regression model is proposed to relate the variation of water well depth with topographic properties (area and slope), the variation of hydraulic conductivity and vertical decay factor. The implementation of this model in GIS environment (ARC/TNFO) based on known water data and DEM is used to estimate the variation of hydraulic conductivity and decay factor of different lithology units in watershed context. The results show that the lower members of the glacial deposits such as Newmarket Till show the lower hydraulic conductivity, whereas the upper members such as ORM Halton, late ponds and recent deposits, show the higher hydraulic conductivity, which is clear evidence to support the geological hypothesis for groundwater assessment in the area that the lower formation Newmarket Till serves as a block layer and the upper layers such as ORM as recharge layers of groundwater systems. The variation of vertical decay factor of the hydraulic conductivity mentioned above is not significant for such a large-scale (basin-scale) study, compared with that of horizontal decay factor.