3-D analysis and interpretation of magnetotelluric data from the Aluto-Langano geothermal field, Ethiopia

Abstract

The Main Ethiopian Rift Valley encompasses a number of volcanoes, which are known to be actively deforming with reoccurring periods of uplift and setting. One of the regions where temporal changes take place is the Aluto volcanic complex. It hosts a productive geothermal field and the only currently operating geothermal power plant of Ethiopia. We carried out magnetotelluric (MT) measurements in early 2012 in order to identify the source of unrest. Broad-band MT data (0.001–1000 s) have been acquired at 46 sites covering the expanse of the Aluto volcanic complex with an average site spacing of 1 km. Based on this MT data it is possible to map the bulk electrical resistivity of the subsurface down to depths of several kilometres. Resistivity is a crucial geophysical parameter in geothermal exploration as hydrothermal and magmatic reservoirs are typically related to low resistive zones, which can be easily sensed by MT. Thus by mapping the electrical conductivity one can identify and analyse geothermal systems with respect to their temperature, extent and potential for production of energy. 3-D inversions of the observed MT data from Aluto reveal the typical electrical conductivity distribution of a high-enthalpy geothermal system, which is mainly governed by the hydrothermal alteration mineralogy. The recovered 3-D conductivity models provide no evidence for an active deep magmatic system under Aluto. Forward modelling of the tippers rather suggest that occurrence of melt is predominantly at lower crustal depths along an off-axis fault zone a few tens of kilometres west of the central rift axis. The absence of an active magmatic system implies that the deforming source is most likely situated within the shallow hydrothermal system of the Aluto-Langano geothermal field.

Magnetotellurics, Hydrothermal systems, Volcano monitoring, Africa

1 INTRODUCTION

The global energy demand is ever rising and renewable energies are considered to be a major contributor to any future energy mix. Geothermal energy is an attractive candidate as it is carbon neutral and readily available in regions that may have no access to conventional energy resources. The advantage of geothermal sources compared to other renewable sources like wind, solar and hydro power is their independence from weather conditions, their constant output and their availability around the clock. Geothermal power generation is most attractive in volcanic regions with ready access to shallow high-enthalpy systems. One of the privileged regions is the East African Rift system (EARS), which possesses a remarkable geothermal potential (e.g. Endeshaw 1988). Our region of interest is Ethiopia's only producing geothermal field, Aluto-Langano, in the central Main Ethiopian Rift Valley (MER). The power plant is located on the Aluto volcanic complex where satellite observations have revealed significant surface deformation with cycles of uplift and subsidence (Biggs et al.2011). The working hypothesis is that the causative source for the deformation is either in the hydrothermal reservoir, in a deeper magmatic system or in coupled magmatic-hydrothermal system. This state of play clearly raises the question of its implication on the plant in terms of productivity and geohazard.

One way to address this problem and to delineate the nature of the deforming source is the magnetotelluric (MT) method. This geophysical method—by deciphering the subsurface electrical conductivity—is especially sensitive to high conducting zones, as hydrothermal and magmatic reservoirs usually are and it easily covers the necessary exploration depth down to approximately 10 km. The objective of this study is to investigate the electrical conductivity structure under the Aluto volcanic edifice in order to put constraints on the source of unrest and to advance knowledge about the distribution of melt by using MT.

2 TECTONIC SETTING AND VOLCANISM

The EARS is an intercontinental active rifting zone that extends for more than 3500 km from the southern coast in Mozambique to the north of Ethiopia. At its northern end, in the so-called Afar depression, it meets the Red Sea Rift and the Gulf of Aden Rift. This complex triple junction zone marks the boundary between three divergent plates, which are the Nubian plate in the west, the Somalian plate in the east and the Arabian plate in the north (Fig. 1). Rifting started in the Afar region and propagated southwards. Along its way from north to south the EARS shows all stages of continental rift evolution: from early stages of rift initiation at the southern branches to incipient continental breakup and final stages associated with the formation of oceanic crust in the Afar depression (Ebinger 2005). The Afar region has experienced extensive volcanism since 45 Ma with the onset of rifting around 30 Ma ago; here the highest spreading rates in the EARS are observed with up to 6 ± 1.5 mm a⁻¹ separation (Chu & Gordon 1999). A remnant of the pre-rift volcanism is the Ethiopian flood basalt province, which covers around 600 000 km². The volume of basalts is estimated to be 300 000–350 000 km³ with local thicknesses of up to 2000 m (Mohr 1983).

Figure 1.

