Assessment of Aquifer Vulnerability Using Integrated Geophysical Approach in Weathered Terrains of South China

3.1 ERT Method

The 2D ERT method can assess areas of the heterogeneous settings for the groundwater evaluation where applications of the other geophysical techniques are not appropriate [35]. It provides more depth penetration by increasing the electrode interval, and gives a 2D subsurface model with vertical and lateral changes in resistivity values. An inversion of the apparent resistivity in the ERT survey gives the systematic measurements that can enhance quality of the subsurface geological model [36, 37, 38, 39, 40]. ERT was conducted using WDJD-4 multi-function electrical instrument produced by the Chongqing Pentium CNC Technology Research Center and WDZJ-120 multi-electrode converter. A layout of pole-dipole configuration was used to acquire the resistivity data in the ERT survey which is a more suitable array for such heterogeneous site to demarcate the sub-surface weathered/fractured zones for assessment of the groundwater reserves [41]. The resistivity data were measured for 25 layers. ERT was carried out along three geophysical profiles including 101 electrodes along each profile, a profile length of 500 m and inter-electrode distance of 5 m (Figure 1). The field measurements were obtained using GPS systems (MAP60, Garmin, Olathe, KS, USA) and other surveying instrumentation. The topographic variations were measured using a clinometer along each profile with the distance interval of 10 m and less when it was necessary. The resistivity data of a single point was collected using maximum 10 stacking to improve the signal to noise ratio. A computer-based multichannel resistivity meter placed at centre of the electrode array was applied to get the electrical resistivity measurements [42]. In the ERT survey, 50 m thick subsurface formation was assessed to delineate the weathered/fractured zone depending on the hydrogeological data and the thickness of the subsurface geological strata in the studied area. For the post-processing of the resistivity data, an inversion program was performed to obtain a 2D resistivity model of the subsurface formation including the topographic relief along each ERT profile [35]. The first step to make up a 2D resistivity model is to make a pseudo section. In this step, each value of apparent resistivity is plotted on a separate section at centre of four electrodes, and at depth equal to the median depth of the investigated array [36, 40]. A least-squares technique was adopted for inversion of the apparent resistivity data [36, 40, 43]. An inversion procedure of RES2DINV software was performed to generate a 2D ERT pseudo-section along each profile [44]. This software automatically inverts the apparent resistivity data to generate a 2D resistivity model including the topographic relief [35]. A least-squares inversion technique of RES2DINV has to apply a smoothness constraint [45, 46]. The inversion procedure of the software can generate a smooth model by fitting the resistivity data to a given error level. In this investigation, RMS (root mean square) error which is the difference in calculated values of apparent resistivity and measured values of apparent resistivity values was less than 5% for all the inversion models. Such model contains a number of rectangular blocks which are equivalent to the data points in the resistivity model. The centre of the depth for the inner blocks is used as the investigation’s median depth [47]. In order to minimize the difference between the measured apparent resistivity and the modeled resistivity values, the inversion program of conventional Gauss-Newton least squares technique was applied [48]. The Gauss-Newton’s modified model [46] follows the equation:

where i shows iteration number, g_i is discrepancy vector which depends on the difference between the logarithms of the measured and calculated values of the apparent resistivity, λ_i is the damping factor, J_i represents Jacobian matrix of partial derivatives, p_i is perturbation vector to the model parameters for the i^th iteration, and C shows 2D flatness filter.

The apparent resistivity is calculated in the first step of the least squares inversion. In the second step, Jacobian matrix J is calculated. All parameters in equation (1) are solved in the third step. The first two steps are completed by applying a finite element or finite difference technique, whereas different methods including Cholesky, the modified Gram-Schmidt, and the singular value decomposition techniques are used to solve the third step [49].

3.2 Magnetic Method

The magnetic field intensity is the magnetic field generated by the North and the South Poles. The main magnetic field is controlled by the factors i.e., magnitude of the field, magnetic inclination (dip of a magnetic compass from horizontal i.e., −90^∘ for south magnetic pole, +90^∘ for north magnetic pole and 0^∘ for magnetic equator) and magnetic declination (angle between geographic north and magnetic north). The electric currents in the earth’s ionosphere produce the diurnal variations also called as magnetic storms varying from 50 to 200 gammas which can be removed by applying diurnal corrections [4].

