2D joint inversion of dc and scalar audio-magnetotelluric data in the evaluation of low enthalpy geothermal fields

Abstract

Audio-magnetotelluric (AMT) and resistivity (dc) surveys are often used in environmental, hydrological and geothermal evaluation. The separate interpretation of those geophysical data sets assuming two-dimensional models frequently produces ambiguous results. The joint inversion of AMT and dc data is advocated by several authors as an efficient method for reducing the ambiguity inherent to each of those methods. This paper presents results obtained from the two-dimensional joint inversion of dipole–dipole and scalar AMT data acquired in a low enthalpy geothermal field situated in a graben. The joint inverted models show a better definition of shallow and deep structures. The results show that the extension of the benefits using joint inversion depends on the number and spacing of the AMT sites. The models obtained from experimental data display a low resistivity zone (<20 Ω m) in the central part of the graben that was correlated with the geothermal reservoir. The resistivity distribution models were used to estimate the distribution of the porosity in the geothermal reservoir applying two different approaches and considering the clay minerals effect. The results suggest that the maximum porosity of the reservoir is not uniform and might be in the range of 12% to 24%.

joint inversion, audio-magnetotellurics, resistivity, geothermics

1 Introduction

Audio-magnetotelluric (AMT) and dc-resistivity are surface geophysical methods that play major roles in shallow studies, environmental investigations and geothermal exploration (Ingham 1992, Jones and Dumas 1993, Fuji-ta et al1999, Meju 2002, Harinarayana et al2004, Monteiro Santos et al1996, 1997a, 1997b, 2006a, 2006b, El-Qady 2006). The limitations of those methods in resolving subsurface structures, mostly due to equivalence problems and incomplete and inaccurate data, are well known (Vozoff and Jupp 1975, Tuberg and Barker 1996, Monteiro Santos et al1997b).

Several authors have noted that the limitations and ambiguity of individual methods can be strongly reduced by adopting joint inversion techniques. There are several examples of works dealing with joint inversion of different measured data, mainly assuming layered earth models. Vozoff and Jupp (1975) were the first to jointly invert magnetotelluric (MT) and vertical electrical soundings (VES) in order to decrease the ambiguity of the models obtained separately from each method. A similar work was presented by Gustafson and McEuen (1987) and Monteiro Santos et al (1997b). Several other authors presented examples of one-dimensional (1D) joint inversion concerning IP, MT, AMT, TEM and dc data: Roy and Elliot (1980), Raiche et al (1985), Sandberg (1993) and, more recently, Meju (1996) and Harinarayana (1999).

Two-dimensional (2D) joint inversion of different data is not so common. Works from Sasaki (1989) inverting dc and MT data, Gallardo and Meju (2004) for seismic and dc data inversion and from Roy et al (2005) using gravity and seismic data deserve mention. The benefits of jointly inverted dc and AMT data arise from the complementarity of the data.

The aim of this work is to present the application of 2D joint inversion of scalar AMT and dc data in geothermal reservoir evaluation. The data set used in this work was collected more than ten years ago, when the use of the tensorial AMT method was not so common. This work presents a reinterpretation of some of the data using new interpretative tools.

2 The Chaves low enthalpy geothermal field—an overview

2.1 Geological setting

In the northern part of the Portuguese mainland, the most important geothermal focus is located along the NNE-trending megalineament of Verin–Chaves–Régua–Penacova. There, hot (Chaves 76 °C) and cold (Vilarelho da Raia, Vidago and Pedras Salgadas 17 °C) CO₂-rich mineral waters are present in an area of approximately 800 km² (figure 1). Most of these CO₂-rich mineral waters are used both as a source of bottled water and a recreational resource (spa facilities, tourism, etc).

Figure 1

(A) Geological sketch of the Chaves graben area (adapted from Geological Map N° 6 B, Serviços Geológicos de Portugal). The most important mineral springs on the Chaves–Verin fault (top right): 1. Vilarelho da Raia; 2. Chaves; 3. Vidago and 4. Pedras Salgadas. Main geological structures in the study area (middle): 1. quaternary sediments; 2. calk-alkaline granite (Hercynian); 3. alkaline granite (Hercynian); 4. schist (Silurian); 5. quartzites; 6. graphitic slates; 7. dykes; 8. normal faults; 9. hot spring. (B) NW–SE geological cross-section of the Chaves graben.

