in Table Mountain Group aquifers of South Africa

Based on the field measurements of the physical properties of fractured rocks, the anisotropic properties of hydraulic conductivity (HC) of the fractured rock aquifer can be assessed and presented using a tensor approach called hydraulic conductivity tensor. Three types of HC values, namely point value, axial value and flow direction one, are derived for their possible applications. The HC values computed from the data measured on the weathered or disturbed zones of rock outcrops tend to give the upper limit values. To simulate realistic variations of the hydraulic property in a fractured rock aquifer, two correction coefficients, i.e. the fracture roughness and combined stress conditions, are adapted to calibrate the tensor model application. The application results in the Table Mountain Group (TMG) aquifers show that the relationship between the HC value and fracture burial depths follows an exponential form with the power hyperbola.


Introduction
Darcy's law is always used to estimate the groundwater flow in both porous and fractured media, depending upon realistic estimates of aquifer hydraulic conductivities (viz.k) and hydraulic gradients (viz.J) at the scale of problem.In the case of fractured rock aquifers, presentation and determination of the hydraulic conductivity prove to be challenges.With respect to a fracture set with a mean aperture of b and a parallel face distance of d, the following classic expression is adopted for flow through the set of conduits (Talobre, 1957;Jaeger, 1972): (1) where: ρ is the density of fluid μ is dynamic viscosity of fluid, which is 10 -6 m 2 /s for water at 20ºC g the acceleration of gravity J hydraulic gradient Eq. ( 1) represents an idealised type of flow behaviour that has been intensively studied, both experimentally and numerically, by many researchers.The term gb 3 /12μd in Eq. ( 1) is usually referred to as the hydraulic conductivity (HC) K for the set of fractures involved.For the determination of K value, many theories and methods have been developed.A series of results for one of the intrinsic properties of fractured rock aquifers have been obtained to various extents for more than 30 years.As a summary, there are thus far three approaches to the estimation of hydraulic conductivity of fractured rock aquifers, namely: • HC tensor approach based on statistic or stochastic methods of in situ fracture geometry and physical measurements • Fracture property field and laboratory tests for the parameter K evaluation • Inverse analysis on continuous or discontinuous problems dependent on numerical models and parametric calibrations.
The estimation of K values using either pumping or packer test is based on the assumption that the groundwater is flowing through a geological continuum.It is often an expensive exercise to estimate and predict the regional aquifer properties (viz.K and J, etc.) from local-scale hydraulic tests.Also the large variation of HC, both along borehole sections and in between holes, usually makes it difficult to determine the representivity of the parameters in terms of groundwater assessments.Even where a representative elementary volume (REV) can be defined, it may not be appropriate to directly apply the local test results to a regional aquifer.In porous media the REV can be very small, whereas in fractured media the REV may be very large or even does not exist in some cases (Kulatilake and Panda, 2000;Wang et al., 2002).The statistic methods for calculating HC tensor were developed in 1980s (Hsieh and Neuman, 1985;Hsieh et al., 1985;Oda, 1985;Tian, 1988).The results from these methods can successfully indicate 3-D principle HC values and directions by means of coordinate rotation of the incorporation of input data that derived from the surface measurements.The basic assumptions of the tensor approach are: • Groundwater flow is exclusively governed by fractures • The fractures through a rock matrix are well-connected • Flows between fracture sets do not interfere, or no deflection flow occurs.
For the ideal flow pattern with M sets of fractures involved in a study area, the hydraulic conductivity of the fracture sets is expressed in the form of matrix which reflects a sort of flow superposition: (2) where: K is the HC tensor matrix accounting for the anisotropic nature of studied media I is the unit matrix n is the direction cosine vector whose components are expressed in terms of the dip azimuth β and dip angle α of the fracture sets in the coordinate system where the x, y axes are pointing to the north and east direction, respectively, while z axis is pointing upward.
The elements of the matrix are dependent on the geometric and physical parameters of fractures, which were studied by many others (Tian, 1989).From Eq. ( 2) two key HC parameters, namely principal K values and orientations and their corresponding composite K value can be achieved by employing the techniques of linear algebra arithmetic.
However, the complexity of K in fractured rocks is far beyond what the existing models could handle since there are many geometrical and mechanical factors that impact on the flow through fracture gaps (Snow, 1969;Peter, 2001).Among them, fracture aperture is known to be most critical in controlling the quantity of flow, and the study of the aperture doubtlessly becomes a main focus in this paper.The factors affecting the aperture of fractures include geometrical and mechanical properties of the fracture walls such as rigidity, roughness, mineral fillings, and stress levels surrounding the fractures, in which the roughness has a crucial impact on the aperture around the surface zone (Lomize, 1951;Louis, 1974;Patir and Cheng, 1978).Therefore, an expression of equivalent aperture due to roughness is suggested to calibrate the original HC tensor model.Furthermore, by taking into account crustal stress, lithostatic and hydraulic pressures that act on fractures, together with the equivalent aperture, an expression of hydraulic aperture is accordingly developed for model calibration purposes.