The Ground high-precision magnetic survey was conducted along three profiles according to the Technical Specifications for High-precision Magnetic Survey on the Ground (industry standard, DZ/T 0071-93). The geomagnetic total field (nT) was observed using a GSM-19T high-precision proton magnetometer manufactured in Canada with <0.7nT accuracy and a high precision GPS system (MAP60, Garmin, Olathe, KS, USA). The sampling interval for the diurnal observations was 20 seconds. The magnetometer automatically sampled and recorded the magnetic data. The magnetic probe height was fixed at 2m to get more accuracy. The data changes were observed every 20 to 30 minutes, and there was no magnetic storm observed during the magnetic survey. However, a base station magnetometer was used to record the time variations of the magnetic field and then the changes were removed from the readings. This survey was conducted for a total of 251 measurements with 5 m station interval along three profiles (Figure 1). The magnetic data were processed to get 2D magnetic model by inverting the data using IX2D Interpex (Golden, USA) with a distance interval of 50 m along each profile. The inversion program records a response from the magnetic model. Afterwards, the response is compared with the measured magnetic data. The geometry of the subsurface geologic layers and their magnetic properties are changed constantly until the model response fits the measured magnetic data convincingly. In this way, a four layered model is generated by the inversion procedure of IX2D. The 2D magnetic model was interpreted for four layers such as the topsoil layer, the highly weathered layer, the partly weathered layer and the fresh bedrock based on magnetic susceptibility and the available upfront hydrogeological information along each profile in the study area. The weathered layer (the highly weathered and the partly weathered layer) underlying the topsoil cover and the fresh bedrock at the bottom in the model were interpreted by low magnetic susceptibility and high magnetic susceptibility respectively.

3.3 Estimation of Hydraulic Parameters

Estimation of hydraulic conductivity and transmissivity of any given aquifer system is essential for delineation of the aquifer potential zones contained within the weathered/fractured zones. Water contained within the fractures/fissures of a hard rock controls flow of the electric current in the electrolytic (mineralized) water through the ions flowing in the same pathways of water [50, 51, 52]. The electrical and hydraulic conductivity of the aquifer system are affected by the similar variables, since both depend on the potential gradients flowing from higher to lower potential [13, 17]. The pumping test naturally causes the occurrence of hydraulic potential gradient, whereas the electrical resistivity measurements generate the electric potential gradient [16]. The aquifer parameters were estimated using the pumping test and resistivity data.

The fractures/fissures are well connected because they may not exist in the real domain. In order to ensure occurrence of the interconnected fractures, the hydraulic parameters of the aquifer system were estimated from a long pumping duration. The aquifer parameters were measured using the pumping test performed at 25 boreholes (Figure 1). The double-porosity model (DP model), also called as the double continuum or overlapping continua, was used to perform the pumping test analysis [53]. This model was originally introduced by Barenblatt and Zheltov [54], and Barenblatt et al. [55] to estimate the hydraulic parameters of the weathered/fractured aquifer system. The model represents the fractured porous medium by two discrete but interacting subsystems in which one consists of the porous blocks and the other contains a network of fractures. Each subsystem is represented by a continuum taking up the entire medium domain. Hence, the interaction phenomena of two continua provide the exchange of fluid between porous blocks and fractures for a pumping test [53].

Using Darcy’s law for horizontal fluid flow and Ohm’s law of current flow in a medium, the following two relations can be derived [24, 56, 57]:

and

Where, S_c = (t/ρ) and T_r = tρ

S_c is longitudinal conductance (in mho), T_r is transverse resistance (in Ωm²), ρ is electrical resistivity of the saturated formation (in Ωm), K is hydraulic conductivity (in m/d), T is transmissivity (in m²/d), and α and β are constant of proportionality. The above equations (2 and 3) show inverse and direct relationship between hydraulic conductivity and electrical resistivity. Equation (2) is valid for an aquifer system with highly resistive basement where electrical currents tend to flow horizontally as in case of the investigated area, whereas equation (3) exists in case of highly conductive basement where electrical currents tend to flow vertically. Derivation of equations (2) and (3) is explained by the following steps:

According to Darcy’s law:

Where, q is the specific discharge (m/s) and dh/t is hydraulic gradient.

Theoretically, the electrical resistivity depends on Ohm’s law and the conservation equation of charges. According to Ohm’s law:

The flow of electric current in the medium (aquifer) is proportional to the potential gradient between the two points (source and sink). Hence mathematically, it can be expressed by:

Where dv is the potential difference, J is current density, is electrical conductivity (reciprocal of resistivity), and t is thickness respectively.

Combining equation (4) and (5) [53]:

Where the constants A_q = q/dh and A_J = dv/J describe the water flow and the electric current flow respectively, and these constants depend on the hydraulic conductivity (K).

A_J in equation (7) depends on the salinity. In the investigated area, the groundwater is fresh and there is no clay content so the constant α replaces the product of A_q and A_J:

or

Equation (9) implies that

The above equation (10) shows that the hydraulic conductivity (K) is inversely proportional to the aquifer resistivity (ρ). The above relation suggests that the low resistivity indicates the presence of fractured/fissured aquifer which is also true (valid) for the investigated area where resistivity decreases with increasing the water content in the weathered/fractured zone. Resistivity decreases with water content as well as clay content, however the studied area has no clay content except a thin topsoil layer, it implies that the low values of resistivity is only caused by the water content in the weathered/fractured a zones.

Multiplying equation (9) by t:

Since T = Kt [58]

Equation (11) is demonstrated as:

It implies that

And

Equation (14) implies that transmissivity is inversely proportional to transverse resistance for the fractured/fissured aquifer system of the investigated area.

Based on equations (10) and (14), empirical relations between the electrical parameters (ρ and T_r) and pumped hydraulic parameters (K_w and T_w) were obtained to estimate the aquifer potential over the entire area using the resistivity measurements.

The rock formation factor (F) is the ratio of aquifer resistivity (ρ_o in Ωm) and groundwater resistivity (ρ_w in m) given by the relation:

The formation factor F was estimated for the selected ERT data points near the boreholes. The groundwater resistivity required in equation (15) was calculated by:

Where, EC is electrical conductivity (in μS/cm).

The rock porosity (Φ) was calculated using Archie’s equation [59]:

Combining equations (15) and (17):

Equation (18) can be simplified as:

The coefficients a and m are related with the lithology of the aquifer system. Depending on the aquifer lithology of the studied area, the values of a and m are assumed as 1 and 2 [60, 61].

4.1 Delineation of Water Resources

The ERT modeled sections were integrated with the boreholes liyhology of the study area to get four different layers such as the topsoil layer, the highly weathered layer, the partly weathered layer and the unweathered layer. Resistivity ranges of the above layerswere obtained after calibrating the resistivity modeled sections with the boreholes data along three profiles as shown in Table 1.

The top layer consists of silt, clay and boulder with the resistivity range of 22-472 m. The next layer underlying the topsoil is highly weathered having resistivity range from 22 to 289 m. Third layer underlying the highly weathered layer is partly weathered having resistivity values between 225 and 472 m. The unweathered layer (fresh basement) is revealed below the partly weathered layer having resistivity values from 402 to 1535825 m.

The magnetic data was processed by IX2D Interpex to get 2D forward magnetic sections. The 2D magnetic modeling has been used in hydrogeology for the groundwater assessment for many years. The magnetic modeling generates a model with an interface between the weathered and unweathered layers depending on their definite geometries and magnetic properties which generate a modeled field analogous to the measured magnetic field. The magnetic data were visualized along each profile to obtain the above step. After that, a four layered model was constructed containing the top layer, the highly weathered, the partly weathered and the unweathered layers based on their magnetic properties. In this inversion program of forward magnetic modeling, a response from the model is recorded and compared with the magnetic data. Geometry of the four layered model with the magnetic properties of the subsurface layers were changing constantly until the magnetic data fitted the model response convincingly. The modeled geometry of the subsurface four layers including their magnetic susceptibility was controlled using the boreholes data that extensively improved the consistency of the subsurface four layered geological model.

A conceptual model of the hydrogeological characteristics of the subsurface in the investigated area is shown in Figure 2. The average thickness of the first three layers is about 5, 20 and 10 m respectively, whereas bottom layer is revealed at an average depth of about 30 m. The topsoil cover consists of the materials such as clay, silt and boulder. The highly weathered layer mostly contains highly weathered tuff and highly weathered fissure tuff. The partly weathered layer contains weathered tuff and weathered fissure tuff with volcano dust, feldspar, quartz, and matrix. The fresh basement mainly consists of tuff, volcano dust, quartz, matrix, feldspar, pyroxene and labradorite etc. 2D models of ERT and magnetic were integrated to delineate seven fractures/faults such as F1, F2, F3, F4, F5, F6, and F7 with NW-SE orientation along three profiles. F1 and F4 are the largest with the extension length more than 500 m, these are compressive-torsional fractures. F2 has the medium length whereas F5, F6 and F7 are small fractures/faults. Small fractures/faults were caused by the upward crust of concealed hard rock such as granite and basalt veins, resulting in the relative fragmentation of the overburden. Most of the groundwater reserves were found along the fractures/faults and the weathered/partly weathered zones.

Figure 2

A conceptual model of the subsurface geologic formations in the weathered/fractured hard rock of the investigated area.