The Chaves region is located in the tectonic unit of Middle Galicia/Trás-os-Montes sub-zone of the Central-Iberian Zone of the Hesperic Massif (figure 1). The main geological formations are (i) Hercynian granites and (ii) Silurian metasediments of the Upper, Intermediate and Lower Peritransmontano Group, which consists of a sequence of quartzites, phyllites and carbonaceous slates (Aires-Barros et al1998). The most recent formations are Miocene–Pleistocene graben filling sediments with variable thickness, showing their maximum development along the central axis of the Chaves graben.

2.2 Thermomineral waters

The thermomineral waters flow from natural springs and drilled wells located either in granitic outcrops or in the peribatholitic boundaries parallel to the main NNE–SSW fault trend. Chaves hot CO₂-rich thermomineral waters emerge within a wide graben, whereas the cold mineral waters (e.g. Vidago/Pedras Salgadas) are found in areas where the NNE–SSW megalineament does not exhibit such an important morphological structure.

These mineral waters have a high carbon content (up to 5000 mg l^-1 of free CO₂ and 5500 mg l^-1 of HCO₃), reflecting the abundant presence of CO₂ in the system (Marques et al2000). Electrical conductivity values range from 2510 to 2890 µS cm^-1 and pH values vary from 6.87 to 7.36.

Using SiO₂ and K²/Mg geothermometers Aires-Barros et al (1998) found that the Chaves CO₂-rich mineral waters (the most representative of the deep fluids in this area) indicate equilibrium temperatures between 100 and 120 °C, which is in agreement with the issue temperatures of Chaves mineral waters. Considering the mean geothermal gradient 30 °C km^-1 (Duque et al1998), we can estimate a maximum depth of about 3.5 km reached by the water system.

3 Geophysical data

3.1 Dipole–dipole data

Various geophysical methods, mainly dc, scalar audio-magnetotelluric (AMT) and tensorial magnetotelluric (MT), were collected between 1990 and 1992, to investigate the shallow and deep structures of the Chaves graben. The dc survey comprised 29 Schlumberger vertical electrical soundings (VES), 5 dipole–dipole lines, 5 pole–dipole lines and several rectangle surveys (Monteiro Santos et al1996, 1997a, 1997b). In this work only three dipole–dipole lines placed together with scalar AMT soundings will be considered (figure 2): the DD#1, DD#2 and DD#4 lines, respectively. DD#1 and DD#2 lines cross the graben, roughly in the NW–SE direction. The E–W DD#4 line was acquired inside the graben. A 300 m spacing between electrodes was used in order to sample the deepest part of the graben. This long spacing reduced the resolution of the dipole–dipole surveys, mainly for shallow structures.

Figure 2

Location of the geophysical surveys carried out in the Chaves graben region. The locations of the sites shown in figure 5 are marked by open circles.

The main characteristic of the dipole–dipole apparent resistivity pseudo-sections (figure 3) is the low resistivity zone (10–20 ohm-m) in the central and eastern parts of the graben, underlying a relatively resistive overburden (100–400 ohm-m). All the pseudo- sections show an intense resistivity gradient limiting the conductive zone in their western part. This gradient is associated with the known Chaves–Verin fault.

Figure 3

Dipole–dipole apparent resistivity pseudo-sections (field) for DD#4, DD#1 and DD#2 lines.

3.2 Dipole–dipole data inversion

The three dipole–dipole lines were inverted separately using a regularized algorithm (Sasaki 1989). The three calculated models are presented in figure 4. The misfit error of the responses is 0.08, 0.11 and 0.08 for models #1, #2 and #4, respectively. The main characteristics of the dipole–dipole models are (1) an overburden with thickness ranging from 200 to 300 m and a resistivity range of 100–400 ohm-m, (2) low resistivity zones (10–20 ohm-m) at the central part of the graben (in the profiles DD#1 and DD#2 this low resistivity structure extends to the eastern border of the graben), and (3) an electrical basement at a depth not well defined. Towards the west the conductive zone is bound by the known Chaves–Verin fault.

Figure 4

2D resistivity models obtained from the inversion of the dipole–dipole data.