Background of the Table Mountain Group aquifer
The Table Mountain Group (TMG) comprises a sequence of sedimentary units that extends from Vanrhynsdorp in the northwest to the Cape Peninsula in the south and then incurves eastward to Port Elisabeth (Fig. 1).As a part of African Craton, the distribution of TMG extends to an area of about 248 400 km 2 with the outcrops of 37 000 km 2 .It forms the backbone of the Cape Fold Belt that was produced in Permo-Triassic period extending from Australia through Antarctica and South Africa to South America (McCathy and Rubidge, 2005).The structural frame was modified as an effect of the break-up of Gondwanaland during the Jurassic to Cretaceous time periods, which led to a series of tensile and dextral displacements.
Due to the experiences of groundwater usage by borehole abstractions, together with the cognitions of lithological characteristics, stratigraphic build-ups and structural fabrics of the Table Mountain Group, it has been concluded that the TMG is a regional aquifer system which may extend to great depths (2 000~5 000 m b.g.l).It has also been recognised that the fractured rocks mainly consist of sandstones, siltstones, and sandwiched shales and mudstones that were formed during the Ordovician to the Silurian period, 500 to 400 Ma ago, with sedimentation along a south-eastward trough (Rust, 1967(Rust, , 1973;;De Beer, 2001).
As underlain by Precambrian metamorphic rocks, overlain by mid-to neopalaeozoic basin deposits and bounded by some regional faults such as Kango and Worcester faults, the southernmost aquifer system on the African continent has the potential to become a major source of bulk water supply for both agricultural and urban requirements in the Western and Eastern Cape Provinces of South Africa.Extensive exploration and exploitation of the groundwater resource in the aquifer system have been done for about 30 years; and more than 45 Mm 3 of groundwater is annually abstracted in about 30 locations for the requirements of municipalities, irrigation farmers and holiday resorts.Minor users are the smaller scale farmers, homesteads and stock farms.Currently, major problems faced are the lack of information on the properties of the huge fractured rock aquifers and shallow and deep groundwater circulation.With regard to the determination of key aquifer parameters, such as hydraulic conductivity, transmissivity and storativity, the estimates at different scales are largely interpreted from borehole tests.In terms of the hydraulic conductivity, a wide range of K values from 1.99 to 1.99×10 -3 m/d has been given in various locations of the TMG area (Rosewarne, 2001).In some well-fields the overestimation of the aquifer parameters including the K value leads to unrealistic recommendations on water supply capacity, causing continuous decreases of borehole water levels or even a water scheme failure (Jolly, 2001).

Adaptation of hydraulic conductivity tensor theory
It is well established that the direction of groundwater velocity (V) and hydraulic gradient (J) is usually not coincident with each other over time in fractured rock media.There is an included angle θ between V and J, and the J component on flow direction is: (3)  In an anisotropic medium, the HC is not a scalar any more compared with the isotropic one.The flow velocity is correspondingly expressed as follows: (4) where: i, j ,k are the unit vectors on coordinates x, y, z respectively.
Its parameter term is generally presented in the form: (5) Eq. ( 5) is the hydraulic tensor matrix, in which K xy =K yx , K yz =K zy , and K xz =K zx .Note that in Eq.( 2) the expression of the row matrix n is: The above K is a symmetric square matrix with three different eigenvalues and corresponding orthogonal eigenvectors satisfying the following relation: (7) where: δ ij is Kronecker's symbol λ i the eigenvalues U i is the corresponding eigenvectors associated to λ i Eq. ( 7) is the representative of a homogeneous equation group and the solutions of λ i and U i (i=1, 2, 3) may be obtained from Eqs. ( 8) and ( 9) respectively. (8) Using U i (i=1,2,3) as the basic vectors, the HC tensor reduces to the diagonal matrix: (10) Hydrogeologically it is easy to understand that the λ i (i=1,2,3) represent three principal HC values, namely K p1 = λ 1 , K p2 = λ 2 and K p3 = λ 3 , and the direction of K pi is given by Ui.Whereby the composite HC value of a given sample site can be represented by the geometric mean of the three principal HC: Furthermore, the relationship between HC components K x , K y and K z along x, y, and z axes defined above and principal HC tensor can be established via U i : (12) Assuming that the flow direction is known (Fig. 2), using Eqs.
(3) and ( 4) we have the following expression for HC along flow direction: (13) Alternatively, using the projective relations of V x /V=cosβ•sinα, V y /V=sinβ•sinα and V z /V=cosα (see Fig. 2), one may also write that: (14) The above analyses indicate that the three types of K values, namely composite HC K comp , axial HC K i (i=x,y,z) and flow direction HC K V , are physically different.However, they are related through the HC tensors and vector projections.The former (HC K comp ) may be accounting for the K value at a measuring site in the form of scalar that is averaged from three anisotropic principal HC.The HC K i (i=x,y,z) are the projections of three principal HC along original or user-defined coordinate system x,y,z.This is important in practice, for it is more convenient to use the K i (i=x,y,z) than to use the three principal HC directly in groundwater modelling processing.Note that the quantity of HC K V is not simply the quadratic sum of K i (i=x,y,z) as it used to be because the flow direction HC is more physically determined by the components of flow velocity and head gradient than geometrically determined by axial HC.There exists a non-linear relationship between K V and K i (i=x,y,z).Accordingly, if there exist not less than two sets of foregone values of K V and K i (i=x,y,z) obtained from hydraulic tests and surface measurements respectively, it is possible to deduce the flow direction by using the following linear equation: (15) where: .
The analysis on the thin sections sampled from TMG sandstones reveals that the porosity of intact rocks is very low or even null.This implies that the motion of groundwater in those aquifers is basically controlled by various types of discontinuities in the form of bedding planes, joints, faults, unconformities and weathered fractures.There are at least three sets of fractures in TMG rocks.Those are bedding fractures, conjugate joints, and the others structural and weathered fractures.
The above-mentioned relations have been programmed in MS Excel Workbook.In datum preprocessing, the basic data were measured at sites of the TMG outcrops from which the aperture (b) and distance (d) are geometric mean values and the angles arithmetic ones.Figure 3 shows the site locations and distributions of fracture orientations in the form of fracture density.The computed results, using the above-mentioned models, are listed in Table 1.The composite HC in Table1 are ranging from 1 to 21 m/d which are much higher than those of borehole tests mostly ranging from 5×10 -2 to 1.5×10 -3 m/d.This is because of the inevitable magnification of results when applying the measured data to the smooth plate model.Particularly, the apertures, measured at road cuttings and open quarries where the fractures that have long been undergoing disturbance and stress release, tend to be dilated compared with their original status.

Determination of hydraulic aperture
Equation (1) for groundwater flow in fractures is derived from the Navier-Stokes differential equation for pressureinduced laminar flow through the gap of two flat parallel plates where the hydraulic aperture is assumed to be uniform.The actual condition of the TMG is quite different, as the rocks including fractures have undergone phases of deformations or even distortions.What can be measured at rock outcrops is the mechanical aperture (Olsson and Barton, 2001) mostly ranging from 1 to 10 -3 mm.The hydraulic aperture discussed here is regarded the effective aperture for groundwater flow and can be obtained or inferred from both tracer tests (Charles, 1988) and laboratory experiments.Taking into account the main influencing factors, i.e. roughness and stress conditions, we rewrite Eq. ( 2) by adding the correction coefficients due to fracture roughness and stress condition respectively: (16) where: C er is the correction coefficient of roughness C es is the correction coefficient of stress condition.

375
From Eq. ( 16), it is easy to understand that the hydraulic aperture is the following: (17) Letting b er =C er 1/3 b, the b er may be regarded as the equivalent aperture due to fracture roughness.Correspondingly, letting b es = C es 1/3 b er , the b es is the equivalent aperture due to stresses acting on the rough fracture surface, namely the hydraulic aperture.Therefore, the following sections are focused on the discussion and determination of the C er and C es , respectively.

Fracture roughness
In the TMG rocks, fractures have rough walls and variable apertures and the internal properties of individual fractures that control groundwater flow are still poorly understood.The influence of roughness on fracture apertures has long been discussed and analysed by many others.Among them the important contributions are those of Lomize (1951), Louis (1974) and Barton et al. (1985).Different expression of the aperture correction coefficient due to roughness listed in Table 2 is a compilation of the results from them.A comparison of HC values obtained from the different models is plotted in Fig. 4, where the C er from Lomize and Louis' models shares the same form where the hydraulic radius is 2b in Louis' model.
Given the parameter measurements, it seems more practicable to apply Barton's JRC model for macroscopic fracture sets.That is: (18) where: b is the mechanical aperture basically measured by feeler gauge JRC is the joint roughness coefficient that can be measured with a profile-meter in a field or laboratory scale b er is the equivalent aperture calibrated by the roughness other than the hydraulic aperture in this study.
It used to be generally accepted that the groundwater flow is dominated more by the bedding fractures within the zone from surface to weakly weathered bedrocks than the structural joints because of the higher mechanical apertures within the bedding fractures.However, in most of the TMG rocks around the measuring sites, the bedding fractures tend to have a relatively higher JRC value in the range of 2~6.Comparatively, the JRC of the structural or conjugate joints is mainly in the range of 1~3.By applying equivalent aperture b er to the tensor model introduced, the difference of HC values between bedding and structural fractures is quite small and the results are listed in Table 3.Note that the HC value for the joint set at Site 5 is about one order of magnitude higher than that of bedding fractures.It is attributed to the relatively rougher surfaces in bedding fractures which produce more impedance to the inner flow on the surfaces than that of structural joints even though bedding fractures are much longer than structural joints.
The comprehensive impact of roughness due to bedding and structural fractures in the TMG sites on the HC value is listed in Table 4. Comparing with the computed results for the 5 sites in Table 1, the difference of the composite and axial HC values is obvious with three orders of magnitude of K value between smooth and rough fractures.And the HC results in Table 4 can also be roughly consistent with those from pumping tests in the TMG aquifers.
For shear tests of fractures, Barton et al. (1985) also develop an empirical formula to express the relationship between normal and shear displacements: where: d m is the mobilised dilatation angle and it can be determined as follows (Olsson and Barton, 2001): where: M is the damage coefficient JRC m the mobilised joint roughness coefficient JCS the compressive strength of the joint wall.
We assume that an arbitrary fracture surface, with a dip azimuth of β and dip angle of α, is subject to a combined three-dimension stress [σ x , σ y , σ z , τ xy , τ yz , τ zx ] T ; the normal stress σ N and shear stress τ s acting on the surface may be expressed as: In the present case, because rocks are normally subjected to compression by lithostatic pressure, it is convenient to denote the normal and shear stresses on the oblique fracture surface in the case of σ x =σ y . (23) Whence, the maximum and minimum values of σ N and maximum value of τ s may be obtained with respect to the derivatives of ∂σ N /∂α=0 and ∂τ S /∂α=0, these yield: (24) Fractures mainly control the deformation of a rock mass, and they may response to the change of normal or shear stresses even no failures happen; these can be formulated by the relations between the stress and strain increments: (25) where: Whence, the increments of strains may be expressed in terms of the increments of displacement: (26) where: b er is the equivalent aperture due to roughness u n the normal displacement u s shear displacement of the fracture.
Using Eqs. ( 25) and ( 26) yields: (27) The first equation in Eq. ( 27) may be integrated with respect to σ N and u n respectively.This yields: (28) In the case of shearing displacement in Eq. ( 27), the integral may be performed by using Eq. ( 23)~ ( 24) and assuming that d m ≈tand m : (29) We obtain: (30) Eqs ( 20) and ( 21) give expressions for the fracture displacement under normal stress and shear stress, respectively.Note that in these equations the normal stress term can be applicable to effective stress, i.e. σ e =σ N -r w h w .Considering the relationship between changed fracture aperture and the displacement under combined stress conditions, we have the hydraulic aperture b h : (31) (32)

Figure 5
Rock mass displacement and fracture closure under normal stress (after Goodman, 1974)

377
Using Eqs. ( 16) and ( 18), ultimate hydraulic aperture expression is: It is obvious from Eq. ( 33) that the hydraulic aperture is affected by both the physical and the mechanical properties of the fractures, particularly the roughness and the stress levels on fracture walls.The importance of this hypothetical model is that most of the fractured rocks in the TMG area with a thickness of 900~4 500 m are untouchable with current investigation tools.For this purpose the predictive models for hydraulic conductivity and flow behaviour are necessary and may become a key tool if validated.

Discussion and conclusion
Fundamental principles of HC tensor and some problems with their application to groundwater flow are discussed.It is observed that the HC value tends to be higher than that of hydraulic tests because the data sets measured on surface are often overestimated with respect to the fracture apertures even by means of scientific sampling and statistical methods.It is also pointed out that different types of HC values derived from HC tensor approach have different meanings physically, although they are related.
For the use of calibration coefficient of roughness C er , it is recommended that Lomize's model be used for microscopic fractures and Barton's model for relative macroscopic ones.
The derivation of negative exponential form of hydraulic aperture is based fundamentally on rock mechanics by considering both fracture roughness and stress conditions.Applying it to the tensor model, Eq. ( 2) can be rewritten as follows and another assumption needs to be added in order that the fracture orientations on the surface are consistent with those at depths: (34) It is possible to calculate the hydraulic conductivity over depths with the above hypothetical model on the site of the TMG, by assuming that the rock masses in an aquifer are subjected to lithostatic pressure.The vertical stress and horizontal stress acting on the rock masses are ρh and νρh respectively.Applying them to Eq. ( 22) to obtain the stress impact on the fracture surfaces, the hydraulic apertures are determined by Eq. ( 33).
The hydraulic aperture model that incorporates the roughness and combined stress condition is applied in three sites of Table Mountain Group area using the basic properties of rock materials as listed in Table 5.The computed results are plotted in Fig. 6 showing that the reductions of HC over depths do not simply follow the negative exponent law.They share a common expression in the form (Fig. 5-d): (35) where: K comp is the composite HC value h the depth of burial a, b and c the site-dependent constants.
In Fig. 6 more sensitive depth of HC inclination is around 1 000~1 200 m and the depth of fracture trending to be closed is 2 500~3 000 m where the HC value is less than 10 -13 m/d.This can be compared with the work of Snow (1969).
The results also show that the relationship between HC and depth of burial follow a type of exponential function where the power is a quadratic function of the depth h.This is attributed to the fact that projected stresses on fracture surfaces are determined by both the depth of burial and the orientations of the fractures and effaced by physical and mechanical properties of fractures.The model provides a scientific tool to simulate the property of TMG aquifers at depths.However, the validity of the model calibration coefficients needs to be further verified through more case studies and more detailed hydraulic tests.

Fig. 1
Fig.1 The TMG outcrop basically stretches along mountain ranges in Western and Eastern Cape of South Africa.

Figure 1
Figure 1The TMG outcrop basically stretches along mountain ranges in Western and Eastern Cape of South Africa

Figure 2
Figure 2Relationship between flow direction and hydraulic conductivity Figure 3 Showing geological settings of the TMG in the study area of structural syntax zone.The sizes of sampling sites are in the range of 20×30m ~ 60×60m.
Figure 4 K versus b er by various roughness models using data measured at the TMG sites Figure 6The variation of composite hydraulic conductivity over depths of burial incorporating with joint roughness.a) Site 1 -Robertson; b) Site 2 -Simon's Town; c) Site 3 -Bot River; d) Relation between lnK and h at Site 1.