Overview map on the Main Ethiopian Rift system. Black lines denote border and Wonji faults based on the fault database from Agostini et al. (2011), faults between the dotted yellow lines are considered to be part of the active Wonji fault belt (Corti 2009). White lines mark centres of Quaternary magmatic segments, the Silti Debre Zeit fault zone (SDFZ) is located west of Aluto. Geothermal fields are marked as red dots: Abaya (Ab), Corbetti (Co), Aluto-Langano (AL), Tulu Moye (TM), Gedemsa (Ge), Kone (Ko), Fentale (Fe) and Dofan (Do). Yellow stars show places with observed surface deformation. The subfigures on the lower right show InSAR interferograms taken from Biggs et al. (2011) for the uplift at Aluto (AL, 17 December 2003 to 18 August 2004) and for subsidence at Corbetti (Co, 23 September 1997 to 13 September 2000). Each fringe represents 2.8 cm of motion in the satellite line of view. The subfigure on the upper left shows an overview map on the EARS. Black lines denote border faults redrawn from Ebinger (2005). Red triangles mark volcanoes listed in the Smithonian volcano database (http://www.volcano.si.edu/). The blue rectangle marks the area within the MER as shown in the main map. Numbers 1, 2 and 3 denote the locations of the Nubian, Somalian and Arabian tectonic plates.

2.1 The Main Ethiopian Rift system

Our study area, the Aluto volcanic complex, is located in the central part of the MER, south of the Afar depression. The MER separates the Nubian and the Somalian plate, its trough is between 35 and 80 km wide and bounded by large normal boundary faults, which developed during the Miocene. The major escarpment between the rift floor and the flanking southern and northern Ethiopian plateaus is built by steep NE striking normal boundary faults and the relative displacement reaches overall differences in height of up to 2000 m (Gebregzabher 1986; Corti 2008). Today the boundary faults are considered to be inactive. A study of their morphology east of Aluto volcano by Pizzi et al. (2006) suggests a tectonically induced rather than a magma-induced faulting. Strain accommodation by mechanical stretching and formation of large border faults is typical for the young stages of rifting and is currently observed in the southern part of the EARS (Kendall et al.2006). The rift floor in the central MER is dominated by a network of active NNE–SSW oriented faults (Fig. 1). This narrow fault system developed in the past 2 Ma and is considered as present-day locus of extension. It is the youngest part of the MER and is referred to as Wonji fault belt (WFB). Main Quaternary volcanism occurs in en echelon disposed segments along the WFB within the rift valley. Seismic measurements have revealed that most of the present-day seismicity coincides with these magmatic segments and using GPS measurements it was shown that up to 80 per cent of the present day strain is accommodated across the WFB. The faults in the rift floor show recent activity with slip rates up to 2 mm a⁻¹ (Bilham et al.1999; Pizzi et al.2006). As a driving force a mixture of tectonic faulting and dyke injection is considered (Corti 2008).

Kendall et al. (2005) use observations of seismic shear-wave splitting to reveal a seismically anisotropic upper mantle with a fast shear wave direction parallel to the Ethiopian Rift. The magnitude of splitting varies significantly within short distances (50 km), which is why its origin is considered to be most likely situated within the upper 100 km of the subsurface. An increased splitting is observed within the Quaternary magmatic segments, where the fast shear-wave polarization is oriented parallel to the faults and dykes within the WFB. Among several mechanisms that could cause the observed anisotropy Kendall et al. (2006) argue in favour of oriented melt pockets like vertically aligned crack-like melt inclusions that arise from magma injection into the lithosphere and dyke injection within the magmatic segments. This is in agreement with the above mentioned dyke injection as driving mechanism for active faulting within the WFB.

A few off axis belts of Quaternary volcanism occur at the western margin of the central MER, for example, at the Silti Debre Zeit fault zone (SDFZ) which comprises the Butajira and Bishoftu volcanic chain west of Aluto. The exact extend of the SDFZ is still under discussion, initially it was considered to be ∼100 km long and 2–5 km wide (Woldegabriel et al.1990). Rooney et al. (2011) propose that the SDFZ is indeed significantly longer extending over more than 300 km from 6.5°N to 9°N.

The Aluto volcanic complex is located within the Aluto-Gedemsa magmatic segment along the WFB in the Lakes District of the MER with its the four major lakes Ziway, Langano, Shala and Abiyata. These lakes have been exposed to significant lake level fluctuations in the past and are believed to be remnants of an ancient mega lake during the Pleistocene, which has separated into the present-day isolated lake system during the late Holocene (Benvenuti et al.2002). Aluto is enclosed by Lake Ziway in the north and Lake Langano in the south, it is set close to the eastern escarpment of the rift along the WFB. Rather than being a volcano with a single eruptive centre, Aluto is considered to be built up by a cluster of overlapping rhyolithic domes, glassy lava flows, pumiceous cones and craters, which are localized along NNE tending normal faults (Benoit et al.2007). Within the Aluto-Langano geothermal field the normal faults are considered to act as an upflow zone for hot geothermal fluids. The radiometric ages of Aluto volcanic products have been determined to range from 155 000 a to 2000 a (Battaglia 2004, and references therein). At present, Aluto is in a stage of fumarolic activity focused in its inner part around the caldera structure, natural hot springs are found at the end of the southern flanks at the shore of Lake Langano.

3 SURFACE DEFORMATION AT ALUTO-LANGANO GEOTHERMAL FIELD AND MOTIVATION FOR MT STUDIES

Biggs et al. (2011) identified four areas of significant deformation in the MER, namely at Haledebi, Bora, Aluto and Corbetti by using Interferometric Synthetic Aperture Radar (InSAR) observations. InSAR is a powerful satellite technique for remote sensing of surface deformation. By comparing the phase of radar images taken at different times surface displacements of >1 cm at a lateral resolution of around 90 m can be detected. InSAR interferograms at Aluto and Corbetti show a symmetric circular deformation pattern which is centred on the volcano (Fig. 1). The largest deformation occurs at Aluto: here two cycles of rapid uplift followed by slow subsidence have been observed since 2004. The two pulses of uplift (10–15 cm) in 2004 and 2008 were followed by subsidence at an averaging rate of 3–5 cm a⁻¹ (Fig. 2). Modelling studies revealed two possible scenarios to explain the observed deformation (Biggs et al.2011): (i) a shallow (2.5 km) penny shaped crack with 3–10 km in diameter or (ii) a deeper (∼5 km) point source. The working hypothesis is that the causative source of unrest is either the shallow hydrothermal reservoir, or a deeper situated magmatic system at depths of greater than 3 km, or a hybrid of the two. However, by using solely information about surface deformation it is not possible to discriminate between a deep and a shallow source of driving force. Further information can be gained by using geophysical deep sounding techniques such as seismic, seismological and electromagnetic (EM) methods like MT. The reason for using MT to identify the source of unrest is twofold: first, MT easily covers exploration depths down to a few kilometres and second, it is especially sensitive to zones of low electrical resistivities, which, in geothermal environments, are typically related to high-temperature hydrothermal or magmatic reservoirs (e.g. Muñoz 2014).

Figure 2.

Time-series of deformation at Aluto taken from Biggs et al. (2011). Aluto shows two pulses of uplift in 2004 and 2008 with intervening phases of subsidence.

4 GEOTHERMAL PROSPECT AREAS IN THE MER AND PREVIOUS STUDIES AT ALUTO

The EARS provides a unique opportunity to meet future growth with renewable sources like hydro and geothermal energy. Geothermal exploration in Ethiopia started in 1969 but so far only limited geothermal development has taken place. In the whole MER about 120 geothermal sites, mainly surface hydrothermal manifestations, have been identified with an estimated overall potential for power generation of around 5000 MWe (Teklemariam & Kebede 2010). The highest geothermal potential is expected in the northernmost part of the MER where it shows the most advanced stage of rift evolution. However, detailed surface exploration and feasibility studies have been focused on the central MER geothermal sites due to easier access to infrastructure. Here, prioritization for geothermal exploration has been given to eight sites, namely Abaya, Corbetti, Aluto-Langano, Tulu Moye, Gedemsa, Kone, Fantale and Gedemsa (Fig. 1). Most advanced studies have been carried out at the Aluto-Langano geothermal field; currently the only producing geothermal power plant of Ethiopia is located here. During the 1980s, eight exploration wells have been drilled reaching depths of 2500 m and temperatures up to ∼350 °C (Endeshaw 1988). Six of the wells are located in the centre of the volcanic complex with two of them, LA3 and LA6, being productive. Four of them, LA4, LA5, LA7 and LA8, have been drilled into the outflow zone and show a temperature reversal with depth. The productive wells are drilled into the upflow zone, which is located under an NNE tending fault system that intersects the volcanic complex. Two other wells located west and south outside of the complex are non-productive (see Fig. 4). A pilot power plant was installed at Aluto in 1999 with an overall capacity of 7 MWe. However, due to operational difficulties it only produced 3 MWe. In 2012, a general overhaul of the power plant has started with a planned capacity of 75 MWe.

5 RESISTIVITY DISTRIBUTION IN HIGH-ENTHALPY GEOTHERMAL SYSTEMS

In search for geothermal reservoirs, in particular high-enthalpy hydrothermal systems, electrical resistivity has proven to be the most useful geophysical parameter for detecting a potential reservoir and for siting initial exploration wells. By using MT one can map the subsurface resistivity distribution easily down to depths of several kilometres. High-enthalpy hydrothermal systems of volcanic origin have temperatures exceeding 150 °C and are typically built up by a deep active or relict magmatic heat source, which is overlain by a convective hydrothermal system (e.g. Muñoz 2014). The electrical resistivities found in a geothermal reservoir are mainly governed by the hydrothermal alteration mineralogy and the reservoir fluids filling the pore space (Ussher et al.2000). The influence of fluids depends strongly on their resistivity, which in turn is determined by their salinity and to a lesser extend by their temperature. In addition, the fluid saturation, porosity and permeability of the system are of great importance as only a well interconnected network of fluids may influence the bulk electrical resistivity. A permeable system saturated with hot geothermal fluids is affected by hydrothermal alteration processes, which are epigenetic modifications of the pre-existing host rock. Common alteration products are clays that have a large effect on the bulk electrical resistivity depending on their type.

In terms of resistivity, the conceptual model of a volcanic high-enthalpy system is mainly characterized by the following features: The shallow part of the reservoir is built by a resistive surface layer depicted as zone 1 in Fig. 3. Underneath is a highly conductive cap (Fig. 3, zone 2), which is located above the more resistive reservoir (Fig. 3, zone 3). The reservoir or upflow zone is underlain by the heat source where resistivities are decreasing again (Fig. 3, zone 4). The shallow surface layer consists mainly of fresh volcanic material like tephra, pyroclastic rock or tuff, here the temperature is T < 70 °C, the water saturation is generally low and thus rocks are unaffected by alteration, which altogether, results in high resistivities in the order of 100 Ωm. The low resistivity in the deeper situated cap layer is mainly attributed to presence of conductive hydrothermal alteration products. At temperatures below 200 °C argillic alteration takes place, which most common products are smectite and illite as temperature increases. Above 200 °C propylitic alteration forms mainly chlorite (e.g. Teklemariam et al.1996). Smectite clays have resistivities below 10 Ωm whereas illite and chlorite clays have resistivities above 10 Ωm. This leads to higher resistivities in the hotter part of the geothermal system (Fig. 3, zone 3), which is covered by the conductive clay cap (Fig. 3, zone 2) in which temperatures are typically between 70 and 200 °C (Ussher et al.2000). The propylitic reservoir in zone 3 is built up by a permeable upflow zone where hot geothermal fluids migrate upwards along fractured formations and/or narrow faults driven by thermal buoyancy. In this upflow zone the temperature generally increases with depth. On top of the upflow zone is the outflow zone with a predominant lateral flow below the less permeable clay cap. A characteristics for the outflow zone is the negative temperature gradient with depth. The isotherm pattern inside such a geothermal reservoir is mainly governed by the flow pattern of the geothermal fluids (Cumming 2009). An active volcanic heat source (zone 4) underneath the reservoir is associated with the presence of partial melt or a magmatic chamber. Melt is known for having a considerably low resistivity (e.g. Pommier & Le-Trong 2011). An active magmatic heat source underneath a geothermal system is therefore identified by another good conductor at greater depths.

Figure 3.

Sketch of the conceptual model of a high-enthalpy geothermal reservoir redrawn from Johnston et al. (1992) and Cumming (2009). 1: resistive shallow layer; 2: conductive clay cap; 3: reservoir; 4: heat source. Red arrows describe the fluid flow and white lines the isotherms.

In the past the resistive propylitic reservoir (Fig. 3, zone 3) beneath the conductive clay cap (Fig. 3, zone 2) has turned out to be the most suitable target for drilling in geothermal exploration. Here temperatures are expected to be maximum and fluid saturation, porosity and/or permeability of the system are high enough to allow for productive exploitation (e.g. Muñoz 2014).

Note that temperatures inferred from conductivities are always maximum values that have been reached in the reservoir during its evolution, actual temperatures can be generally lower. The electrically high conducting properties of clay caps make them ideal targets for MT soundings. But care has to be taken as high conducting anomalies in geothermal environments are not always correlated with clay caps sitting on top of high-enthalpy reservoirs (Cumming 2009).

The described resistivity distribution in a geothermal system appears to be mainly controlled by the temperature dependent hydrothermal alteration clay products. Although fluids may contribute to the low resistivities in the smectite clay cap in high-salinity reservoirs. Towards greater depths the capability of fluids to reduce resistivities in the resistive propylitic zone is suppressed by lower porosities. Here the upflow of fluids is mainly controlled by narrow vertical fractures and faults, rather than by a broad and well interconnected network, thus there is no significant affect on the resistivity obtained by MT soundings, which average over a large volume. The hottest part of a reservoir is furthermore often steam dominated as fluids undergo a phase change with increasing temperatures resulting in the presence of a significant steam fraction as dominating phase within the large fractures and voids of the reservoir. Hereby the capability of the fluid to lower bulk resistivity is reduced as the electrical resistivity of a fluid increases dramatically with increasing steam fraction (Milsch et al.2010).

6 MT METHOD AND SKETCH OF STANDARD MT RESPONSES

MT is an EM sounding technique, which makes use of the horizontal time varying natural magnetic (B) and electric (E) field to infer the electrical conductivity structure of the Earth. Depending on their period natural EM variations originate from different external ionospheric and magnetospheric sources. Up to periods of 3 hr, the source is assumed to be laterally uniform over large regions and is consequently described as a plane wave field. This assumption can be violated in polar and equatorial regions where the uniform source field can be distorted by ionospheric electro-jet currents.

The basic MT response, Z, is defined as the linear relationship between horizontal magnetic and electric field variations at a specific site at the Earth's surface. It is complex-valued since it relates the fields in the frequency domain

\begin{equation} {\left(\begin{array}{ccc}E_x\\ E_y\\ \end{array}\right)} = \frac{1}{\mu _0} {\left(\begin{array}{ccc}Z_{xx} & \quad Z_{xy}\\ Z_{yx} & \quad Z_{yy}\\ \end{array}\right)} {\left(\begin{array}{ccc}B_x\\ B_y\\ \end{array}\right)}. \end{equation}

(1)

E_i and B_i (i ∈ [x, y]) are the electric and magnetic field variations observed along the north (x) and east (y) directions, μ₀ is the magnetic permeability of free space. Z, which is also known as the impedance tensor, is frequency dependent and it contains the information about the subsurface conductivity structure. The tensor elements Z_ij are usually displayed in terms of the apparent resistivity ρ_a

\begin{equation} \rho _{a_{ij}}=\frac{1}{\mu _0\omega }|Z_{ij}(\omega )|^2, \end{equation}

(2)

where ω is the angular frequency and phase Φ

\begin{equation} \Phi _{ij}=\arctan \left(\frac{\rm{Im}(Z_{ij}(\omega ))}{\rm{Re}(Z_{ij}(\omega ))}\right). \end{equation}

(3)

For a comprehensive introduction into MT see, for example, Simpson & Bahr (2005), Berdichevsky & Dmitriev (2008) or Chave & Jones (2012).

Additional and independent information about the conductivity structure can be derived from so-called tipper transfer functions, T, if besides the horizontal B-field variations B_x and B_y also vertical magnetic variations B_z are measured:

\begin{equation} B_z=T_x B_x+T_y B_y. \end{equation}

(4)

As a consequence of the planar wave geometry, variations in B_z (and thus T) are non zero only above non 1-D conductivity structures. T is complex-valued and frequency dependent like the impedance Z, it describes the tipping of the magnetic field out of the horizontal plain (B_z ≠ 0). Tippers are often displayed as arrows with their real parts pointing away from the conductive zone in the Wiese-convention (Wiese 1965). The length of real tipper arrows is given by

\begin{equation} \rm{Re}(T)=\sqrt{\rm{Re}(T_x)^2+\rm{Re}(T_y)^2}, \end{equation}

(5)

their direction is determined by

\begin{equation} \alpha =\arctan {\frac{\rm{Re}(T_y)}{\rm{Re}(T_x)}}, \end{equation}

(6)

where α is clockwise positive starting from North. Imaginary tipper arrows are derived analogously. As B_z is purely related to induced currents in the subsurface tipper arrows are often referred to as induction arrows. In situations with a 2-D horizontally striking conductivity contrast, tipper arrows are oriented orthogonal to the local striking direction.

7 THE MT SURVEY

The MT survey was carried out in January and February 2012 using broad-band Metronix instruments of type ADU06, ADU07 and ADU07e. Data have been acquired at 45 survey sites (ADU07, ADU07e) with a spatial resolution of around 1 km × 1 km. In addition one reference site R1 (ADU06) was set up in the south of Aluto close to the Lake Langano, it was relocated after a few days to the new position R2. The site locations are shown in Fig. 4. A good site coverage was achieved in the area surrounding the power plant in the central region of the volcanic complex, whereas it was not possible to install any sites on the inaccessible steep western and southern slopes of Aluto. The survey area has a pronounced topography with altitudes reaching from 1600 m at the Lake Langano to more than 2300 m at the peak of Aluto. The maximum altitude difference between the sites is 659 m.

Figure 4.

Map showing the elevation along with the survey site locations M1 to M45 and the remote reference sites R1 and R2. The location of the power plant is indicated as yellow star. Red dots denote the locations of wells.

8 DATA ACQUISITION, PROCESSING AND RESULTS

At the reference site time-series were permanently recorded with a sampling frequency of 128 Hz. Each of the survey sites M1 to M45 was set up for approximately 22 hr and time-series were recorded with a sampling frequency of 4096 and 128 Hz.

In order to estimate MT responses (impedances and tippers), we used the EMTF code described in Egbert & Eisel (1998) and Eisel & Egbert (2001), which allows for robust single site and remote reference processing. The responses in the ‘4096 Hz band’ have been estimated using robust single site processing for the period range 0.00098–0.00272 s. Since data in the ‘128 Hz band’ have been recorded simultaneously at the reference and the field sites we used both robust single site and remote reference processing. For periods 0.00391–0.08695 s the quality of the responses is best with single site processing, whereas for periods 0.12500–1024 s remote reference processing gained the best quality. A significant improvement could be achieved in the dead band between 0.1–10 s. Fig. 5 shows a representative plot with a comparison of different processing methods for site M12. Whereas the phases are comparable for both processing methods the apparent resistivities differ in the interval 0.2–20 s at ρ_a,yx. Here the single-site processed responses give no satisfactory results, which can be concluded from a clear violation of the Kramers-Kronig relationship (e.g. Berdichevsky 1999). With remote reference processing the unrealistic behaviour of ρ_a,yx in the dead band disappears.

Figure 5.

Apparent resistivity (top) and phase (bottom) at site M12 estimated by different processing methods. Data in the 4096 Hz band (green symbols) have been processed using a single site robust algorithm. Data in the 128 Hz band are shown for single site (blue) and remote reference processing (red).

Two further examples of ρ_a and Φ curves are presented in Fig. 6. At site M26 the curves are smooth over the whole period range and apparent resistivities of the main diagonal elements Z_xx and Z_yy are significantly smaller than ρ_a of the off-diagonal elements, which is typical for non-3-D conductivity distributions (e.g. Simpson & Bahr 2005). At most sites we obtained good responses, an example with bad responses is site M16 where quality of the responses decreases significantly for periods T > 1 s. This could be caused by, for example, insufficient electric coupling between electrodes and soil or wind induced noise due to movement of roots by a nearby tree.

Figure 6.

Apparent resistivity (top) and phase (bottom) at sites M26 (left) and M16 (right). Phases have been shifted into the [0°, 90°] quadrant.

An overview on all observed off-diagonal apparent resistivities and phases at the survey sites is shown in Fig. 7. It it obvious that both ρ_a and Φ curves are significantly scattered, which could be caused by (i) static shift, (ii) target anomalies or (iii) topography effects.

Figure 7.

Cloud plot of off-diagonal ρ_a (top) and Φ (bottom) at all survey sites. Red dots denote xy elements and blue dots denote yx elements. Responses are plotted without uncertainties. Phases have been shifted into the [0°, 90°] quadrant.

Static shift is common problem in MT studies, it is caused by shallow small scale electrical heterogeneities, so-called galvanic scatterers. The effect of this scatterers on the observed electric field E^obs is described by a real-valued frequency-independent tensor C, which acts on the undistorted horizontal electric field E^r

\begin{equation} {\bf E}^{{\rm obs}}={\bf C}{\bf E}^{r}. \end{equation}

(7)

Note that the magnetic field is assumed not to be influenced by these scatterers. As C is real-valued, the distortion of galvanic scatterers is considered to have weak effect on MT phases, but apparent resistivities are shifted to higher or lower values, which might be a possible explanation for the scattering of ρ_a in Fig. 7. However, scattering is also observed in Φ, thus it might be related to large scale geological structures with strong conductivity contrasts or topography effects (Jiracek 1990).

The tipper responses are shown in arrow representation for four different periods in Fig. 8. At shorter periods (T ≤ 32 s) they show no regular behaviour. Towards longer periods (T ≥ 128 s), they stabilize with their real parts pointing continuously to an eastward direction. At T = 1024 s, real tippers show a similar behaviour over the whole study area with their imaginary part approximately aligned to the real part. This is a clear hint on a north–south elongated conductor west of Aluto, which is sensed by the tippers as the most dominating feature at longer periods. Another contributor to the regular behaviour of tippers at longer periods might be a resistive N–S oriented block at crustal to lower crustal depths east of Aluto.

Figure 8.

Tipper responses in arrow representation at four different periods: (a) 0.0313 s, (b) 32 s, (c) 128 s and (d) 1024 s. Red arrows denote real tippers and blue arrows denote imaginary tippers. Towards longer periods the real tippers point continuously eastwards away from a good conductor in the west.

9 SOURCE EFFECT OF THE EQUATORIAL ELECTROJET ON RESPONSES

The underlying principle of MT is a laterally uniform primary source field, which satisfies the plane wave criterion. This is generally true in regions that are sufficiently far away from polar and equatorial zones. However at high and low latitudes ionospheric electrojets may produce variations that invalidate the plane wave assumption. EM variations at the Earth's surface in equatorial regions might be affected by the equatorial electrojet (EEJ), a strong electric current flowing eastwards in ionospheric altitudes along the magnetic dip equator. A well constrained and well known source signal is crucial for estimating response functions. Otherwise it may result in a so-called source effect, when transfer functions do not only contain the information about the conductivity inside the Earth but are also significantly affected by the source field geometry different from plane wave. Interpretation of source field affected transfer functions may consequently lead to serious misinterpretation (Mareschal 1986; Viljanen et al.1999).

This clearly raises the question about source field effects in MT data collected at the Aluto-Langano geothermal field located at 7°45′N, 38°45′E, which is in immediate proximity 0.5 ° south of the dip equator (Finlay et al.2010). However, Padilha (1999) have shown that MT data collected along the geomagnetic dip equator are not affected by the source field effect in a region ±3° of dip latitude and for periods up to ∼2000 s. This is valid for MT transfer functions but not necessarily for tipper transfer functions, which are known to be more sensitive to source effects (Sokolova & Varentsov 2007). To test the influence of the EEJ on tippers we estimated transfer functions for local day and night time separately in the considered period range of 0.001–1000 s. As the EEJ is a phenomenon that only occurs on the day side of the Earth, there should be significant differences in transfer functions estimated for daytime data and nighttime data separately if a source effect arising from the EEJ exists. We estimated transfer functions using single site processing at site R2 for daytime data between 9:00 and 12:00 local time and for nighttime data between 21:00 and 3:00 local time between 6 February 2012 and 13 February 2012, and for nighttime and daytime data simultaneously, the results are shown in Fig. 9. Slight differences in ρ_a and Φ between different time selections occur in the dead band around 10 s and at long periods, where the uncertainties increase due to the reduced number of available time windows. Tippers for different time selections are almost identical up to a few seconds with slight differences in the amplitude of the maximum in T_yr. Towards longer periods all time windows show a dominated T_yr-component and also a clearly increased scattering most likely due to the previously mentioned reduced number of available windows. As the observed differences in tippers cannot be definitely ascribed to variations in the power and quality of the source signal, we also estimated transfer functions at site M44 for different time windows. This site was measuring from 2012/2/27 16:00 local time for 17.5 hr. Since this site was not operating over noon we analysed nighttime data between 21:00 and 3:00 local time and daytime data from 16:00 to 19:00 local time and 6:00 to 9:00 local time. The results are shown in Fig. 10. The responses were again estimated using single site processing, the overall data quality at M44 is very good considering its relatively short recording time. Again slight differences occur in the dead band in ρ_a and Φ and towards longer periods due to the reduced amount of data, but the general agreement is good. The tippers agree very well up to 100 s, but are more scattered towards the longest analysed periods. However, their general behaviour is clearly identical over the whole period range and the dominating T_yr-component is clearly seen for all time selections. From this analysis we conclude that impedances and tippers measured at Aluto are not influenced by any source effects.

Figure 9.

Transfer function ρ_a, Φ and tippers (vtf) from top to bottom at site R2 estimated for nighttime and daytime data (left), for daytime data (middle) and for nighttime data only (right).

Figure 10.

Transfer function ρ_a, Φ and tippers (vtf) from top to bottom at site M44 estimated for the whole operating time (left), for daytime data (middle) and for nighttime data only (right).

10 3-D INVERSION OF THE MT DATA

In order to get 3-D conductivity models of the subsurface, we inverted the MT data in the period range of 10⁻³–10³ s at 41 logarithmically uniform distributed discrete periods using the ModEM code (Egbert & Kelbert 2012; Kelbert et al.2014). This modular system comprises forward and inversion schemes for frequency-domain EM geophysical data. The 3-D forward scheme is based on finite differences and the inverse scheme exploits the non-linear conjugate gradient method. ModEM can invert for a variety of MT, magnetic and magnetic intersite transfer functions, and responses can be chosen arbitrarily with respect to sites and periods. The advantage of an individual selection of responses for each site is that obviously bad data can be omitted without excluding a complete set of responses over the whole period range. Model regularization is implemented in a flexible manner, which allows us to choose different sets of smoothing and trade-off parameters. The implemented spatial smoother is a recursive autoregressive covariance operator with a quasi Gaussian smoothing kernel (Purser et al.2003). The smoothing parameter has to be chosen from the interval [0, 1] whereas higher values result in smoother models. The impact of the smoothing operator on model smoothness is illustrated in Fig. 11 for different smoothing parameters. In ModEM the model domain significantly exceeds the area of interest, which is necessary to obtain accurate forward solutions. In addition the forward and inverse model domains coincide. The prior model corresponds to the starting model. Note also that including topography is allowed in ModEM.

Figure 11.

Plot showing smoothed model vectors m(α) of a starting model m_initial using the smoothing operator implemented in ModEM for different smoothing parameters α.

The overall modelling area for our studies is chosen to be 406 km (N–S or x) by 402 km (E–W or y) by 271 km (depth or z). Embedded into the large model domain is the area of interest, or core, with 50 cells in x and 38 cells in y direction, each cell has a horizontal size of 400 m to 400 m, thus resulting in a core with 20 km in the x and 15.2 km in the y directions. We implemented topography by vertically discretizing the topography data for the area of interest into 14 layers with 50 m thickness each. To greater depths cells are continued with a subsequently increasing thickness by a factor of 1.2 until a total cell number in vertical direction of 54 (271 km). In the horizontal direction towards the margins cell sizes increase by a factor of 1.3. The total number of model parameters sum up to 86 × 74 × 54 = 343656. All inversions were run on 2 * N_per + 1 processors, where N_per denotes the number of periods, in our case this results in a total of 83 processors.

In order to evaluate the quality of the final recovered model, we analyse the root mean square misfit (RMS) and the distribution of the residuals r_i. The overall RMS is given by ModEM for each iteration and is calculated as follows

\begin{equation} {\rm RMS}=\sqrt{\frac{1}{N}\sum _{i=1}^N\frac{\left(d_i^{{\rm obs}}-d_i^{{\rm pred}}\right)^2}{e_i^2}}, \end{equation}

(8)

where |$d_i^{{\rm obs}}$| and |$d_i^{{\rm pred}}$| are observed and predicted responses and e_i are the errors of the observed responses, N counts the responses at all considered sites and at all considered periods. The residuals are estimates of the errors derived from the differences between observed and predicted responses and are calculated as follows

\begin{equation} r_i=\frac{d_i^{{\rm obs}}-d_i^{{\rm pred}}}{d_i^{{\rm obs}}}. \end{equation}

(9)

The distribution of residuals provides important insight about the observed data, which are not explained by the fitted model. Recovered models are considered to be appropriate if the residuals are at least distributed with a mean of 0 and if their variance is small.

In order to derive a subsurface resistivity model we inverted the full impedance tensor. We did not apply any corrections for static shift effects as ModEM accounts for these effects by introducing a scattered conductivity distribution in the near surface layers to model distortion effects (Tietze & Ritter 2013; Meqbel et al.2014). As prior and starting model we used a homogeneous half-space. The resistivity of the half-space (20 Ωm) was chosen to have minimum starting RMS, compared to alternative starting models. We chose a smoothing parameter of 0.3 in the horizontal and vertical directions, smoothing between land and air layers is turned off in the model covariance file. We run inversions with different smoothing parameters reaching from low smoothing of 0.1 to high smoothing of 0.7. Hereby we figured out that 0.3 is the best choice to allow for a balanced interplay between model roughness and smoothness (compare Fig. 11).

The data used for inversion were first visually inspected and obviously bad responses were omitted. Plots of the apparent resistivity and phase of the selected data at all sites are shown in the Appendix A (Figs A1–A6). The data were not rotated as the coordinate system of the measurement already corresponds to the geoelectric strike direction inferred from tippers at long periods (see Fig. 8). We used the same data selection regarding the period range for Z_xx as for Z_xy and the same for Z_yy as for Z_yx, the assigned errors are 5 per cent of |Z_ij| and higher if the true errors exceed 5 per cent. The starting RMS is 8.28, it decreased to 1.43 within 52 iterations. The predicted responses are shown in the Appendix A (Figs A1–A6) at all sites in terms of apparent resistivity and phase.

The distribution of the normalized residuals for the starting and final inversion run is shown in Fig. 12 separately for diagonal and off-diagonal elements of Z. It can be seen that the variance of the off-diagonal residuals for the final iteration is significantly smaller and the residuals are symmetrically distributed around zero. The distribution of the diagonal residuals shows one peculiarity: for the initial iteration the residuals are centred around 1. This is due to the fact that our starting model is basically 1-D except for 3-D effects arising from topography. The way to calculate the residuals (eq. 9) consequently implies a starting residual of 1 as the |$d_i^{{\rm pred}}$| of the starting model are zero for the diagonal elements except where we have 3-D topographic effects. The final distribution of diagonal residuals is less distinct centred around zero compared to the off-diagonal residuals. This reflects the fact that the diagonal elements are orders of magnitude smaller than the diagonal elements, respectively close to zero. Thus even a small deviation between observed and predicted data result in a larger absolute misfit. However, we conclude from the residuals distribution that the final model gives a plausible representation of the subsurface's conductivity distribution.

Figure 12.

Distribution of initial and final residuals for off-diagonal elements of Z (left) and diagonal elements (right).

The horizontal slice of the final model (Fig. 13) clearly shows the influence of the non-uniform site distribution especially at shallow depths, where sensitivity is only given close to the site locations. To avoid misinterpretations we show vertical slices of the 3-D model only along sections with a good site coverage denoted as AA′, BB′ and CC ′ in Fig. 13. Slices are shown to a maximum depth of 8 km, since there are no significant structures towards greater depths. The vertical slices along AA′, BB′ and CC ′ are shown in Figs 14(a)–(c). The vertical slice along AA′ shows a low resistive zone (ρ < 5 Ωm) at depths down to ∼1.5 km in the central region of Aluto (7 km < y < 14 km), the shallow layers in the upper few hundred meters are resistive (ρ ≥ 100 Ωm). At greater depths (below 2 km) resistivity increases again reaching maximum values of ρ > 100 Ωm at depths of around 5 km at 11 km east (Fig. 14a). The E–W slice along BB′ south of Aluto shows a pronounced shallow conductor in the central part between 7–12 km east, which is only around 100 m thick and situated at depths between 100–200 m below surface. Between 2–6 km east the shallow conductor appears more blurred, which could be an artifact of the inversion as the MT sites are 1–2 km north of the section BB′ in this part. Below the shallow conductor resistivities are between 10–50 Ωm (Fig. 14b). The low resistive zone identified in AA′ is also found in the N–S slice along CC ′ (Fig. 14c), here it starts at the central part of Aluto at 13 km north and ends at around 17 km north, as in the case of section AA′ it reaches depths down to ∼1.5 km. At greater depths this zone is followed by a resistive body that clearly exceeds resistivities of 100 Ωm. The recovered model shows no significant structure at depths below 8 km, here the resistivity is close to the 20 Ωm we have chosen in the starting model.

Figure 13.

Horizontal slice at z = 1455 m.a.s.l.. Projected surface site locations are indicated as black squares. The depth of the horizontal slice corresponds to the first subsurface layer without topography.

Figure 14.

Vertical slices along AA′ (a), BB′ (b) and CC ′ (c).

Fig. 15 shows a 3-D volume plot of the central part of the volcanic complex where we have a good site density. The conductive cap with ρ ≤ 5 Ωm extends over the the central area and overlies the resistive deeper structure with ρ ≥ 80 Ωm. The approximate location of the productive wells LA3 and LA6 are marked in the figure and it can be seen that they are located above the conductive cap and the deep roughly N–S oriented resistive body. Well LA3 is drilled down to a depth of 2144 m below the surface and LA6 to a depth of 2201 m. Both wells are thereby tapping the upper part of the resistive body (see also Fig. 16).

Figure 15.

3-D volume plot showing the central part of the volcanic complex. As values for the shown isosurfaces we chose 80 Ωm for the resistive deep structure and 5 Ωm for the conductive shallow structure. The locations of wells LA3 and LA6 are marked as triangle, white for LA3 and LA6 with a checkerboard pattern.

Figure 16.

EW section with a overlay plot redrawn from fig. 3 in Gianelli & Teklemariam (1993). The isotherms are based on interpolated well data. Locations of wells are shown by the black triangles. The white lines mark the location of major faults.

Note that the recovered resistivity structures appeared to be stable also using different starting models, grid discretizations and model constraints. Besides the starting model with a homogeneous resistivity of 20 Ωm we tested starting models with 50  and 100 Ωm. We also used a starting model with a 1-D resistivity section inferred from inverting the median of ρ_a,xy,yx at all sites, however this model gave no improvements to the fit compared to the inversion using a 20 Ωm starting model. None of the models we obtained in various inversion runs could recover any distinct structures below the resistive body in the centre of the volcanic complex at depths greater than 8 km, where resistivities stay at moderate values close to the chosen starting model. For comparison, at the Tendaho geothermal field in NE Ethiopia a distinct conductor at depths z > 5 km was recently interpreted by Didana et al. (2015) as being correlated with partial melt.

11 INTERPRETATION OF THE CONDUCTIVITY STRUCTURE

The recovered 3-D resistivity model shows a highly conductive zone with a diameter of around 6 km in the centre of the Aluto volcanic complex at depths between ∼500–1500 m and a resistive body at depths below 2000 m underneath the conductive cap. The identified resistivity distribution corresponds very well to the conceptual reservoir model of a high-enthalpy geothermal system with its mushroom shaped zoning comprising a resistive upflow zone overlain by a conductive clay cap as discussed in Section 5 and illustrated in Fig. 3. From well data analysis it is known that the conductive cap at Aluto is related to an argillic hydrothermal alteration zone with temperatures below 200 °C and marked by the abundance of calcite and clay minerals with an estimated overall amount of 20 per cent of smectite in the illite/smectite mixed layer (Teklemariam et al.1996). The clay cap follows mainly the fluid flow pattern of the outflow zone. Fig. 16 shows the details of the EW section along AA′ as recovered by the inversion with an overlay plot taken from Gianelli & Teklemariam (1993) showing the isotherm pattern as interpreted from interpolated well data. It can be seen that the conductive clay cap is situated in the temperature zone below 200 °C. Minimum resistivities are found at LA3 were fluids migrate upwards along a NNE oriented fractured fault zone which is part of the WFB. By gravity data analysis Saibi et al. (2012) could clearly identify the WFB as highly fractured area intersecting the volcanic complex at the location of the productive wells LA3 and LA6. However, for a clear identification of the geometry of individual faults in our model a higher resolution would be necessary, the site spacing of ∼1 km is too large and not sufficient to resolve such narrow vertical features.

The two productive wells LA3 and LA6 are drilled into this zone, they show a very small temperature gradient for depths below 1000 m, which is a typical feature of upflow zones, see plots for LA3 in Fig. 17. In LA3 the temperature profile exceeds the boiling point pressure curve showing that the upflow zone is clearly steam dominated with a dominating gas phase at the well head as reported by Gizaw (1993). All other wells exhibit temperature inversions with depth indicating that they were drilled into the outflow zone, see plots for LA5 in Fig. 17. Note that the shallow conductor that appears at LA3 in 200–1000 m depths does not appear as distinct at LA5. This is in accord with the lacking abundance of conductive smectite and illite clay minerals in LA5 (Teklemariam et al.1996).

Figure 17.

Left: resistivity-to-depth sections from 3-D model recovered by the inversion at locations of wells LA3 and LA5. Central: temperature versus depth distribution in wells LA3 and LA5, data were taken from Gizaw (1993). Right: abundance of clay minerals in LA3 (Teklemariam et al.1996) (ILL: Illite, SM: Smectite, CHL: Chlorite).

The deep resistive body could be identified as zone of upflow and hottest part of the geothermal field characterized by maximum temperatures up to 363 °C and the presence of propylitic alteration. Here most common alteration products are chlorite, epidote, calcit and quartz as major phases. We did not find any deep low resistivity structure as one would expect in case of a hot extended magma reservoir under Aluto. The heat source seems therefore to be a post-magmatic relic of a former active magmatic system.

The clay cap seems not to extend to the south of Aluto, where only a thin shallow conductive layer is found, which is most likely related with a groundwater aquifer fed by highly saline waters of Lake Langano. Well LA1, which was drilled in this area, revealed only one permeable zone at depths corresponding to the shallow conductor. The maximum temperatures in this well are slightly exceeding 100 °C and the pressure profile is close to cold-hydrostatic (Gizaw 1993).

As we found no evidence for an extended active magmatic system under Aluto we conclude that the source of unrest obtained by InSAR studies is most likely related to changes in the hydrothermal regime of the geothermal system in the shallow argillic alteration zone and in the deeper propylitic alteration zone in the central part of Aluto. Possible explanations and mechanisms will be addressed in Section 14. In the next Section we will provide a forward model in order to explain the observed tippers and to put constraints on the possible occurrence of melt within in MER. It is remarkable that the tippers reveal an N–S elongated highly conductive structure west of the geothermal field. The distance to this conductor can be restricted to several tens of kilometres, based on the penetration depth for periods between 100 s and 1000 s.

12 MODELLING THE TIPPERS

Inverting impedances and tippers simultaneously revealed that the regular eastward pointing behaviour of the real tippers at long periods (Figs 8c and d) cannot be explained by any structure inside the modelling region with its horizontal extent of 20 km in N–S to 15 km in E–W. However, it is clear that the regular behaviour of the tippers towards long periods is caused by an N–S elongated conductive structure at lower crustal depths west of Aluto, respectively a resistive structure east of Aluto or a combination of both. The only known structure coming into question as potential conductor is the SDFZ 40 km west of Aluto at the western rift margin (see Fig. 1). Here seismic studies analysing the crustal S-wave velocity structure have revealed the presence of a significant amount of melt. In contrast to the WFB magma bodies at the SDFZ extend throughout the crust down to depths of >20 km (Kim et al.2012). MT studies by Whaler & Hautot (2006) along an N–W S–E profile crossing the MER ∼120 km N–W of Aluto revealed a conductor (ρ = 1.5–2.5 Ωm) at depths between 20 and 25 km apart from the western rift margin. They related it with the Bishoftu volcanic chain, which is not directly intersected by their profile that crossed a few tens of kilometres north of it. In the same study Whaler & Hautot (2006) found out that resistivities in upper and lower crust beneath the southern plateau east of the rift are with about one order of magnitude considerably higher than beneath the northern plateau west of the rift, this might also contribute to the observed behaviour of the tippers.

In order to analyse the influence of a possible N–S elongated conductor underneath the SDFZ and a resistive lithosphere underneath the southern Ethiopian plateau on the tippers at Aluto, we constructed a minimum structure forward model referred to as ALTIP with a modelling core that spans 100 km in E–W and 80 km in N–S direction. ALTIP consists of a conductive structure representing the SDFZ embedded into a 25 Ωm homogeneous half-space. The vertical extent of the conductor reaches from 5 to 28 km depth, in E–W direction it is 4 km wide at its top and with 12 km wider at its lower end, in N–S it is elongated over the whole core. For its resistivity we assigned 5 Ωm, a typical value which was found for magmatic bodies in the Afar region (Desissa et al.2013). Additionally we added a resistive lithosphere with ρ = 250 Ωm in the eastern part of the model to mimic high resistivities beneath the southern plateau as obtained by Whaler & Hautot (2006). A vertical E–W slice of the ALTIP model core is shown in Fig. 18.

Figure 18.

Vertical slice showing the core of the ALTIP forward model to recover the influence on tippers of a N–S elongated conductor at the SDFZ 40 km west of Aluto and a resistive lithosphere east under the southern plateau (SP) of Aluto. The resistive part extends to the model boundary. The location of Aluto is marked as red triangle.

Observed and predicted tippers of the model are shown in Fig. 19 for the period of 1024 s. It can be seen that the predicted tippers have a T_x component of zero as a result of the purely 2-D conductivity distribution of ALTIP. Structures are only varying along the y direction resulting in strictly E–W oriented tipper arrows. Regarding their amplitude at least the predicted real tippers agree well with the observed tippers, whereas imaginary parts are slightly underestimated. Fig. 20 shows a comparison of observed and predicted tippers for periods T > 1 s. Again it is clear that the T_x component of predicted tippers is zero over the whole period range due to the above mentioned reasons. Thus the predicted T_x component does generally not fit the slight variations in the corresponding observed component. The T_y component of the predicted tippers matches the overall behaviour of observed tippers quite well for periods T > 30 s. Especially the agreement in the real parts is encouraging, whereas ALTIP does not reproduce the slightly undulating observed imaginary parts at longer periods. At shorter periods the forward model does not explain the observed tippers as we omitted to include any small scale structures in the vicinity of Aluto in order to focus on the analysis of effects arising from the proposed large scale structures under the SDFZ and under the Somalian plateau. These structures affect the predicted responses only at periods T > 30 s as can be seen from Fig. 20, at shorter periods both imaginary and real parts of T_x and T_y are zero.

Figure 19.

Observed (left) and predicted tippers of model ALTIP shown in Fig. 18 at a period of 1024 s.

Figure 20.

Comparison of observed and predicted tippers of model ALTIP with respect to period at four representative sites reaching from north to south: M22, M44, M27 and M42. The plots on the left show the T_y component, the plots on the right show the T_x component. See Fig. 4 for site locations.

The tipper modelling studies show that the minimum structure model ALTIP featuring the N–S elongated conductor west of Aluto and the resistive lithosphere under the southern plateau is able to explain the main behaviour of the tippers at long periods in a large part. Thus we think that model ALTIP is, despite its simplicity, a step in the right direction as it is well constrained by data of previous magnetotelluric and seismic studies (Whaler & Hautot 2006; Kim et al.2012). However, more regional data would be necessary in order to recover the local 3-D conductivity distribution in the MER by investigating the extent, geometry and true conductivities of the conductor under the SDFZ and the lithosphere under the southern plateau.

13 INFLUENCE OF THE CONDUCTOR UNDER THE SDFZ ON IMPEDANCES AT ALUTO

The tipper modelling study presented in the previous section showed that the tippers behaviour is described by a N–S elongated conductor under the SDFZ. The influence of this offside conductor on the impedances at Aluto remains to be analysed. In order to do so we embedded the conductor into the starting model that we already used for the inversion described in Section 10. The conductor was fixed during the inversion, in other words, no smoothing was allowed between the conductor and the surrounding homogeneous 20 Ωm half-space. For this inversion we used the same data selection as for the previous one. The inversion with the fixed SDFZ conductive feature converges slightly faster, after 39 iterations to an final RMS of 1.34, which is an improvement compared to the inversion without offside conductor. The offside conductor slightly improves the fit of the impedances at long periods as can be seen in Fig. 21, whereas no significant changes occur in the shallow part of the recovered resistivity model.

Figure 21.

Observed and predicted apparent resistivity and phase at site R2 for the inversion run with fixed conductor under the SDFZ. In comparison to the inversion run without offside conductor it can be seen that the fit at long periods T > 100 s slightly improves, especially in ρ_a,xy (compared with Fig. A1).

14 CONSTRAINTS ON THE SOURCE OF UNREST

In this section we will present two possible explanations for the source of the observed unrest, which are both based on major changes over time in the hydrothermal regime of the geothermal field and complex water-rock interaction processes. It is very likely that fluids act as causal agent driving kinematic mechanism that finally results in the observed ground level oscillations. The first mechanism is related to the hydro-mechanical behaviour of clay minerals and their tendency to swell and shrink when exposed to changes of water saturation and pore water chemistry. The second mechanism we refer to is thermoelastic expansion of fractured rock consequent to forced advection of hot fluids. Hereafter we first introduce the two mechanisms and subsequently discuss possibilities how they might initiate the kinematic mechanisms which lead to the observed unrest at Aluto.

14.1 Clay swelling

Illite and smectite clays, which account for the low resistivities observed within the clay cap, are layered aluminosilicates that are built up by negative charge layers and hold together by positively charged interlayer counterions. These counterions have a great tendency to hydrate and to include water molecules between layers inducing the swelling of clay aggregates, which consequently exhibit a macroscopic swelling response. Studies of clay swelling in clay-water-cation systems under crustal pressure and temperature conditions showed that the layer spacing of dry smectite-rich clay samples increases by a factor of two due to hydration under pressure conditions corresponding to depths down to 6 km (de Siqueira et al.1999). Clay swelling is a reversible process and several cycles of swelling and shrinking can occur whereas the swelling ability may decrease or increase depending on whether a clay sample is fully hydrated and dehydrated or not (Basma et al.1996; Tambach et al.2006). Clay swelling is not only controlled by the water content itself but also by the pore water chemistry such as the pH and ionic strength of the water (Laird 2006). Xu et al. (2006) developed a simulation program to study the impact of reinjected waters on clays in fractured geothermal reservoirs. They report that clays may swell to a multiple of their original volume when reinjected geothermal water is diluted with freshwater. However, it is possible that clays shrink back when the ionic strength of the diluted water increases again due to evaporation, decreasing freshwater supply or manipulation of the reinjected water. Clay swelling under freshwater exposure leads to a modest decrease in porosity as minerals precipitate and form a porosity reducing gel, this results in obstruction of flow paths and a large reduction of the overall permeability.

14.2 Thermoelastic expansion

Bonafede (1991) suggests a mechanism driven by hot geothermal fluids to explain periods of uplift and setting at the Campi-Flegrei volcano in Italy. A sudden formation of flow paths by opening faults allows hot pressurized geothermal fluids from greater depths to migrate upwards. Fluids raise relatively fast into the shallower and cooler low-pressure part of the reservoir below the clay cap. This efficient heat transfer increases the host rock temperature, which leads to thermoelastic expansion. Bonafede (1991) showed that forced advection of hot fluids is an efficient source of ground deformation as the associated thermal and pressure uplift easily explains a ground deformation in the order of a few tens of centimetres. Lima et al. (2009) provide a thermodynamic model for the uplift and deflation at Campi-Flegrei by episodic upflow of magmatic fluids that have been exsolved from cooling magma during crystallization, in their model fresh recharge of magma is not a prerequisite. In a more recent study Troiano et al. (2011) succeeded to fit ground deformation data with consecutive cycles of uplift and subsidence at Campi-Flegrei by numerical modelling of deep geothermal fluid injection into a shallow geothermal system. Their numerical approach allows for modelling multidimensional flow of multiphase geothermal fluids in a porous media initiating thermoelastic deformation. The numerically recovered inflation/deflation curves for Campi-Flegrei adequately match the observed deformation pattern using realistic permeability values and injection rates consistent with observed microgravitiy changes. Troiano et al. (2011) related small uplift events in the range of 10 cm, comparable to the ones observed at Aluto, to occasional permeability enhancements by fracture opening in the rock matrix forcing deep hot fluids to migrate upwards. This advective heating events are relatively brief and followed by longer cooling events characterized by setting. Cooling starts either due to exhaustion of the deep source or consequent to obstruction of conduits that allow hot fluids to migrate upwards. However, in their study the question in how far magmatic processes are involved in the genesis of unrest remains unanswered.

14.3 Mechanisms of unrest at Aluto

The introduced mechanisms could give rise to the unrest at Aluto being triggered by stress accommodation along the active faults of the WFB that intersect the volcanic complex, this is discussed in the following. From seismic studies it is known that the region around Aluto shows a diffuse seismicity pattern. Epicentres are widely distributed within a ∼30 km wide zone and down to depths of 20 km along the Aluto-Gedemsa magmatic segment (Keir et al.2006). Increasing tensile stress leads to seismic activity and the opening of fractures and flow paths. First, this allows groundwater and lake water to infiltrate laterally into shallower part of the geothermal system. Inflowing freshwater and enhanced circulation of fluids in the shallow part at depths of the smectite/illite cap will provoke clay swelling with major affects on the ground level depending on the affected volume. With time consequent swelling reduces the permeability and flow paths are finally closed again. This would correlate with the end of an uplift event. Subsequently evaporation and fumarolic activity leads to dehydration and an increase in ionic strength of the pore water ultimately resulting in clay shrinking. This period would correspond to the observed periods of setting in the InSAR time-series of deformation until the next rifting event allows new flow paths to open within the clay cap. Second, the opening of deep faults would simultaneously enforce hot geothermal fluids that have been exsolved from crystallizing magma to move upwards from a proposed deeper reservoir within the fractured upflow zone leading to thermally induced elastic expansion in the upper part of the reservoir. The end of an uplift in this scenario would correspond to the exhaustion of the deep source or the obstruction of flow paths that are, for example, resealed due to mineral precipitation.

Evidence for both scenarios is given in previous studies about the Aluto-Langano geothermal field. The first one of clay cap swelling by freshwater infiltration is supported by studies of Ayenew (2003), who found evidence for a strong hydraulic link between the lakes Ziway and Langano based on isotope studies. The groundwater flow follows active NE-SW tending faults and fractures that partially intersect the Aluto volcanic complex and is driven by the potential gradient between Lake Ziway (1636 m.a.s.l.) and Lake Langano (1582 m.a.s.l.). Isotope studies of rift waters by Rango et al. (2010) confirm that the faults of the WFB facilitate a fast water circulation that feeds hot springs, groundwater and geothermal wells around the Aluto volcanic complex. They found direct evidence for lake water mixing in shallow wells that were drilled for geothermal gradient studies around Aluto, which supports the subsurface hydraulic connection between the lakes. A study of hot spring waters further reveals mixing of deep geothermal fluids with lake water and meteoric water. The reported mixing of fluids in the shallower part of the geothermal system reveals a dynamic hydrologic regime where changes in water chemistry are likely to occur and will lead to clay swelling and shrinking.

The second scenario of thermoelastic expansion is evidenced by analyses of well fluids and fluid inclusions, which report a continuous and ongoing dynamic evolution of the hydrothermal system with important mixing processes of fluids from shallow and deep sources that even affect the deepest part of the system. Teklemariam et al. (1996) found out that argillic alteration in the lateral outflow zone actually overprints and replaces a previous zone of propylitic alteration. Trapping temperatures recovered from fluid inclusions exceed recent temperatures in the outflow zone by up to 170 °C indicating that a major change with significant cooling happened during the evolution of the geothermal system. This is confirmed by the occurrence of epidote in the outflow zone at temperatures below 250 °C in the wells LA2, LA4, LA5 and LA7, which is only formed at higher temperatures that must have existed previously (Gianelli & Teklemariam 1993). On the opposite direct evidence for recent heating by fast upflow of hot fluids is found in the upflow zone tapped by wells LA3 and LA6. Here the in-hole temperatures are almost 100 °C hotter than Valori et al. (1992) concluded from studies of fluid inclusion trapped in hydrothermal alteration products in the same wells. Such an increase in temperature is very likely to implicate thermoelastic expansion within the volume penetrated by the hot fluids. Valori et al. (1992) explain the observed heating by recent intake of a high-enthalpy fluid rising from greater depths with a temperature of at least 350 °C. Well LA7, which is nowadays used as an injection well, is around 1.1 km west of the identified upflow zone. This well is affected by the inflow of cooler waters from a shallow aquifer that result in a temperature reversal and a depression from 225 °C at 800 m depth to 125 °C at around 2000 m depth (Gianelli & Teklemariam 1993). The huge differences in the hydrothermal regime between neighboring wells LA3, LA6 and LA7 illustrate the complexity of the flow pattern and the dynamics of the hydrological system at Aluto. Complex fluid mixing processes have been reported in numerous hydrogeochemical studies (Gianelli & Teklemariam 1993; Teklemariam et al.1996; Pürschel et al.2013): the common conclusion is that all of the present day well fluids are immature waters with none of them being in full equilibrium with the host rock. Significant mixing takes place over a wide area in the geothermal field with dilution factors defining the ratio between shallow and geothermal waters of up to 10 in the outflow zone. Mixing even affects the deep hottest part of the reservoir where a deep running fault system allows migration and interaction of fluids from a deep thermal reservoir with meteoric water from shallower sources. The described circulation of fluids geothermal fluids is illustrated in Fig. 22.

Figure 22.

Sketch illustrating the origin of fluids, fluid circulation and mixing that may lead to the observed unrest. Red arrows indicate upward flow of deep hot fluids, blue arrows show the lateral infiltration of water from Lake Ziway into the clay cap along opening faults. Hot fluids lead to thermoelastic expansion in the upper part of the reservoir, cooler fluids lead to clay cap swelling in the outflow zone. Well LA7 is drilled into the outflow zone where hot fluids mix with cooler fluids of the upper aquifer, LA3 is drilled into the hot upflow zone.

We believe that the reported fluid driven dynamics in the hydrologic regime of the geothermal field might trigger the kinematic mechanisms of clay swelling and thermoelastic expansion respectively shrinking and will consequently result in the observed cyclic surface deformations.

15 CONCLUSIONS AND OUTLOOK

The MT measurements at the Aluto-Langano geothermal field revealed a 3-D electrical resistivity distribution corresponding to expected resistivities of the conceptual model for high-enthalpy volcanic geothermal systems. The conductive clay cap is overlying the resistive propylitic upflow zone as confirmed by the well data. The clay cap could not be identified south of Aluto, where only a thin shallow conductive layer is found, which is most likely related to a groundwater aquifer fed by saline water of Lake Langano. The recovered 3-D resistivity model shows no distinct electrically conductive structures below the resistive upflow zone, as one would expect if a magmatic chamber with the presence of melt would exist under Aluto. Based on these results we conclude that the presence of a hot extended magma reservoir under the Aluto volcanic complex is unlikely and that the source of unrest is most likely related to changes in the hydrothermal system rather than to changes within an active magmatic system. As possible mechanisms resulting in cyclic surface deformations we suggest processes of thermoelastic expansion and clay swelling triggered by the complex and dynamic fluid flow within the reservoir with mixing of shallow and deep fluids. The dynamics of the fluid flow are driven by seismic activities with permanent formation and obstruction of flow paths related to the active fault system of the WFB intersecting the volcanic complex. Major changes over time in the hydrothermal regime of the geothermal system have been reported before, mainly from studies during the 1990s (e.g. Valori et al.1992; Teklemariam et al.1996). However, we have no recent data that could be compared with the observed uplift events in 2004 and 2008.

A forward model to predict the tippers at Aluto suggests a conductor at lower crustal depths along the SDFZ 40 km west of Aluto and a resistive lithosphere under the Somalian plate. This conductor corresponds to a seismic low velocity zone (Kim et al.2012) and evidence regarding its existence was already found in a previous MT study (Whaler & Hautot 2006). However, so far no MT data have been measured at the SDFZ itself which is why there is no information about the local resistivity distribution and the exact location, shape and extend of the proposed conductor. Nevertheless, we conclude from the MT data that the main occurrence of magma and melt in the Lake District of the MER is prevalently at lower crustal depths along the SDFZ west of the main rift axis and not at the seismically active faults along the central rift axis.

Our study shows the importance of making a regional resistivity study in the rift valley. Understanding the nature of the proposed deep conductor underneath the SDFZ could play a major role in getting a clear picture of the magmatic evolution and the thermal regime in the MER. Future studies should also concentrate on analysing the nature and origin of deep well fluids that might be responsible for the unrest. Understanding the driving source is important for realistic evaluations of volcanic hazard in the MER.

We thank M. Becken, H. Brasse, E. Schill and Y. Abdelfettah for providing MT stations for the survey. We sincerely thank M. Kendall and J. Biggs who inspired us to launch the project. Thanks to U. Kalberkamp, H. Eysteinsson, T. Finnsson and K. Whaler for sharing their experiences with field work in Ethiopia. For their contribution to the successful accomplishment of the survey we sincerely thank E. Lewi, F. Gebre and B. Gebresilassie of the Institute for Geophysics, Space Science and Astronomy, Addis Ababa University, S. Kebede, B. Woldesmayet and Y. Lemma from the Geological Survey of Ethiopia and G. Andarge, N. Abera and the always helpful staff from the Aluto-Langano geothermal power plant. We thank the reviewers S. Hautot and G. Muñoz for their valuable comments which helped to improve our manuscript. Our research was supported by ETH grant No. ETH-3010-3 and the ETH scholarship for doctoral students.

APPENDIX A

The following Figs A1 to A6 show the selected apparent resistivities and phases of the observed and predicted responses as recovered by the inversion presented in Section 10.

Figure A1.

Observed and predicted responses at site R1 to M8.

Figure A2.

Observed and predicted responses at site M9 to M18.

Figure A3.

Observed and predicted responses at site M19 to M27.

Figure A4.

Observed and predicted responses at site M28 to M36.

Figure A5.

Observed and predicted responses at site M37 to M45.

Figure A6.

Observed and predicted responses at site M46 (reoccupation of site M1).