Several shallow wells were drilled into the graben, for freshwater evaluation and exploitation. None of the wells reaches the basement. The maximum depth (225 m) was attained in well AC75. The geological information from two of the wells (AC75 and ACP1) is compared with the 2D resistivity models in figure 4 and led us to conclude that the conductive zone is associated with the Pliocene–Miocene formations. The correlation is quite good in the model from the DD#1 profile and not so good in the DD#2 profile. The 300 m dipole spacing used in the dc survey is too large to allow a good resolution of the graben overburden. The low resistivity structure at the central part of the graben might be responsible for a significant skin effect limiting the resolution of the deepest structures.

3.3 Scalar AMT data

The scalar AMT survey comprises more than 100 soundings acquired in the frequency range from 2300 to 4.1 Hz (Monteiro Santos et al1996). The electric field was recorded using a 50 m dipole. The magnetic field was recorded using an induction magnetometer CM 216. The superficial geology suggests that the graben is a 2D structure approximately in the NNE–SSW direction. For this reason, in the majority of the AMT soundings the measurement direction of the electric field was N25°–30° E, whereas the measurement direction of the magnetic field was N115°–120° E. These data correspond to the called TE-mode AMT data the most adequate mode for 1D interpretation in accordance with the initial objectives of the survey.

The static-shift distortion is one of the most important problems in using MT data. The static-shift effect corresponds to the translation of the apparent resistivity curves originated by small scale and superficial structures. The static-shift effect, if present, must be corrected before the inversion of the data. Regarding the data used in this work, the agreement between the AMT and Schlumberger apparent resistivity curves is excellent (figure 5), showing that AMT data were not affected by static-shift distortions (Monteiro Santos et al1997a).

Figure 5

Schlumberger (VES) and AMT apparent resistivity curves and respective 1D joint models obtained at two sites in the Chaves graben (Adapted from Monteiro Santos et al1997b). The symbols represent the data and the lines represent the model responses.

3.4 1D and 2D inversion of the scalar AMT data

The AMT soundings collected inside the graben were inverted assuming three-layer models (Monteiro Santos et al1997b). The main characteristics of the models were an overburden layer (recent sedimentary formations, Pleistocene) with thickness varying from 50 to 300 m and resistivity ranging from 60 to 400 Ω m; an intermediate conductive layer (Pliocene–Miocene?) with resistivity varying from 10 to 27 ohm-m with thickness ranging from 200 to 900 m, and the basement at a depth range of 450 to 800 m and resistivity ranging from 130 to 600 ohm-m. The contour map of the conductance values in the intermediate conductive layer, derived from 1D models, is shown in figure 6. The conductance is a well-resolved parameter. The map in figure 6 shows that the relatively high values are located in the central part of the graben coinciding with the low apparent resistivity in dipole–dipole pseudo-sections.

Figure 6

Contour map of electrical conductance in the low resistivity layer associated with the geothermal reservoir. The conductance values were obtained from the 1D models of the AMT soundings carried out in graben.

The AMT soundings carried out coincidently with the dipole–dipole profiles were inverted using a 2D algorithm similar to the one used in the inversion of the dipole–dipole lines. Nowadays 2D inversion of both AMT data modes is quite common. Having only one of the data modes the 2D inversion should be performed using the TM-mode. However, only TE-mode AMT data were available for this work: the MT#1, MT#2 and MT#4 lines.

The calculated models are shown in figure 7. The misfit for lines MT#1, MT#2 and MT#4 is 0.12, 0.28 and 0.04, respectively. All the 2D models display a low resistivity zone approximately in the middle of the graben. These results are similar to those obtained from dipole–dipole data inversion. The more significant differences are related to (1) the dimension of the conductive zones and (2) more extreme resistivity values for conductive (<10 Ω m) and resistive structures (>1000 Ω m).

Figure 7

2D resistivity models obtained from the inversion of the AMT data.

4 The joint inversion results

The application of the nonlinear inversion technique (Tikhonov and Arsenin 1977, De-Groot and Constable 1990) to the 2D joint inversion of dipole–dipole and AMT data was presented by Sasaki (1989). In the nonlinear smoothness-constrained inversion adopted in this work the optimization equations are represented as  

1