The incorporation of 2D ERT and magnetic models along profile 1 identified four different layers including the top soil layer with resistivity range of 32-468 m, the highly weathered layer having resistivity from 32 to 289 Ωm, the partly weathered layer with resistivity range of 289-468 m and the unweathered bedrock for resistivity range from 468 to 8372 m (Figure 3a). The magnetic anomaly varies from -170 to 80 nT along this profile (Figure 3b). The subsurface resistivity and magnetic intensity show high values for the unweathered rock, and represent low values for the highly weathered/partlyweathered rock and fractures/faults saturated with water. Groundwater reserves in the Basement Complex are located in the weathered or fractures/faults zones of the hard rock system [62]. Four fractures/faults namely F1, F2, F3, and F4 were identified by the incorporation of ERT with magnetic models (Figure 3). The appropriate drilling places along profile 1 are found from 20 to 40 m, 60 to 200 m, 230 to 280 m and 330 to 360 m mainly along the faults. The integrated results are highly correlated with the boreholes data along this profile.

Figure 3

(a) 2D modeled ERT section obtained by the inversion of resistivity data along profile 1; (b) 2D forward magnetic model along the same profile.

Figure 4 of ERT and magnetic sections along profile 2 reveals four distinct layers ranging from top to the bottom as the topsoil cover having resistivity from 22 to 472 m, the highly weathered layer with resistivity values from 22 to 283 m, the partly weathered layer for resistivity values from 283 to 472 m and the unweathered layer with resistivity values as high as 153582 m and as low as 472 m (Figure 4a). The magnetic variations were measured from -60 to 80 nT (Figure 4b). Four fractures/faults (F1, F2, F4 and F5) were revealed by the integrated approach (Figure 4). The results obtained were included with hydrogeological data which shows good matching. This profile offers good drilling locations from −40 to 50 m, 100 to 230m and 300 to 380 m especially along the faults.

Figure 4

(a) 2D modeled ERT section obtained by the inversion of resistivity data along profile 2; (b) 2D forward magnetic model along the same profile.

The first layer revealed along profile 3 is the topsoil cover with resistivity ranged from 37 to 402 m, second layer is highly weathered with resistivity values from 37 to 225 m, partly weathered is the third layer with resistivity from 225 to 402 m and fresh basement is the fourth layer with resistivity range of 402 to 2520 m (Figure 5a). The magnetic anomaly varies from −27 to −12 nT along profile 3 (Figure 5b). The magnetic intensity and resistivity values along this profile are lower than other two profiles because there is no granite or basalt identified along profile 3. The integrated results of ERT and magnetic show that four fractures/faults (F1, F4, F6 and F7) exist along this profile (Figure 5). The suitable locations for drilling are suggested from -40 to 10 m, 40 m to 160 m, 220 to 250 m and 300 to 380 m along this profile.

Figure 5

(a) 2D modeled ERT section obtained by the inversion of resistivity data along profile 3; (b) 2D forward magnetic model along the same profile.

4.2 Estimation of Water Resources

In order to estimate the groundwater resources contained within the weathered and fractures/faults zones revealed by the integrated geophysical approach, aquifer resistivity (ρ_o) and transverse unit resistance (T_r) were calculated for all resistivity measurements. Based on the values of aquifer resistivity, aquifer potential was divided into four zones i.e., the high potential aquifer with ρ_o < 100 m, the medium potential aquifer with ρ_o from 100 to 200 m, the poor potential aquifer with ρ_o from 200 to 300 m and the negligible potential aquifer with ρ_o > 300 m [63]. Aquifer potential was also differentiated by transverse unit resistance (T_r). A careful observation of T_r values shows that T_r < 3300 m² reveals the high potential aquifer, T_r from 3300 to 4000 m² indicates the medium potential aquifer, T_r from 4000 to 4500 m² identifies the poor potential aquifer, and T_r > 4500 m² represents the negligible potential. The distribution of ρ_o and T_r over the entire studied area with the aquifer potential zones is shown in Figure 6. The results suggest that the high potential aquifer is contained within the fractured/fault zones (Figure 6).

Figure 6

(a) contour map of aquifer resistivity and (b) contour map of transverse unit resistance.

The groundwater resources were estimated as a function of hydraulic conductivity and transmissivity. Initially, hydraulic conductivity (K_w) and transmissivity (T_w) were determined at the specific locations of the boreholes using pumping tests analysis. In order to calculate the effective parameters for all resistivity stations to estimate the aquifer potential over the entire area, an empirical relation between aquifer resistivity (ρ_o) and pumped hydraulic conductivity (K_w) was obtained to estimate hydraulic conductivity (K’) for all resistivity measurements, and another relation was obtained between transverse resistance (T_r) and pumped trnasmissivity (T_w) to estimate transmissivity (T’) for all stations. In this way, entire area was covered for the estimation of aquifer potential based on the aquifer parameters. The empirical equations obtained from the graphical plots shown in Figure 7 are given by:

Figure 7

(a) Relation between aquifer resistivity and hydraulic conductivity, (b) relation between transverse unit resistance and transmissivity.

The values of K’ and T’ estimated from above equations (20 and 21) at the selected resistivity data points near the boreholes are given in Table 2. The pumped hydraulic conductivity (K_w) and transmissivity (T_w) measured from the pumping test are also given in Table 2. The comparison between the estimated aquifer parameters (T’ and K’) and pumped aquifer parameters (T_W and K_w) shows good matching (Table 2). The contour maps of K’ and K_w in Figure 8 provides distribution of hydraulic conductivity values estimated by pumping test and resistivity measurements. The contour maps of the estimated transmissivity (T’) and pumped transmissivity (T_w) are shown in Figure 9. The maps of estimated and pumped hydraulic parameters show good correlation.

Figure 8

(a) contour map of estimated hydraulic conductivity and (b) contour map of pumped hydraulic conductivity.

Figure 9

(a) contour map of estimated transmissivity and (b) contour map of pumped transmissivity.

The groundwater reserves estimated by hydraulic conductivity and transmissivity were characterized into four different zones. The results of hydraulic conductivity and transmissivity reveal that T’>200m²/d and K’>5 m/day delineate the high potential aquifer, T’ from 75 to 200 m²/d and K’ from 3 to 5 m/d show the medium potential aquifer, T’ from 25 to 75 m²/d and K’ from 1 to 3 m/d identify the poor potential aquifer, and T’<25 m²/d and K’<1 m/d represent the negligible potential aquifer (Figure 8 and 9). The estimation of aquifer parameters suggests that high potential groundwater reserves occur along the fractured/ fault zones (Figure 8 and 9).

The groundwater potential zones were also delineated by rock formation factor (F) and rock porosity (Φ). The estimated values of F and Φ obtained for the selected ERT data points and the nearby boreholes are given in Table 2.Acareful observation suggests that the high potential aquifer is revealed with Φ > 0.6 and F < 2.8, the medium potential aquifer is delineated with Φ from 0.54 to 0.6 and F from 2.8 to 3.5, the low potential aquifer zone is identified with Φ between 0.42-0.54 and F between 3.5-6, and the negligible potential aquifer is mapped with Φ < 0.42 and F < 6 as shown in Figure 10. Wright [63] suggested that the aquifer potential can be expressed as a function of aquifer resistivity. In the investigated area, the aquifer potential was estimated as a function of aquifer resistivity, transverse unit resistance, hydraulic conductivity, transmissivity, rock formation factor and rock porosity as shown in Table 3.

Figure 10

Delineation of aquifer potential zones based on (a) rock formation factor and (b) rock porosity

4.3 Analysis of Groundwater and Rock Samples

In order to assess to quality of groundwater contained within the weathered/fractured zones, the physicochemical analysis was performed. The physicochemical parameters of groundwater samples taken from the boreholes were analytically analyzed by the World Health Organization [64]. The results of physicochemical analysis for twenty five groundwater samples are summarized in Table 4(a). The physicochemical analysis was performed for the main anions, the cations and the parameters such as pH, total dissolved solids (TDS) and electrical conductivity (EC) as per standard procedures [65, 66, 67]. The results show that all the parameters lie within the limit suggested by WHO. The physicochemical analysis revealed that the groundwater quality is good in the study area. The groundwater quality was also assessed by geophysical analysis as shown in Table 4(a). The aquifer resistivity obtained from all ERT data points was analyzed to evaluate the groundwater quality. Generally, low values of aquifer resistivity (i.e., <25 m) suggest the saline/brackish water [1]. The geophysical analysis shows that aquifer resistivity values fall within the limit of fresh water. Hence, based on the physicochemical and geophysical analysis, groundwater quality of the investigated area is good.

The mineral analysis of twenty five rock samples taken from the borehole sites using the X-ray Diffraction Technique (XRD) was performed. The analyzed minerals were interpreted as major minerals with >50%, secondary minerals with 10-30% and minor minerals with 5-10% as shown in Table 4(b). The results suggest quartz as the major mineral, whereas kaolinite, zinnwaldite, microcline and albite as the secondary minerals, and halloysite, sodalite, perlite, and ankerite as the minor minerals in the study area.