where δp is the vector of the corrections applicable to the parameters (resistivity) of an initial model, δd is the vector of the differences between calculated and measured logarithm of the apparent resistivity (δd = y^c - y^m). J represents the Jacobian matrix whose elements are given by J_ij = ∂y_i/∂p_j. The accuracy of the measured data is taken into account by the diagonal matrix W. The superscript T denotes the transpose operation. The quantity λ is the damping factor controlling the amplitude of the parameter corrections, whose best value is empirically determined. The elements of the matrix C are the coefficients of the values of the roughness in each parameter, which is defined in terms of the neighbours.

An iterative process allows us to obtain the final model, whose response fits the data set in a least-squares sense. The misfit between data and model response is given by  

2

where N represents the number of data points.

Two questions should be clarified before the presentation of the results. The first question is related to the use of the TE-mode data in the 2D joint inversion. In his work, Sasaki (1989) used the TM-mode to perform the joint inversion of dc and MT data. In the TM-mode the electrical component of the magnetotelluric field is measured perpendicularly to the strike, which is also the direction in which the electrical potential in the dipole–dipole survey is measured. Therefore, the TM-mode data are more sensitive to lateral contrasts than the TE-mode data and seem to be the more adequate to perform dipole–dipole/AMT joint inversion. Nevertheless, in this work, only TE-mode AMT data will be used because they are the only available data. The second question is related to the fact that in line#4 the dipole–dipole profile is not exactly perpendicular to the direction of the telluric field. The angle is about 20°. As the scalar AMT data cannot be rotated, a simple theoretical study was performed to estimate the error in the model obtained using the biased dc data (see the  appendix). The result indicates that the error is compatible with the data error.

The models obtained from the joint inversion of each profile are shown in figure 8. The misfit error between experimental data and calculated responses of the models are 0.15, 0.17 and 0.14 for lines #1, #2 and #4, respectively.

Figure 8

2D resistivity models obtained from the joint inversion of dipole–dipole and AMT data.

5 Discussion

Comparing the joint inverted models with those shown in figures 4 and 7, obtained using dc and AMT data, one notes that the major differences are observed in line #4. The contribution of the AMT sounding for the joint model is more significant in this line, where the AMT soundings are more regularly spaced along the line. The conductive zone appears shallower than in the dc and AMT models, except in the Chaves–Verin fault zone, between coordinates 600 and 1200, e.g., beneath sites 38 and 45. A better definition of the resistive ‘bedrock’, between coordinates 1800 and 3600, is achieved in the joint inverted model.

A comparison of the results obtained from dc and joint inversions for lines #1 and #2 shows that the differences are less extensive and are mainly concentrated in the shallow structures. The transition between Quaternary and Pliocene–Miocene formations beneath site 12, in the line #1, is now better defined. The western limit of the conductive zone in the line #2, which, according to the geology, must be done by a vertical to sub-vertical fault, is better defined in the joint model.

Electrical conductivity is a useful parameter for characterizing porous media. The conduction of electricity through porous media occurs mainly by two mechanisms: (a) by the movement of ions through the electrolyte and (b) along the surface of pores. Thus, the conductivity of porous media depends on porosity, pore geometry, fluid saturation, fluid conductivity and surface morphology of the mineral grains. For saturated porous media, Archie (1942) proposed the following relationship between bulk (σ_b) and fluid (σ_f) conductivities:  

3

or  

3a

where a and m are parameters that are assumed constants for a certain type of rock, ϕ is the porosity and F is the formation factor. Equation (3) assumes that the surface conduction is negligible when compared with electrolyte conduction. If that is not the case, equation (3) has to be modified to include those effects. Waxman and Smits (1968) proposed the first model shaly sand formations, assuming that clay particles contribute exchange cations to the electrolyte increasing the conductivity of the formation. The Waxman and Smits (WS) model is expressed by  

4

Here, Q_v is the clay charge contribution per unit pore volume and is expressed in meq ml^-1.

Sen et al (1988) proposed an empirical law relating porosity to bulk conductivity:  

5

where μ_DL is the mobility of the ions and m_s is a surface tortuosity factor; C is a constant that depends on rock geometry and the effective mobility of cations near the surface (Sen et al1988). Sen and his co-authors found that their core data could be well fitted with μ_DLm_s = 1.93 × m Sl mol^-1, CQ_v = 0.7 S, EQ_v ≈ 0 and m ≈ 2. Thus, these values were adopted in this work.

Equations (4) and (5) were applied to the Chaves graben using the joint AMT/dc^-1 inverted models shown in figure 8 in order to estimate the distribution of the porosity in the survey zones. The calculations were performed considering an average fluid conductivity of 2600 µS cm^-1 and m = 2. The value of Q_v depends on the amount of clay present in the formation. As there is no information about this parameter, the calculations were performed for two extreme values (0.5 and 2.04 meq ml^-1). The porosity values obtained assuming the lowest value of the Q_v parameter are higher than those obtained using the highest value (88% and 75% higher for the WS and Sen et al models, respectively). Therefore, the value of 2.04 meq ml^-1 was assumed for the calculations presented in this work.

The results obtained for the three lines are shown in figure 9. The porosity distribution shown on the left was obtained using equation (5). The application of the WS model produced the distribution shown on the right. The maximum value of the porosity varies between 12 (from the WS model) and 24% (from the Sen et al model). The porosity values estimated by the WS model are lower (30% to 50%) than those calculated by the Sen et al model. The relative error in the porosity evaluation (ε_ϕ) originated by an incorrect value of the formation factor is approximately ε_ϕ = ε_F/m, where ε_F is the relative error in the formation factor. Therefore, an error ranging from 8 to 10% in the porosity evaluation is expected, considering only the error on the resistivity models.

Figure 9

2D porosity distributions obtained from the dc/AMT models using the Waxman and Smits (right) and the Sen et al (left) models. The dashed lines represent porosity values of 5% and 10%.

The relatively high porosity is restricted to the central part of the graben decreasing to south and northwards. This result agrees with the geological observations that support the existence of a secondary graben in the centre of the graben, bounded by faults roughly oriented NW–SE. Such faults, and mainly their intersections with the NNE–SSW Chaves–Verin fault, would provide an effective conduit system for fluid ascending towards the reservoir in the deep part of the graben.

6 Conclusions

AMT and dc surveys are commonly used to investigate environmental, hydrological and geothermal structures. Frequently, the separate interpretation of each geophysical data set produces ambiguous results. The combination of AMT and dc data can be used in an attempt to reduce the ambiguity inherent to each method. This paper presents results obtained from the 2D joint inversion of dc and scalar AMT acquired along three profiles crossing the Chaves graben. The joint inverted models show a better definition of shallow and deep structures, when the number and spacing of AMT sites is similar to those of dc surveys.

These inverse models of the resistivity distribution were used to estimate the distribution of the porosity in the geothermal reservoir located inside of the graben. These calculations have been done using the Waxman and Smits (1968) and the Sen et al (1988) models. Both models take into account the influence of the clay minerals over the bulk conductivity. The results suggest that the porosity of the reservoir is not uniform and might be in the range of 12% to 24%.

Appendix. Dipole–dipole data across a 2D conductive structure

A 3D simulation of a 2D conductive structure was used to investigate the differences between the apparent resistivity obtained from two dipole–dipole profiles: the first profile carried out perpendicularly to the strike of the conductive body (profile AA′ in figure A1) and the second one crossing the body in a direction of 110° (profile BB′ in figure A1). The model consists of a 2D body with resistivity of 20 ohm-m embedded in a 100 Ω m medium. An error of 1.0% is estimated for numerical calculations, which were performed using a finite element code based on Pridmore (1978). The relative difference between the two apparent resistivity pseudo-sections is shown in figure A2. The maximum difference is of about 2–4% which is less, or at least of the same order of the data errors. Assuming that the Chaves graben is a two-dimensional dominant structure, this result supports the conclusion that the main characteristics of the 2D models obtained from joint inversion of the dc and AMT data acquired along the line #4 are basically the same as those main characteristics of the models that could be obtained from the joint inversion of the AMT data and the dc data acquired along a N115-120° E profile.

Figure A1

(A) Plan view of the 3-D model used in the calculation explained in the appendix. (B) Vertical cross section of the 3D model. (C) Pseudo section of the differences between apparent resistivity obtained in the two dipole–dipole profiles marked AA′ and BB′.

Acknowledgments

The authors thank the Editor in Chief and the anonymous reviewer for their critical and helpful reviews.