Abstract

The Okavango Delta of northern Botswana is one of the world's largest inland deltas or megafans. To obtain information on the character of sediments and basement depths, audiomagnetotelluric (AMT), controlled-source audiomagnetotelluric (CSAMT) and central-loop transient electromagnetic (TEM) data were collected on the largest island within the delta. The data were inverted individually and jointly for 1-D models of electric resistivity. Distortion effects in the AMT and CSAMT data were accounted for by including galvanic distortion tensors as free parameters in the inversions. By employing Marquardt–Levenberg inversion, we found that a 3-layer model comprising a resistive layer overlying sequentially a conductive layer and a deeper resistive layer was sufficient to explain all of the electromagnetic data. However, the top of the basal resistive layer from electromagnetic-only inversions was much shallower than the well-determined basement depth observed in high-quality seismic reflection images and seismic refraction velocity tomograms. To resolve this discrepancy, we jointly inverted the electromagnetic data for 4-layer models by including seismic depths to an interface between sedimentary units and to basement as explicit a priori constraints. We have also estimated the interconnected porosities, clay contents and pore-fluid resistivities of the sedimentary units from their electrical resistivities and seismic P-wave velocities using appropriate petrophysical models. In the interpretation of our preferred model, a shallow ∼40 m thick freshwater sandy aquifer with 85–100 Ωm resistivity, 10–32 per cent interconnected porosity and <13 per cent clay content overlies a 105–115 m thick conductive sequence of clay and intercalated salt-water-saturated sands with 15–20 Ωm total resistivity, 1−27 per cent interconnected porosity and 15–60 per cent clay content. A third ∼60 m thick sandy layer with 40–50 Ωm resistivity, 10–33 per cent interconnected porosity and <15 per cent clay content is underlain by the basement with 3200–4000 Ωm total resistivity. According to an interpretation of helicopter TEM data that cover the entire Okavango Delta and borehole logs, the second and third layers may represent lacustrine sediments from Paleo Lake Makgadikgadi and a moderately resistive freshwater aquifer comprising sediments of the recently proposed Paleo Okavango Megafan, respectively.

INTRODUCTION

Geological and hydrological setting

The Okavango Delta occupies a trough within the Kalahari Basin of northern Botswana close to the border with Namibia (Fig. 1). It is one of Earth's largest inland deltas or megafans consisting of a multitude of channels, islands and both perennial and seasonal swamps that make it a unique ecosystem.

Figure 1.

Map of northern Botswana showing the Okavango Delta and investigation sites. White outline: area of Okavango Delta covered by regional helicopter transient electromagnetic (HTEM) measurements; HR1 and HR2: areas of high-resolution HTEM measurements; Wheat Fields (WF): site of present investigation on Chief's Island; HR2, Jao and Golf Course (GC): previously investigated sites using seismic, ERT and TEM methods (Meier et al.2014; Reiser et al.2014; Podgorski et al.2015); blue lines: major rivers; thick solid lines: prominent faults where G-Gumare, K-Kunyere, T-Thamalakane (approximate positions after Shaw & Thomas 1992; Gumbricht et al.2001); bh: borehole. Thata Island is located within HR1. Landsat 8 satellite image courtesy of the U.S. Geological Survey.

Figure 1.

Map of northern Botswana showing the Okavango Delta and investigation sites. White outline: area of Okavango Delta covered by regional helicopter transient electromagnetic (HTEM) measurements; HR1 and HR2: areas of high-resolution HTEM measurements; Wheat Fields (WF): site of present investigation on Chief's Island; HR2, Jao and Golf Course (GC): previously investigated sites using seismic, ERT and TEM methods (Meier et al.2014; Reiser et al.2014; Podgorski et al.2015); blue lines: major rivers; thick solid lines: prominent faults where G-Gumare, K-Kunyere, T-Thamalakane (approximate positions after Shaw & Thomas 1992; Gumbricht et al.2001); bh: borehole. Thata Island is located within HR1. Landsat 8 satellite image courtesy of the U.S. Geological Survey.

About 6 km3 of water enter the Okavango Delta through direct rainfall during the wet season in January and February (McCarthy 2006), and the Okavango River originating in the Angolan highlands contributes roughly 10 km3 of water per year. Most of the riverine waters enter the Okavango Delta during annual floods that occur, with some delay, after the rainfalls in the Angolan highlands. The peak of the flood affects the Panhandle region (Fig. 1) in late April to early May and then propagates slowly through the delta reaching its southeastern margin in late July to late August. Because of the hot semi-arid climate, about 95 per cent of the water entering the delta is lost through evapotranspiration; only a small volume of water percolates into the aquifer systems or leaves the delta by surface flow (Gieske 1996; McCarthy 2006; Milzow et al.2009). As a result of evapotranspiration, concentrations of dissolved solutes in groundwater are expected to increase with time. Different mechanisms have been proposed during the last two decades to explain how the shallow aquifers can maintain freshwater by effectively removing minerals from the shallow subsurface. First, uptake of soluble ions in peat may remove solutes from swamp and channel waters (McCarthy et al.1989). Second, both seasonal flooding and rainfall may lead brines to sink into deeper parts of the aquifers (McCarthy & Metcalfe 1990). Third, salt that accumulates beneath many deltaic islands as a result of evapotranspiration (McCarthy et al.1991) may be transported to greater depths by density-driven flow (so called ‘salt fingering’) when the brine density exceeds the density of the fluids in the underlying aquifer (Gieske 1996; Bauer et al.2006; Zimmermann et al.2006; Bauer-Gottwein et al.2007; Milzow et al.2009). At Thata Island within the HR1 area of Fig. 1, electrical resistance tomographic (ERT) and ground-based transient electromagnetic (TEM) methods have been successfully used to outline increased brine concentrations resulting from density-driven flow (Bauer et al.2006; Bauer-Gottwein et al.2010).

Most of the Okavango Delta lies between the northeast–southwest trending normal faults (Gumare and Thamalakane faults in Fig. 1) of the asymmetric Okavango Rift, which is commonly interpreted as a southwestern branch of the East-African Rift System (Reeves 1972; Scholz et al.1976; Modisi et al.2000; Gumbricht et al.2001; Kinabo et al.2008). The width of the Okavango Rift is approximately 150 km between the Gumare and Thamalakane faults (Fig. 1), and the downthrow is more than 300 m along its southeastern margin (fig. 5 in Podgorski et al.2013b).

Tectonic activity has impacted the channel morphology within the delta (McCarthy et al.1993; McCarthy 2006) and has led to repeated re-directions of major river systems throughout the Kalahari Basin (Thomas & Shaw 1991; McCarthy et al.1993; McCarthy 2006; Moore & Larkin 2001; Burrough et al.2009; Moore et al.2012). As a result of tectonic movements, hydrological conditions and sedimentation cycles within the Kalahari Basin are assumed to have varied between fluvial, lacustrine, megafan and aeolian (Podgorski et al.2013b).

An area of 66 000 km2 to the south and the east of the Okavango Delta is understood to have been covered by a mega-lake in the Pleistocene epoch (Burrough et al.2009). Thomas & Shaw (1991) have speculated that this Paleo Lake Makgadikgadi (PLM) may even have covered an area of 120 000 km2 that included the entire present Okavango Delta. The age and evolution of the mega-lake are not well established (Ringrose et al.2008; Burrough et al.2009; Moore et al.2012). According to the age of pollens from beneath the Makgadikgadi Pans, formation of PLM may have commenced in the Miocene (Moore et al.2012).

Previous geophysical studies

A helicopter transient electromagnetic (HTEM) survey was commissioned in 2007 by the Department of Geological Survey of Botswana to study the division between freshwater and salt-water aquifers and the depth to basement (Figs 1 and 2a). The HTEM survey covered the whole Okavango Delta with 2 km line spacing and two smaller areas with 50 m line spacing, one at its western edge (HR2) and the other in the northeast (HR1). The editing, processing and inversion methodologies applied to the HTEM data have been described in detail by Podgorski et al. (2013a). After processing, station spacing along lines was 25 m for all data sets.

Figure 2.

(a) Okavango Delta and Wheat Fields (WF) investigation site. (b–f) Resistivity model derived from a laterally constrained inversion of the regional HTEM data depicted as slices at different depths below the surface. Justification for the paleo megafan is briefly provided in the text and fully explained in Podgorski et al. (2013b). Figure modified from Meier et al. (2014).

Figure 2.

(a) Okavango Delta and Wheat Fields (WF) investigation site. (b–f) Resistivity model derived from a laterally constrained inversion of the regional HTEM data depicted as slices at different depths below the surface. Justification for the paleo megafan is briefly provided in the text and fully explained in Podgorski et al. (2013b). Figure modified from Meier et al. (2014).

Although nearly all of the HTEM data can be explained by a simple 3-layer resistive–conductive–resistive model, their joint interpretation with borehole logs and seismic data necessitates a 4-layer model beneath a large region of the delta (Podgorski et al.2013b); essentially, the basal layer of the 3-layer model requires two components. The surface resistive layer (Fig. 2b) is interpreted as a freshwater aquifer with limited salinity. The underlying conductive layer (Figs 2c–e) is explained as lacustrine sediments consisting of clay units intercalated with sandy units containing saline pore water. In a regional context, Podgorski et al. (2013b) propose that this layer was once part of PLM, suggesting that PLM may have covered an area of more than 90 000 km2, thus providing support for Thomas & Shaw's (1991) earlier speculation. Starting at a depth of 100 m, a more resistive basal structure is observed under the Panhandle and at greater depth also beneath the western, central and northeastern parts of the delta. Podgorski et al. (2013b) interpret the western part of this resistive structure to be basement, whereas it is interpreted to be a freshwater aquifer (i.e. upper component of the basal resistive layer) underlain by basement (i.e. lower component of the basal resistive layer) in the central and northeastern parts of the delta. In the HTEM models, the poorly defined resistivity and thickness of the proposed lower freshwater aquifer may range from 40 to 300 Ωm and from 20 to 50 m, respectively. Its semi-circular outline in Fig. 2(f) suggests that it may represent a paleo megafan that is referred to as the Paleo Okavango Megafan (POM) by Podgorski et al. (2013b).

Since the resistivity of the presumed POM probably lies between the resistivities of the overlying lacustrine PLM layer and the underlying basement and it is thin compared to its burial depth, we face a classic suppression or equivalence issue in the inversion of the electromagnetic data (e.g. Raiche et al.1985; Albouy et al.2001). To resolve this issue and, hence, to make sound judgements on the suggested paleo megafan, seismic data can be utilized to distinguish basement with high seismic P-wave velocities (vp ∼ 4500 m s−1) from the overlying sedimentary units with low seismic P-wave velocities (vp ∼ 1800 m s−1). On the western side of the delta (at HR2 in Fig. 1), a comparison of seismic reflection images, tomographic seismic refraction models and borehole information (Reiser et al.2014) with inversions of ERT and TEM data (Meier et al.2014) demonstrates that the deep resistive zone observed in the HTEM models is basement at that location. Underneath the north-central and central areas of the delta at Jao and Golf Course in Fig. 1, comparisons of seismic and ERT/TEM models (Reiser et al.2014; Meier et al.2014; Podgorski et al.2015) tend to support the interpretation of Podgorski et al. (2013b) that at these locations the resistive zone is the POM freshwater aquifer lying on top of basement.

Further tests of the POM hypothesis required

The results of the geophysical surveys in the north-central and central areas of the delta are somewhat inconclusive for several reasons: (i) absence of the characteristic shallow resistive layer and failure to demonstrate compatibility of the resistivity and seismic models at the Jao site, (ii) inadequate penetration of electromagnetic signals at the Golf Course site and (iii) an inability to determine confidently the resistivities of units underlying the conductive layer at both sites.

To test further the POM hypothesis, seismic, ERT and diverse types of electromagnetic measurement were made within an area called the Wheat Fields on Chief's Island in the north-central part of the Okavango Delta (Fig. 1). The electromagnetic methods included audiomagnetotelluric (AMT), controlled-source audiomagnetotelluric (CSAMT) and central-loop TEM. Compared to the 1-D models derived from the HTEM data, it was anticipated that the combination of electromagnetic data generated by both inductively and galvanically coupled sources would yield improved constraints on the resistivities of the more resistive structures. Furthermore, joint interpretation of the new electromagnetic models, seismic reflection images and tomographic seismic refraction models was expected to provide a more convincing test of the POM hypothesis. The processing scheme applied to the seismic reflection data and the tomographic method used to invert the seismic refraction data have been described in Reiser et al. (2014) and the results presented in Podgorski et al. (2015). Since the ERT data are heavily affected by coupling problems, they contain information only for the upper 25 m of the sedimentary sequence (Podgorski et al.2015) and, therefore, are not considered here. 1-D inversions of the TEM data by Podgorski et al. (2015) suggest compatibility of the TEM and seismic models, but the TEM data recorded at the Wheat Fields site lack sufficient depth penetration and rigorous assessment.

Outline

In the following, we briefly describe the data acquisition, the relevant theory of the electromagnetic methods and our inversion approaches. We introduce a new formulation for explicitly including constraints offered by seismic data and/or borehole information in joint inversions of multiple types of electromagnetic data. The relevant data processing and editing steps are then summarized and the single- and joint-inversion results are described in some detail. Using the seismic depths to two layer interfaces as constraints on the joint inversion of the electromagnetic data, we evaluate the compatibility of the electromagnetic data with the seismic reflection images and tomographic seismic refraction models by Podgorski et al. (2015) and examine the validity of the POM hypothesis. We relate the electrical resistivities and seismic P-wave velocities of the sedimentary layers to consistent ranges of porosity, clay content and pore-fluid resistivity using suitable petrophysical models for mixtures of sand, clay and pore fluid. Finally, we provide geological and petrophysical interpretations of the individual layers.

DATA ACQUISITION

The locations of the AMT and CSAMT receiver stations, grounded CSAMT transmitter cables, central-loop TEM soundings and the HTEM line within the Wheat Fields site are displayed in Fig. 3. The most notable surface feature of geoelectric significance is the river north of the CSAMT transmitter electrodes.

Figure 3.

Wheat Fields site showing recording lines A and B, locations of AMT and CSAMT receivers (Rx, tilted black rectangles), long grounded cable sources of CSAMT transmitter (Tx, inverted triangles at grounding points), central-loop TEM sounding sites (red squares for transmitter loop), seismic lines (dark blue lines), the nearest HTEM profile (orange line) and the location of a single HTEM sounding site (orange cross) for which a 1-D resistivity model is reproduced in many subsequent figures. AMT and CSAMT data were recorded and inverted in the co-ordinate system depicted by the x and y axes. Note the system of river channels located to the north and east of the CSAMT transmitter. Satellite image provided by the European Space Agency.

Figure 3.

Wheat Fields site showing recording lines A and B, locations of AMT and CSAMT receivers (Rx, tilted black rectangles), long grounded cable sources of CSAMT transmitter (Tx, inverted triangles at grounding points), central-loop TEM sounding sites (red squares for transmitter loop), seismic lines (dark blue lines), the nearest HTEM profile (orange line) and the location of a single HTEM sounding site (orange cross) for which a 1-D resistivity model is reproduced in many subsequent figures. AMT and CSAMT data were recorded and inverted in the co-ordinate system depicted by the x and y axes. Note the system of river channels located to the north and east of the CSAMT transmitter. Satellite image provided by the European Space Agency.

AMT and CSAMT horizontal electric field and horizontal and vertical magnetic field data were recorded at six stations along profile A in the frequency range 1 Hz–10 kHz using a Phoenix V8 receiver. The electric and magnetic field sensors were arranged in the x- and y-directions oriented along and perpendicular to the profile length (A in Fig. 3). For the CSAMT survey, two almost perpendicularly oriented long grounded cables were used as sources to generate two sets of vectorial electromagnetic fields. For both source cables, the amplitudes of the injected currents were about 3 A, varying somewhat with frequency. Given the frequency range, the expected depth of investigation for the AMT and CSAMT measurements ranged from a few tens of metres to about 1000 m.

TEM data were acquired using a Geonics PROTEM-47D system with a 100 m×100 m transmitter loop and a receiver coil for vertical magnetic field transients centred within the transmitter loop. Transient decays were recorded at repetition rates of 237.5, 62.5 and 25 Hz with currents of 1, 3 and 3.3 A, respectively. For these repetition rates, a total time range of 6.8 μs to 6.9 ms yields expected investigation depths between a couple of tens of metres and about 150 m. In total, eight TEM soundings were made along profiles A and B with a station spacing of ∼240 m.

THEORY

AMT method

The sources of AMT fields are thunderstorms located mostly in the tropics. Electromagnetic fields generated by thunderstorms propagate within the lossy waveguide defined by the Earth's surface and the ionosphere (Vozoff 1991). A small part of the energy propagates into the subsurface and is dissipated as heat through the electromagnetic skin effect. A characteristic quantity is the skin depth:  

(1)
\begin{equation} \delta _{{\rm AMT}}=\sqrt{\frac{2}{\omega \mu _0\hat{\sigma }}}, \end{equation}
where the amplitudes of the AMT fields are diminished to 1/e of their values at the surface, ω = 2πf is the angular frequency of the electromagnetic field, μ0 is the magnetic permeability of free space and $$\hat{\sigma }$$ is an average electrical conductivity (Spies 1989). Skin effect depth is only an approximation however for non-uniform subsurfaces, in particular the electric and magnetic fields attenuate differently for layered strata (e.g. Jones 2006). At sufficiently large distances from thunderstorms, the electromagnetic fields can be considered to be plane waves. At the Earth's surface, the horizontal electric fields Ex and Ey (in V m−1) and the horizontal and vertical magnetic fluxes Bx, By and Bz (in nT) are measured. Subsequently, the magnetic fluxes are converted to horizontal and vertical magnetic fields Hx, Hy and Hz (in A m−1). To make the obtained quantities independent of the unknown amplitude and phase factors that they have in common, the electromagnetic field components are connected through transfer functions in the frequency domain (i.e. through the impedance tensor Z and the vertical magnetic transfer function or tipper tensor T).

The impedance tensor Z relates the horizontal magnetic to the horizontal electric fields as (Berdichevsky & Dmitriev 2008):  

(2)
\begin{equation} \left[\begin{array}{c}E_x\\ E_y \end{array}\right] = \left[\begin{array}{c@{\quad}c}Z_{xx} & Z_{xy}\\ Z_{yx} & Z_{yy} \end{array}\right] \left[\begin{array}{c}H_x\\ H_y \end{array}\right], \end{equation}
where  
(3)
\begin{equation} \bf {Z} = \left[\begin{array}{c@{\quad}c}Z_{xx} & Z_{xy}\\ Z_{yx} & Z_{yy} \end{array}\right]. \end{equation}
On the surface of a 1-D Earth, Zxy = −Zyx and Zxx = Zyy = 0 V A−1. To provide a more obvious link to the resistivity distribution of the subsurface, impedance tensor elements Zij can be transformed to apparent resistivities ρa, ij and phases φij through:  
(4)
\begin{eqnarray} \rho _{a,ij} & = & \frac{1}{\omega \mu _0} \left|Z_{ij}\right|^2, \end{eqnarray}
 
(5)
\begin{eqnarray} \varphi _{ij} & = & \arg {Z_{ij}}, \end{eqnarray}
which on the surface of a homogeneous half-space would equal the half-space resistivity and 45°, respectively.

The vertical magnetic transfer function tensor T = [TxTy] relates the vertical and horizontal components of the magnetic field as:  

(6)
\begin{equation} H_z = T_x \cdot H_x + T_y \cdot H_y. \end{equation}
On the surface of a 1-D model, Tx = Ty = 0, because electromagnetic fields entering the Earth at AMT frequencies are refracted towards the normal, independent of the angle of incidence.

In a 1-D subsurface all induced currents flow horizontally, such that thin resistors are hard to impossible to detect, because their large skin depths yield only small amplitude and phase changes in Z. For perfect data at all frequencies there exists a uniqueness theorem in MT for a 1-D Earth such that only one model will precisely fit the data (Bailey 1970). Non-uniqueness in the MT method in 1-D is thus due to data insufficiency, data inadequacy and data error, and that non-uniqueness particularly is reflected in poor definition of the actual resistivity of resistive layers. Typically a minimum value can be set, but the layer has to be of reasonable extent to determine the true resistivity (e.g. Jones 1999).

An inhomogeneity in the vicinity of a receiver results in both electric and magnetic field anomalies. If an anomalous body has a depth of burial and dimensions much smaller than the skin depths inside and outside the anomaly, the deviations of the Z and T tensors from their background values Z0 and T0 are referred to as galvanic distortions. They can be expressed through real-valued frequency-independent distortion parameters P and Q (Wannamaker et al.1984; Garcia et al.2003; Jones 2012), where P and Q are given by:  

(7)
\begin{eqnarray} \bf {Z} & = & \left(\bf {I} + \bf {P}_h\right) \bf {Z}_0 \left(\bf {I}+\bf {Q}_h\bf {Z}_0\right)^{-1}, \end{eqnarray}
 
(8)
\begin{eqnarray} \bf {T} & = & \left[\bf {T}_0+\bf {Q}_v\bf {Z}_0\right] \left(\bf {I} + \bf {Q}_h\bf {Z}_0\right)^{-1}, \end{eqnarray}
where I is the identity matrix,  
(9)
\begin{eqnarray} \bf {P}_h & = & \left( \begin{array}{c@{\quad}c}P_{xx} & P_{xy}\\ P_{yx} & P_{yy} \end{array} \right) \end{eqnarray}
consists of the distortion parameters of the horizontal electric field,  
(10)
\begin{eqnarray} \bf {Q}_h & = & \left( \begin{array}{c@{\quad}c}Q_{xx} & Q_{xy}\\ Q_{yx} & Q_{yy} \end{array} \right) \end{eqnarray}
is the distortion tensor of the horizontal magnetic field components and  
(11)
\begin{eqnarray} \bf {Q}_{v} & = & \left( \begin{array}{c@{\quad}c}Q_{zx} & Q_{zy} \end{array} \right) \end{eqnarray}
is the distortion tensor of the vertical magnetic field. Since the distance to the source is large and relatively low frequencies are employed in AMT surveys, galvanic distortion of the magnetic field components in AMT measurements is usually negligible (Wannamaker et al.1984; Qian 1994). Kalscheuer et al. (2012) provide relevant details on the computation of AMT forward responses and sensitivities.

CSAMT method

In the CSAMT method (Zonge & Hughes 1991), active sources (e.g. long grounded cables or horizontal magnetic dipoles) are used to generate electromagnetic signals. To obtain full Z and T tensors, pairs of source cables or dipoles, preferably oriented perpendicular to each other are employed in the field (Li & Pedersen 1991). If the source-receiver separation r is much smaller than (near-field zone) or comparable to (transition zone) the plane-wave skin depth δAMT in eq. (1), the recorded Z and T tensors depend on r, on conductivity σ and in the transition zone additionally on frequency. Long grounded cable sources rely on both inductive coupling (through temporal variations of the source current) and galvanic coupling (through injection of the source current at the grounding points) to the subsurface. Hence, the currents in a 1-D Earth have both horizontal and vertical components in the near-field and transition zones. The vertical current component generates surface charges on layer interfaces that are proportional to the conductivity contrast. In turn, these surface charges modify the electric field and, hence, may allow for better detectability of resistive layers than is possible with purely inductively coupled methods (such as the AMT method) alone. Previously published 1-D inversion schemes for CSAMT data include Routh & Oldenburg (1999), Garcia et al. (2003) and Kalscheuer et al. (2012), the latter two offering the advantage that the layer and distortion parameters can be inverted for simultaneously. Kalscheuer et al.'s (2012) scheme offers the extra advantage of modelling fields from long grounded cable sources that fully account for the volumes and positions of the sources. Further details on how we compute CSAMT forward responses and sensitivities are provided by Kalscheuer et al. (2012).

Compared to AMT data, distortions of CSAMT magnetic fields due to small inhomogeneities in the vicinity of a receiver may be significant, being strongest for short transmitter-receiver separations (Qian 1994). In contrast, small localized 3-D anomalies close to the transmitter usually lead to changes of the effective source moment that may hardly affect the Z and T tensors, even though the distortion on the individual electric and magnetic fields may be substantial (Qian 1994; Kalscheuer et al.2012). But large extended anomalous bodies close to the transmitter may change the effective source moment and redirect current flow over large volumes, thus markedly distorting the Z and T tensors. Galvanic distortion effects due to sizable inhomogeneities in the vicinity of the transmitter are referred to as ‘source overprint’ (Zonge & Hughes 1991).

TEM method

Through the relatively sudden extinction of a constant current I in a transmitter loop of area A, eddy currents are induced in the subsurface (Nabighian & Macnae 1991). The transient decay of the time-derivative of the vertical magnetic flux component ∂bz/∂t is recorded as a voltage induced in a receiver coil over a certain time range after current shut-off. In a 1-D layered Earth, the induced currents flow purely horizontally, diffusing with time to greater depths and larger horizontal distances. The diffusion depth  

(12)
\begin{equation} \delta _{{\rm TEM}}=\sqrt{\frac{2t}{\hat{\sigma }\mu _0}} \end{equation}
is the depth at which the induced electric field attains its maximum at a given time t after current extinction (Spies 1989; Nabighian & Macnae 1991). In central-loop measurements, the maximum depth of investigation depends on the background noise level, the source moment I · A and the diffusion depth.

We summarize our approach for computing TEM forward responses and sensitivities with due regard for instrumental effects in Appendix A.

Inverse modelling

An appropriate procedure is required to model AMT, CSAMT and TEM data with a minimum number of layers in 1-D inversions. The electromagnetic data are first inverted for models consisting of numerous layers of constant thickness with resistivities varying gradually with depth and (possibly tensorial) galvanic distortion parameters in preliminary inversions. An adequate model regularization may consist of up to three terms: (i) smoothness constraints on the layered model (Constable et al.1987), (ii) minimum solution length constraints on electric and magnetic galvanic distortion parameters for AMT and CSAMT data (Kalscheuer et al.2012) and (iii) additional Marquardt–Levenberg damping on all model parameters to enforce convergence (Lines & Treitel 1984). The layered model is typically divided into 30–50 layers with the layer thicknesses increasing by factors of 1.15 to 1.3 from one layer to the next deeper layer. From the resulting smoothness-constrained layered model, resistivities and thicknesses of simplified models with reduced numbers of layers may then be deduced (e.g. Kalscheuer et al.2007). Finally, these layer parameters and the distortion tensors from the smoothness-constrained inversion of the AMT and CSAMT data may be used as a start model in an inversion where Marquardt–Levenberg damping is applied to layer resistivities and thicknesses and distortion parameters are regularized with minimum solution length constraints. Details on the choice of Lagrange multipliers for the model regularization and on the computation of root-mean-square (RMS) misfits of data and model responses normalized by absolute data uncertainties are presented in Appendix B.

In inversions of electromagnetic data, suppression and equivalence problems (Raiche et al.1985; Albouy et al.2001) can be significantly mitigated by including depths to lithological interfaces as determined from borehole logs or seismic sections as a priori information (e.g. Auken & Christiansen 2004). For inversion schemes that employ logarithms of resistivities ρ and thicknesses t of NL layers as primary model parameters, such as the previously described hybrid Marquardt–Levenberg and minimum solution length inversion, a priori information on depths to layer interfaces can conveniently be implemented as an additional term in the cost function. To demand that the depths z of a number of Nic ≤ NL − 1 layer interfaces be close to Nica priori or reference values zr, the additional term takes the form:  

(13)
\begin{eqnarray} \gamma Q_m^z & = & \gamma \left(\bf {z}_r-\bf {z}\right)^T\left(\bf {W}^{z}_{m}\right)^T\bf {W}^{z}_{m}\left(\bf {z}_r-\bf {z}\right)\nonumber \\ & = & \gamma \sum _{i=1}^{N_{ic}}\left(\frac{z_{r,i}-\sum _{j=1}^{l(i)-1}e^{\log t_j}}{\sigma _{z_{r,i}}}\right)^2, \end{eqnarray}
where we use summation indices l(i) to retrieve the depths of the correct layer interfaces. The diagonal matrix $$\bf {W}^{z}_{m}$$ contains the reciprocals of the uncertainties (standard deviations) $$\sigma _{z_{r,i}}$$ of the a priori depths, and the factor γ is a Lagrange multiplier that expresses the importance of the depth constraints relative to the data and other model constraints. Further details on the implementation of the constraints for depths to layer interfaces are given in Appendix B.

We assess model uniqueness and stability using the truncated singular value decomposition (TSVD) scheme of Kalscheuer & Pedersen (2007). In contrast to linearized approaches, this model error and resolution analysis scheme has the advantage that the nonlinearity of the inversion problem is partly accommodated via the computation of model errors with nonlinear most-squares inversion (Jackson 1976; Meju & Hutton 1992; Meju 1994). Typically, most-squares inversion yields more reliable error estimates than linearized schemes, which often underestimate true model uncertainty (Kalscheuer & Pedersen 2007; Kalscheuer et al.2010). Compared to Kalscheuer & Pedersen (2007), we consider models with very limited numbers of parameters (layer and distortion parameters). Hence, we include the entries of the Jacobian matrix of all model parameters in the TSVD analysis of a single parameter and employ the same truncation level in the analyses of all model parameters.

PROCESSING AND INVERSION

Processing

The AMT time series were processed using a standard single-station processing algorithm that included (i) division of the time series into sections, (ii) tapering and Fourier transformation of the individual sections and (iii) computation of transfer functions at individual frequencies using cross-spectra and auto-spectra that involve either the local H or E field (Simpson & Bahr 2005; Chave 2012). Standard deviations of all AMT impedances were adjusted to reflect an uncertainty floor of 2 per cent. Excessively noisy data were deleted in the AMT dead band between 800 and 3000 Hz (Vozoff 1991; Garcia & Jones 2002; Sternberg 2010); it was necessary to remove all data within the dead band recorded at station 9 (Fig. 3). Due to the relatively short acquisition time of roughly 1 hr per station, the lowest usable frequency was in the 3–5 Hz range.

For the CSAMT data at each frequency, we estimated the Z and T tensors by combining the electric and magnetic fields generated by the individual source cables as described in Li & Pedersen (1991). Standard deviations of the Z and T tensors were computed using the covariance approach of Bastani (2001); the off-diagonal elements of the co-variance matrices were disregarded. We assumed uncertainty floors of 3 per cent for all Z tensor elements and 0.03 for the T tensor elements. Because of the large distance to the source, the CSAMT data for station 10 were relatively noisy.

To remove instrumental effects from the AMT and CSAMT transfer functions, the magnetic field components were divided by the system responses of the induction coils in the frequency domain. Since the system responses of the receiver filters affected the electric and magnetic field measurements to the same extent, the estimated transfer functions did not contain the system response effects of the receiver filters.

For the TEM data, outliers were removed manually considering both apparent resistivity and $$\frac{\partial {}b_z}{\partial {}t}$$ transients. Subsequently, individual records for the same station and repetition rate were stacked. In stacking, uncertainty estimates were determined using a noise model similar to that of Auken et al. (2008); we combined a uniform noise contribution of 5 per cent related to possible instrumental noise with an absolute background noise contribution estimated from the field data. To account for instrumental effects, the forward responses and sensitivities were convolved in the time-domain with the system responses for the transmitter current form (Geonics Ltd. 1998), the receiver coil and the amplifiers in the receiver (Geonics Ltd. 2011) during inversion.

3-Layer versus 4-layer resistivity models

From our previous experiments in the Okavango Delta (Podgorski et al.2013b; Meier et al.2014), it is known that many of the ground and airborne electromagnetic data sets can be explained with a simple 3-layer model, in which a thick layer of conductive material (i.e. the lacustrine sediments of the PLM) is sandwiched between two resistive units (i.e. the upper shallow freshwater aquifer and basement). Podgorski et al. (2013b) interpret an additional comparatively thin layer of moderate to high resistivity (i.e. the freshwater aquifer of the POM) to lie between the electrically conductive layer and the underlying resistive basement in the north-central part of the Okavango Delta (including the Wheat Fields area). However, this thin layer is difficult to detect using electromagnetic methods alone, because it represents a typical case in which suppression or equivalence occurs in the modelling and inversion of electromagnetic data (e.g. Albouy et al.2001), seemingly permitting interpretations with just three layers.

After the initial smoothness-constrained inversions, we applied Marquardt–Levenberg inversions for 3-layer models to the individual data sets and to the joint AMT, CSAMT and TEM data sets at each station (not presented here). In these models, basement resistivities were usually in the expected range of 3000–4000 Ωm, whereas the basement depths of 150–180 m were much shallower than the 205 ± 10 m depth in the seismic models of Podgorski et al. (2015). Although we briefly discuss 3-layer models in the section on joint inversion results (e.g. Table 1), because it was not possible to derive 3-layer resistivity models that were compatible with the seismic reflection images and tomographic seismic refraction models, we focus in the following on Marquardt–Levenberg inversion models with four layers that reproduce basement depths in the correct range. For the joint data sets, we additionally present 4-layer models for which the seismic depths to the top of the POM and basement were used as constraints.

Table 1.

Overview of CSAMT/AMT and TEM data sets used for the joint inversions and resulting data fits (RMS misfits) for smooth layered, 3-layer and 4-layer resistivity models. Station numbers are shown in Fig. 3. For one set of 3-layer models, seismic constraints on depth to basement were included during the inversion; these models have distinctively higher RMS misfits than the other models indicating incompatibility of the seismic and 3-layer resistivity models. In contrast, the 4-layer resistivity models that were constrained by seismic depths to the POM and basement have data fits as good as the best electromagnetic inversion results without seismic constraints.

Data sets with RMS misfits for different inversion approaches 
station labels with seismic depth constraints in brackets 
AMT/ TEM Smooth 3-Layer 3-Layer ML model 4-Layer 4-Layer ML model 
CSAMT  R1 model ML model (seismic basement) ML model (seismic POM and basement) 
1.0 1.0 1.3 1.0 1.0 
1.0 1.0 1.5 1.0 1.0 
12 1.2 1.3 1.6 1.3 1.3 
1.1 1.2 1.7 1.2 1.2 
1.9 2.0 2.4 2.0 2.0 
Data sets with RMS misfits for different inversion approaches 
station labels with seismic depth constraints in brackets 
AMT/ TEM Smooth 3-Layer 3-Layer ML model 4-Layer 4-Layer ML model 
CSAMT  R1 model ML model (seismic basement) ML model (seismic POM and basement) 
1.0 1.0 1.3 1.0 1.0 
1.0 1.0 1.5 1.0 1.0 
12 1.2 1.3 1.6 1.3 1.3 
1.1 1.2 1.7 1.2 1.2 
1.9 2.0 2.4 2.0 2.0 

For each suite of inversions, we first show models and data for a typical single observation station (e.g. station 3 or station 4) and then show compilations of the models for all stations.

Inversion of AMT data

At all sites, the AMT transfer functions indicate 1-D conditions with very small diagonal values and nearly identical off-diagonal values in the Z tensor and values below 0.1 in the T tensor. Therefore, we limit further interpretations of the AMT data to the impedance tensor elements Zxy and Zyx that only involve sets of scalar distortion parameters Pxx, Pyy, Qxy and Qyx. Although Qxx and Qyy may deviate from zero when inverting the off-diagonal Z tensor elements, there is very little sensitivity to these distortion parameters in our AMT data.

In Fig. 4, we present the layered models and distortion parameters P and Q together with the apparent resistivities and phases of Zxy and Zyx for station 3. The R1 model with 45 layers was derived using first-order smoothness constraints for layer resistivities. The overall data fit at a normalized RMS misfit of 1.05 is very good. A similarly good data fit (RMS misfit of 1.06) is obtained by adjusting the resistivities and thicknesses of 4 layers using Marquardt–Levenberg regularization (model ML). The R1 and 4-layer ML models have similarly small distortion parameters and are in good agreement with the HTEM model of Podgorski et al. (2013b).

Figure 4.

(a) 1-D inversion models R1 and ML and distortion parameters P and Q for AMT station 3 (cf. Fig. 3) and (b) data fit of the response of model R1 to the off-diagonal elements of the AMT impedance tensor. Inversion models were computed using first-order smoothness constraints for R1 and Marquardt–Levenberg regularization for ML for layer parameters and a minimum solution length criterion for galvanic distortion parameters. Since only the off-diagonal impedance tensor elements were inverted, only four distortion parameters were involved.

Figure 4.

(a) 1-D inversion models R1 and ML and distortion parameters P and Q for AMT station 3 (cf. Fig. 3) and (b) data fit of the response of model R1 to the off-diagonal elements of the AMT impedance tensor. Inversion models were computed using first-order smoothness constraints for R1 and Marquardt–Levenberg regularization for ML for layer parameters and a minimum solution length criterion for galvanic distortion parameters. Since only the off-diagonal impedance tensor elements were inverted, only four distortion parameters were involved.

A compilation of all 4-layer ML models and distortion parameters is presented in Fig. 5. The depths to the resistive half-space vary between 185 and 230 m, which are in general agreement with the 205 ± 10 m seismic basement depth, which is practically horizontal beneath the entire investigation site. The thickness of the third layer varies between 30 and 110 m. The relatively large variability of the thicknesses and resistivities of the second and third layers is a result of equivalence in the AMT data. Moreover, there appears to be very little effect of galvanic distortion on the AMT data (Fig. 5b).

Figure 5.

1-D inversion models for Zxy and Zyx impedances of all AMT stations along line A (cf. Fig. 3) computed using Marquardt–Levenberg and minimum-solution length regularization for (a) layer and (b) galvanic distortion parameters, respectively. Although data at all stations can be explained by 3-layer models, constraints provided by seismic reflection images and refraction tomograms necessitate 4-layer models (e.g. as shown here).

Figure 5.

1-D inversion models for Zxy and Zyx impedances of all AMT stations along line A (cf. Fig. 3) computed using Marquardt–Levenberg and minimum-solution length regularization for (a) layer and (b) galvanic distortion parameters, respectively. Although data at all stations can be explained by 3-layer models, constraints provided by seismic reflection images and refraction tomograms necessitate 4-layer models (e.g. as shown here).

Inversion of CSAMT data

In the frequency range from 500–100 Hz, the CSAMT impedances of station 3 (Fig. 6b) show a rapid transition to the near-field zone with the imaginary parts of the Z tensor dropping to zero. Such behaviour is typical of geological settings where a conductive layer overlies resistive basement (Zonge & Hughes 1991; Boerner et al.1993; Wannamaker 1997). At all stations, the absolute values of the CSAMT responses Zyx and Zyy (i.e. impedance tensor elements related to the Ey field) are almost equal and about half that of Zxy (Fig. 6b). Two observations suggest, that energy leaked from Zyx into Zyy as a result of severe source overprint (Zonge & Hughes 1991). First, forward modelling tests based on the layered model derived from the HTEM survey (Podgorski et al.2013b) and the source-receiver geometry used in the field could not reproduce comparable responses. Second, there is no hint in the AMT impedances that there are significant 2-D or 3-D inductive effects or significant galvanic distortion below the receiver sites. A physical explanation for the suspected source overprint is substantial current channelling through the river adjacent to the northern transmitter electrode (Fig. 3).

Figure 6.

(a) 1-D inversion models R1 and ML and distortion parameters P and Q for CSAMT station 3 (cf. Fig. 3). (b) Fit of responses of model R1 to the CSAMT data. Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML) on layer parameters and a minimum solution length criterion for galvanic distortion parameters. Strong galvanic distortion Pyx and Pyy of Ey—and hence the Zyx and Zyy impedances—is evident at all receiver sites. This is likely an effective distortion of the transmitter current caused by current channelling through the river immediately north of the transmitter (cf. Fig. 3).

Figure 6.

(a) 1-D inversion models R1 and ML and distortion parameters P and Q for CSAMT station 3 (cf. Fig. 3). (b) Fit of responses of model R1 to the CSAMT data. Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML) on layer parameters and a minimum solution length criterion for galvanic distortion parameters. Strong galvanic distortion Pyx and Pyy of Ey—and hence the Zyx and Zyy impedances—is evident at all receiver sites. This is likely an effective distortion of the transmitter current caused by current channelling through the river immediately north of the transmitter (cf. Fig. 3).

In inverting the CSAMT data, the strong source overprint in Zyx and Zyy results in Pyx and Pyy values of about 0.5 and −0.6, respectively (Fig. 6a). As expected (Qian 1994; Kalscheuer et al.2010), both P and Q tensors are generally larger than their respective AMT counterparts. The fit of the observed and model predicted CSAMT data (RMS misfit of 0.98) in Fig. 6(b) is very good for all tensor components except for Zxx and Ty, which are rather low in amplitude and, hence, more affected by noise. The 4-layer ML model for station 3 has a data fit (RMS misfit of 0.95) comparable to that of model R1. There is an excellent correspondence between the CSAMT models of Fig. 6(a) and the independently derived AMT models of Fig. 4(a).

Fig. 7 displays all 4-layer ML models and distortion parameters derived from the CSAMT data acquired along line A. Because of the suppression problem, the resistivities and thicknesses of the third layer vary considerably from 20–150 Ωm and from 30–100 m, respectively. Similarly, the depth to the resistive half-space lies in the 185–270 m range. At stations 7 and 9, the second layer has resistivities and thicknesses of <1 Ωm and <10 m, respectively. By allowing the inversion to attain the minimum misfit at these two stations, the thicknesses of the second and third layer were decreased and increased, respectively, effectively transforming the suppression problem of the third layer (i.e. A-type equivalence) to a thin conductive layer problem of the second layer (i.e. H-type equivalence; Raiche et al.1985).

Figure 7.

1-D inversion models for data of all CSAMT stations along line A (cf. Fig. 3) computed using Marquardt–Levenberg and minimum-solution length regularization for (a) layer and (b) galvanic distortion parameters, respectively. As for the AMT models in Fig. 5, constraints provided by seismic data require 4-layer models. Strong similarities of distortion parameters P and Q along the line hint at significant source overprint.

Figure 7.

1-D inversion models for data of all CSAMT stations along line A (cf. Fig. 3) computed using Marquardt–Levenberg and minimum-solution length regularization for (a) layer and (b) galvanic distortion parameters, respectively. As for the AMT models in Fig. 5, constraints provided by seismic data require 4-layer models. Strong similarities of distortion parameters P and Q along the line hint at significant source overprint.

Inversion of TEM data

Since we did not record TEM data at station 3, we consider the closest TEM sounding at station 4 (Fig. 3). The normalized induced voltages in the last few time gates of the 237.5 Hz transient are slightly smaller than the corresponding voltages at 62.5 and 25 Hz as a result of the periodicity of the transmitter signal (Fig. 8a). At the lowest repetition rate of 25 Hz, excessively noisy values in the last 4 time gates were eliminated. Nevertheless, the accelerated decay rate of the TEM data in the last few time gates suggests the presence of a deeper resistive unit. Except for the half-space resistivity of the TEM R1 model, the distributions of resistivities in the R1 and 4-layer ML TEM models of Fig. 8(a) are similar to those of the AMT and CSAMT models of Figs 4(a) and 6(a). There are two reasons for the anomalously low half-space resistivity of the R1 model in Fig. 8(a). First, the combination of inductively coupled sources and magnetic field sensors employed in the central-loop TEM measurements does not constrain the resistivities of resistive units as well as electric field measurements of signals generated by galvanically coupled sources. Second, a comparison of models derived from first and second order smoothness constraints (not presented here) demonstrates that the maximum depth of investigation provided by the TEM data is 180–200 m suggesting that the basement is barely reached through the TEM soundings.

Figure 8.

(a) 1-D inversion models R1 and ML and (b) fit of the forward response of model R1 to the central-loop TEM ∂bz/∂t data for station 4 (cf. Fig. 3). Transients for all repetition rates (237.5 , 62.5 and 25 Hz) are normalized by the respective transmitter current strengths. Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML).

Figure 8.

(a) 1-D inversion models R1 and ML and (b) fit of the forward response of model R1 to the central-loop TEM ∂bz/∂t data for station 4 (cf. Fig. 3). Transients for all repetition rates (237.5 , 62.5 and 25 Hz) are normalized by the respective transmitter current strengths. Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML).

Fig. 9 summarizes the 4-layer models derived from Marquardt–Levenberg inversion of TEM data acquired along profiles A and B. As for the 4-layer ML models of the AMT and CSAMT data, the top of the resistive half-space is at roughly 180–230 m. The thickness of the overlying unit of intermediate resistivity varies from 30–100 m.

Figure 9.

1-D inversion models for the data of all central-loop TEM soundings computed using Marquardt–Levenberg inversion for (a) profile A and (b) profile B. As for the AMT models in Fig. 5, constraints provided by seismic data require 4-layer models.

Figure 9.

1-D inversion models for the data of all central-loop TEM soundings computed using Marquardt–Levenberg inversion for (a) profile A and (b) profile B. As for the AMT models in Fig. 5, constraints provided by seismic data require 4-layer models.

Joint inversion of AMT, CSAMT and TEM data

The AMT, CSAMT and TEM data were jointly inverted for simple layered models using the same procedures as employed for the individual data sets. AMT and CSAMT distortion parameters were computed separately during the inversions. Because the AMT/CSAMT and TEM recording stations were offset from each other (Fig. 3), we jointly inverted data from adjacent locations. A list of the five jointly inverted data sets, the inversion schemes and resulting RMS misfits is given in Table 1.

Typical R1 and ML joint inversion results for the AMT and CSAMT data of station 3 and the TEM data of station 4 are shown in Fig. 10. Since the individual inversions of the AMT, CSAMT and TEM data sets produce similar R1 and 4-layer ML models, it is not surprising that the jointly inverted models resemble the individual models of Figs 4(a)6(a) and 8(a). Note, that the total RMS misfits of 1.0 and the RMS misfits of the individual data sets (AMT: 1.1; CSAMT: 1.0; TEM: 0.7) for the jointly inverted R1 and 4-layer ML models for stations 3 and 4 in Table 1 are comparable to the RMS misfits for the individual inversion models (AMT: R1 and ML—1.1 and 1.1; CSAMT: R1 and ML—1.0 and 1.0; TEM: R1 and ML—0.6 and 0.8). In addition, the AMT and CSAMT distortion parameters for the jointly inverted models (those for the R1 model are displayed in Fig. 10) are very similar to those of the individual inversions, indicating that there are no serious trade-offs between any incompatibilities between the different data sets and the distortion parameters.

Figure 10.

(a) 1-D joint inversion models R1 and ML and distortion parameters P and Q and fit to (b) AMT, (c) CSAMT and (d) central-loop TEM data for stations 3/4 (cf. Fig. 3). Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML models for fourth and 25th iterations are shown) for layer parameters and a minimum solution length criterion for galvanic distortion parameters. Distortion parameters and data fits are for the R1 model.

Figure 10.

(a) 1-D joint inversion models R1 and ML and distortion parameters P and Q and fit to (b) AMT, (c) CSAMT and (d) central-loop TEM data for stations 3/4 (cf. Fig. 3). Inversion models were computed using first-order smoothness constraints (R1) and Marquardt–Levenberg damping (ML models for fourth and 25th iterations are shown) for layer parameters and a minimum solution length criterion for galvanic distortion parameters. Distortion parameters and data fits are for the R1 model.

Figure 10.

(Continued.)

Figure 10.

(Continued.)

The two 4-layer ML models for stations 3 and 4 in Fig. 10 were obtained during the same inversion run, one after the 4-th iteration with an RMS misfit of 1.01 and one after the 25th iteration with a slightly lower RMS misfit of 0.98. As a consequence of equivalence, the main differences between the two ML models are in the resistivities and thicknesses of the second and third layers. In progressing from the fourth to the 25th iteration, ρ3 becomes almost as small as ρ2 while t3 becomes substantially larger than t2.

All five jointly inverted 4-layer ML models along line A are presented in Fig. 11. As for the models displayed in Fig. 10, the variability of the second and third layer parameters are mainly due to equivalence rather than differences in actual layering. Along the profile, RMS misfits of the best fitting models mostly range from 1.0 to 1.3 (Table 1). The large RMS misfit of 2.0 for station 9 is caused by the poor quality of the AMT data at that location.

Figure 11.

1-D joint inversion models with four layers for the AMT, CSAMT and central-loop TEM data along line A (cf. Fig. 3) computed using a hybrid Marquardt–Levenberg and minimum-solution length inversion scheme. The layer model parameters, galvanic distortion parameters P and Q for AMT stations and galvanic distortion parameters P and Q for CSAMT stations are displayed in panels (a), (b) and (c), respectively. Seismic depths to the POM (picked as the uppermost seismic reflector from a package of reflectors between 145 and 170 m depth, see the text for more detail) and basement from Podgorski et al. (2015) are displayed for comparison. RMS misfits are listed in Table 1.

Figure 11.

1-D joint inversion models with four layers for the AMT, CSAMT and central-loop TEM data along line A (cf. Fig. 3) computed using a hybrid Marquardt–Levenberg and minimum-solution length inversion scheme. The layer model parameters, galvanic distortion parameters P and Q for AMT stations and galvanic distortion parameters P and Q for CSAMT stations are displayed in panels (a), (b) and (c), respectively. Seismic depths to the POM (picked as the uppermost seismic reflector from a package of reflectors between 145 and 170 m depth, see the text for more detail) and basement from Podgorski et al. (2015) are displayed for comparison. RMS misfits are listed in Table 1.

Joint inversion with seismic constraints

To help resolve the equivalence problem, we have included the seismic depths to the proposed POM and basement as additional constraints in the inversions. In Podgorski et al.'s (2015) seismic reflection images and tomographic refraction models (Fig. 12), the boundary between PLM and POM at 145 ± 10 m is defined by a continuous subhorizontal reflector that separates a zone of discontinuous weak sub-horizontal reflectors or no reflectors from an underlying zone of conspicuous continuous sub-horizontal reflectors. The lack of significant reflectors within the PLM deposits is consistent with seismic images of lacustrine units elsewhere (Büker et al.1998, 2000; Nitsche et al.2002, and references therein) and could be a consequence of fine layering of the lacustrine units that is below the resolution limits of the seismic images or insufficient seismic impedance contrasts between the units. The interface between unconsolidated sediments and basement at 205 ± 10 m is well defined by: (i) a 1800 to 4500 m s−1 increase in P-wave velocity, (ii) a prominent continuous reflector and (iii) a change from subhorizontal reflectors above to mostly reflector free below.

Figure 12.

Reflection seismic image overlain by tomographic seismic velocity model (colours) computed from high-resolution seismic data acquired along profiles A and B (modified after Podgorski et al.2015). Dashed line marks the interpreted boundary between POM sediments and basement.

Figure 12.

Reflection seismic image overlain by tomographic seismic velocity model (colours) computed from high-resolution seismic data acquired along profiles A and B (modified after Podgorski et al.2015). Dashed line marks the interpreted boundary between POM sediments and basement.

We have assumed that the uncertainties $$\sigma _{z_{r,i}}$$ in eq. (13) are equal to 10 m for the two seismically constrained depths (based primarily on uncertainties in the seismic velocities) and adjusted the Lagrange multiplier γ to explore how closely the seismic constraints can be met. With γ = 10, the two lower boundaries of the resistivity models in Fig. 13 correspond to the seismically defined depths to the POM and basement to within a couple of metres. It is notable that the depths to the two lower boundaries in the models of Fig. 13 are uniform along the line, which contrasts markedly with the erratic models in Fig. 11. Most importantly, the RMS misfits for the seismically constrained 4-layer ML models are identical to the RMS misfits for the 4-layer ML models not constrained by the seismic depths and the same or only slightly larger than the misfits for the R1 models (Table 1), indicating that the electromagnetic and seismic data are fully compatible.

Figure 13.

1-D joint inversion models for AMT, CSAMT and central-loop TEM data along line A (cf. Fig. 3) that are constrained by seismic depths to POM and basement in a modified Marquardt–Levenberg inversion scheme. The layer model parameters, galvanic distortion parameters P and Q for AMT stations and galvanic distortion parameters P and Q for CSAMT stations are depicted in panels (a), (b) and (c), respectively. RMS misfits listed in Table 1 are equal to those of the 4-layer models in Fig. 11 without seismic constraints. Hence, the electromagnetic data can be equally well explained by incorporating the seismic depths to the POM and basement as a priori information.

Figure 13.

1-D joint inversion models for AMT, CSAMT and central-loop TEM data along line A (cf. Fig. 3) that are constrained by seismic depths to POM and basement in a modified Marquardt–Levenberg inversion scheme. The layer model parameters, galvanic distortion parameters P and Q for AMT stations and galvanic distortion parameters P and Q for CSAMT stations are depicted in panels (a), (b) and (c), respectively. RMS misfits listed in Table 1 are equal to those of the 4-layer models in Fig. 11 without seismic constraints. Hence, the electromagnetic data can be equally well explained by incorporating the seismic depths to the POM and basement as a priori information.

Although jointly inverted 3-layer ML models without seismic constraints explain the AMT, CSAMT and TEM data with the same misfits as the jointly inverted 4-layer ML models and the same or slightly greater misfits as the jointly inverted R1 models (Table 1), the depths to the lower resistive layer in these 3-layer models are uniformly much shallower than the basement depth. Once the seismically determined basement depth of 205 ± 10 m is imposed as a constraint on the joint 3-layer ML inversion, the RMS misfits increase by 20–50 per cent. In these 3-layer ML inversions, the largest deterioration in data fit is observed in the low-frequency AMT data and in the CSAMT data between 100–300 Hz (i.e. in the transition zone). We conclude that a minimum of four layers is required to explain the three types of electromagnetic data and at the same time satisfy the seismic constraint on the basement depth.

Model analyses

We have applied the TSVD model error and resolution analysis of Kalscheuer & Pedersen (2007) to the 4-layer ML models for stations 3 and 4 that were inverted without and with seismic constraints on depths to the tops of the POM and basement. In the analyses of both models, the truncation level was set to 21, the number of model parameters (i.e. 7 layer parameters, 10 CSAMT distortion parameters and 4 AMT distortion parameters). To fulfil the prerequisites of the most-squares inversion, we analysed the inversion model that yielded the minimum misfit in both cases (i.e. the model from the 25th iteration in Fig. 10). Approaching the inversion model with minimum misfit, the damping factors of the final Marquardt–Levenberg iterations are much smaller than the singular values of the Jacobian matrix. As a consequence, all seven layer parameters are nearly perfectly resolved. Hence, in assessing how well the layer parameters are constrained, only the model parameter uncertainties need to be further considered.

The most-squares errors $$1/f^-_{{\rm MSQ}}$$ and $$f^+_{{\rm MSQ}}$$ in Table 2 correspond to parameter ranges $$\rho /f^-_{{\rm MSQ}},\ldots ,f^+_{{\rm MSQ}}\cdot \rho$$ and $$t/f^-_{{\rm MSQ}},\ldots ,f^+_{{\rm MSQ}}\cdot {}t$$ of layer resistivities and thicknesses, respectively. Without seismic constraints, resistivity and thickness of the second layer are relatively unstable, such that ρ2 can vary from 0.53 · ρ2 to 1.09 · ρ2 and t2 can vary from 0.56 · t2 to 1.40 · t2. Including seismic constraints on the depths to POM and basement leads to a significant reduction of the uncertainties for layer thicknesses and most layer resistivities (uncertainties of no more than 11 per cent, with most uncertainties less than 4 per cent). With $$1/f^-_{{\rm MSQ}}=0.91$$ and $$f^+_{{\rm MSQ}}=1.11$$ for ρ3, even the resistivity of the third layer (POM) is tightly constrained; the range is 38–47 Ωm for ρ3 = 42 Ωm, such that the resistivity of the POM can be clearly distinguished from the resistivities of the overlying PLM and the underlying basement.

Table 2.

Results of error analyses for resistivities ρ and thicknesses t of the 4-layer models for station 3 in Figs 10 (model ML from the 25th iteration) and 13. Error estimates $$1/f^-_{{\rm MSQ}}$$ and $$f^+_{{\rm MSQ}}$$ were computed using most-squares inversions involving truncated singular value decomposition with a truncation level p = 21, the effective number of model parameters.

Parameter Without seismic constraints With seismic constraints for 
  depths to the POM and basement 
 $$1/f^-_{{\rm MSQ}}$$ $$f^+_{{\rm MSQ}}$$ $$1/f^-_{{\rm MSQ}}$$ $$f^+_{{\rm MSQ}}$$ 
ρ1 0.98 1.02 0.99 1.02 
ρ2 0.53 1.09 0.99 1.01 
ρ3 0.97 1.05 0.91 1.11 
ρ4 0.96 1.04 0.96 1.04 
t1 0.97 1.04 0.99 1.01 
t2 0.56 1.40 0.97 1.03 
t3 0.93 1.06 0.93 1.07 
Parameter Without seismic constraints With seismic constraints for 
  depths to the POM and basement 
 $$1/f^-_{{\rm MSQ}}$$ $$f^+_{{\rm MSQ}}$$ $$1/f^-_{{\rm MSQ}}$$ $$f^+_{{\rm MSQ}}$$ 
ρ1 0.98 1.02 0.99 1.02 
ρ2 0.53 1.09 0.99 1.01 
ρ3 0.97 1.05 0.91 1.11 
ρ4 0.96 1.04 0.96 1.04 
t1 0.97 1.04 0.99 1.01 
t2 0.56 1.40 0.97 1.03 
t3 0.93 1.06 0.93 1.07 

ESTIMATION OF PETROPHYSICAL PROPERTIES

We have estimated plausible ranges of average porosity, clay content and fluid resistivity for the sedimentary units underlying the Wheat Fields site using (i) geological logs of the sediments in boreholes located in southwestern and southeastern regions of the present delta, (ii) layer resistivities of our seismically-constrained 4-layer ML models, (iii) P-wave velocities of the tomographic seismic refraction models (Podgorski et al.2015) and (iv) relevant petrophysical models. Boreholes along the southwestern and southeastern margins of the delta (MMEWR – Ministry of Minerals, Energy and Water Resources 2004) intersect sediments varying from sand to mixtures of sand and clay and units of thick clay. Podgorski et al. (2013b) show the general locations of many boreholes and discuss the implications of the geological information they provide. Fig. 14(a) shows a typical geological borehole log obtained from the southwestern region of the delta (see Fig. 1 for location).

Figure 14.

(a) Geological log of borehole bh in the southwestern part of the Okavango Delta (see Fig. 1 for location) depicting depositional units varying from sand to mixtures of sand and clay and units of thick clay. (b) Electromagnetic joint inversion and seismic models at the Wheat Fields site and their geological interpretation: aquifers of the present Okavango Delta (OD), lacustrine sediments of Paleo Lake Makgadikgadi (PLM), a freshwater bearing aquifer of the Paleo Okavango Megafan (POM) and basement. Note, borehole bh is about 110 km from the Wheat Fields site.

Figure 14.

(a) Geological log of borehole bh in the southwestern part of the Okavango Delta (see Fig. 1 for location) depicting depositional units varying from sand to mixtures of sand and clay and units of thick clay. (b) Electromagnetic joint inversion and seismic models at the Wheat Fields site and their geological interpretation: aquifers of the present Okavango Delta (OD), lacustrine sediments of Paleo Lake Makgadikgadi (PLM), a freshwater bearing aquifer of the Paleo Okavango Megafan (POM) and basement. Note, borehole bh is about 110 km from the Wheat Fields site.

Petrophysical models

To explain the resistivities of the sedimentary layers, we have attempted to relate total resistivity to interconnected porosity ϕi, clay content c, surface conduction of clay minerals and pore fluid resistivity ρf using the differential effective medium model introduced by Revil et al. (1998). This model includes contributions from both anionic and cationic conduction. Whereas anions and cations migrate through the electrolyte of interconnected pore space, only cations also migrate through the electrical double layer covering clay mineral surfaces. Electrolytic conduction is dependent on the porosity ϕi and cementation constant m of the interconnected pore space and pore fluid (electrolyte) conductivity. Surface conductivity through the electrical double layers is determined by the surface mobility of the cations and by the cation exchange capacity (CEC) of the clay minerals. These are well-determined material properties (e.g. Revil et al.1998). Following Patchett (1975) and Revil et al. (1998), the effective CEC is the product of the mass fraction of clay minerals and the arithmetic mean of the different CECs weighted by the relative fractions of the different clay types. We consider only one clay type at a time for simplicity. For estimating pore fluid and surface conductivities, we assume a surface temperature of 25 $$^{\circ }{\rm C}$$ and a temperature increase of 2.5 $$^{\circ }{\rm C}$$/100 m. The effects of temperature change are negligible for the limited depth range that we need to consider. Interconnected porosity ϕi, cementation constant m, clay type, clay content c and fluid resistivity ρf, which is a measure of salinity, need to be provided for each layer.

For the seismic velocities, we adopt Marion et al.'s (1992) model for unconsolidated clay and sand mixtures. Several assumptions need to be made in applying their model. Bulk and shear moduli of clayey sands (c < ϕs, where ϕs is pure sand porosity) and sandy clays (c > ϕs) are estimated via Gassmann's relations and Reuss’ averages, respectively. Well-constrained estimates of elastic moduli and densities of the individual constituents (i.e. sand grains, clay minerals and pore fluid) are provided by Marion et al. (1992) and Mavko et al. (2009). Total porosity ϕt is computed as a function of pure sand porosity ϕs , pure clay porosity ϕc and c. For ϕs and ϕc, we assume exponential decays with depth and use the laboratory values at different pressures determined by Marion et al. (1992). In Marion et al.'s (1992) clayey sand model, we use dry-frame moduli for unconsolidated sands based on the modified Hashin-Shtrikman lower bound (Dvorkin & Nur 1996). According to Dvorkin & Nur (1996), sand mineralogy (i.e. relative content of quartz, mica and feldspar) has a negligible effect on P-wave velocities, and we found that clay type did not change the modelled P-wave velocities of clayey sands by more than a few per cent. Hence, for clayey sand, fractional clay content c is the main parameter to adjust. For sandy clay, P-wave velocity can vary by a factor of two by changing the clay type, such that both clay type and clay content are important for this type of sediment. Seismic wave velocities and total porosities ϕt are output values of Marion et al.'s (1992) model. Besides the adopted model, clay type and clay content, the porosities of the sand and clay play key roles in determining the P-wave velocities. Indeed, the porosities are the largest potential sources of error in our petrophysical interpretation. To obtain a reasonably good fit to the tomographic P-wave velocities, we only need to vary sand and clay porosities by a few per cent from the values suggested by Marion et al. (1992).

The following relationships need to be considered when comparing and evaluating the assumed and derived petrophysical properties listed in Table 3:

  1. Whereas electrical resistivity is only dependent on ϕi, velocity is dependent on ϕt. Estimated porosities are consistent only if ϕi < ϕt.

  2. For clayey sands with roughly c < 0.15, small increases in clay content or changes in clay mineralogy may result in marked reductions in electrical resistivity, but only small increases in velocity. Hence, for such sediments the upper limit of clay content is primarily determined by electrical resistivity.

Table 3.

Ranges of petrophysical parameters that explain (a) the electrical resistivities ρ and (b) the seismic P-wave velocities vp of the sedimentary units using Revil et al.'s (1998) and Marion et al.'s (1992) petrophysical models for sand and clay mixtures. c, clay content; CEC, cation exchange capacity; ϕi, interconnected porosity; m, cementation constant; ρf, pore fluid resistivity; ϕt, total porosity.

(a) Electrical resistivity 
Layer Geophysical model Petrophysical model 
 ρ [Ωm] ρ [Ωm] c CEC [meq g−1ϕi m ρf [Ωm] 
Shallow aquifer 85–100 85–100 <0.13 0.023–0.05 0.10–0.32 1.3–1.9 20–200 
PLM 15–20 20–22 0.15–0.27 0.023–0.05 0.15–0.27 1.3–1.9 10–200 
  14–20 0.50–0.60 0.023 0.01–0.05 1.9 >0.3 
POM 40–50 40–50 <0.15 0.023–0.05 0.10–0.33 1.3–1.9 10–200 
(b) Seismic P-wave velocity 
Layer Geophysical model Petrophysical model 
 vp [ms− 1vp [ms− 1ϕt c    
Shallow aquifer 800–1800 1630–1770 0.29–0.35 <0.13    
PLM ∼1800 1780–1840 0.22–0.27 0.15–0.27    
  1770–1820 0.27–0.32 0.50–0.60    
POM 1700–2100 1650–1850 0.26–0.35 <0.15    
(a) Electrical resistivity 
Layer Geophysical model Petrophysical model 
 ρ [Ωm] ρ [Ωm] c CEC [meq g−1ϕi m ρf [Ωm] 
Shallow aquifer 85–100 85–100 <0.13 0.023–0.05 0.10–0.32 1.3–1.9 20–200 
PLM 15–20 20–22 0.15–0.27 0.023–0.05 0.15–0.27 1.3–1.9 10–200 
  14–20 0.50–0.60 0.023 0.01–0.05 1.9 >0.3 
POM 40–50 40–50 <0.15 0.023–0.05 0.10–0.33 1.3–1.9 10–200 
(b) Seismic P-wave velocity 
Layer Geophysical model Petrophysical model 
 vp [ms− 1vp [ms− 1ϕt c    
Shallow aquifer 800–1800 1630–1770 0.29–0.35 <0.13    
PLM ∼1800 1780–1840 0.22–0.27 0.15–0.27    
  1770–1820 0.27–0.32 0.50–0.60    
POM 1700–2100 1650–1850 0.26–0.35 <0.15    

Results of petrophysical modelling

Petrophysical property ranges that are compatible with our preferred P-wave velocity and resistivity models (Figs 12 & 13 and Table 3) are explained for each layer in the following.

The 85–100 Ωm resistivity of the shallow aquifer can be explained by a clay content c < 0.13, CEC values of 0.023–0.05 meq g−1 (representative of kaolinite, the predominant clay mineral currently being transported into the delta; McCarthy 2013), interconnected porosity ϕi = 0.10–0.32, cementation factor m = 1.3–1.9 (typical of well-rounded sand grains to more disk-shaped particles; Mavko et al.2009; Schön 2011) and pore fluid resistivity ρf = 20–200 Ωm (Table 3a). The higher values of ρf are typical of lacustrine and riverine waters. For pure unconsolidated sand (c = 0 and m = 1.3), ϕi = 0.32 combined with ρf = 20 Ωm yields a resistivity value similar to that of our preferred inversion models. Clay contents larger than 0.13 reduce the total resistivity below the range observed in these models. We set the maximum clay content to 0.13 for our petrophysical analysis of P-wave velocity. This choice results in total porosities ϕt = 0.29–0.35 and P-wave velocities vp = 1630–1770 m s−1 that are in good agreement with the upper end of the range observed in the tomographic seismic refraction models (Table 3b).

The borehole geological logs suggest that PLM comprises intercalated clayey sands and thick units of clay. Hence, we consider both clayey sands and sandy clays in the petrophysical modelling. For clayey sand, c = 0.15–0.27, CEC values of 0.023–0.05 meq g−1, ϕi = 0.15–0.27, m = 1.3–1.9 and ρf = 10–200 Ωm give reasonable fits to the inversion model resistivities of 15–20 Ωm. Because of the relatively high clay content, the choices of ϕi and m have less effect on the petrophysically modelled resistivity than for the shallow aquifer. For sandy clay, c = 0.5–0.6, ϕi = 0.01–0.05, m = 1.9 and ρf > 0.3 Ωm can explain the inverted resistivity range when a low CEC of 0.023 meq g−1 (corresponding to kaolinite) is assumed. For such high clay contents, changes of ρf have little effect on the total resistivity. The assumption that both the clayey sands and clay PLM units have resistivities of roughly 15–20 Ωm would explain why there is no evidence for intra-PLM layering or anisotropy in our electromagnetic data. The tomographic P-wave velocity of PLM (vp ∼ 1800 m s−1) can be explained by either Marion et al.'s clayey sand model with c = 0.15–0.27 and ϕt = 0.22–0.27 or their sandy clay model with c = 0.5–0.6 and ϕt = 0.27–0.32. Note that also brackish water saturated sands with c < 0.05, CEC ∼ 0.05 meq g−1, ϕi ∼ 0.3, m = 1.3 and ρf ∼ 4 Ωm would have a resistivity of 15–20 Ωm. The resulting P-wave velocities ∼1700 m s−1 would be only slightly smaller than the ∼1800 m s−1 of the tomographic seismic models.

For the proposed POM, c < 0.15, CEC values of 0.023–0.05 meq g−1, ϕi = 0.10–0.33, m = 1.3–1.9 (representative of loose sediments) and ρf = 10–200 Ωm yield resistivities compatible with the 40–50 Ωm resistivities in our preferred inversion models. Assuming pure unconsolidated sand (c = 0 and m = 1.3), a fluid resistivity of 10 Ωm and a porosity of 0.33 would give a total resistivity comparable to that of these models. An average clay content much larger than 0.15 is not feasible, because it would reduce the electrical resistivity below the range observed in the inversion models. Setting the upper limit of clay content to 0.15 yields ϕt = 0.26–0.35 and vp = 1650–1850 m s−1 using Marion et al.'s clayey sand model.

Due to compensating effects, there exists considerable ambiguity in the choice of material parameters. Potentially meaningful parameter ranges are larger than the ones given here.

DISCUSSION

Resolving equivalence in multiple electromagnetic data sets acquired within the Okavango Delta

It has long been known that ambiguities in electromagnetic, geoelectric and seismic models can be markedly reduced by combining information from multiple types of data. Traditional approaches for mitigating equivalence and suppression in resistivity models have involved joint inversions of electromagnetic data generated by inductively coupled sources and geoelectric data generated by galvanically coupled sources (Vozoff & Jupp 1975; Raiche et al.1985; Albouy et al.2001; Kalscheuer et al.2010). Unfortunately, the geoelectric data collected at the Wheat Fields site was compromised by coupling problems and by current channelling in the conductive second layer that severely limited depth penetration. Although inductive and galvanic sources were used to generate our CSAMT data, no significant reduction in equivalence for the interpreted POM was observed in inversions of the CSAMT data or in joint inversions of the AMT, CSAMT and TEM data.

Rather than jointly inverting our multiple electromagnetic and coincident seismic data sets, we developed and applied a joint inversion scheme to the electromagnetic data that produced resistivity models with depths to layer boundaries required to mimic the depths to well-determined seismic interfaces. The resulting 4-layer resistivity models with layer boundaries compatible with the seismic reflection images and tomographic seismic velocity models explain the observed electromagnetic data as well as the best resistivity models not constrained by the seismic data. Furthermore, the resistivities of the two deepest layers in the new 4-layer models appeared to be better constrained than in any other models.

To test the sensitivity of the joint inversion scheme to the interpreted depth to the POM, we jointly inverted the electromagnetic data using other depth values. For example, a 135 m deep seismic reflector that is less continuous and less prominent than that at 145 m depth could be the PLM–POM boundary. By assuming the POM to be at 135 m depth in the joint inversions for stations 3 and 4, the resistivity of the third layer decreased from 40–50 Ωm to 30–40 Ωm whereas other layer resistivities were largely unaffected remaining within statistical error of the model with the depth to the POM at 145 m. At the other extreme, by setting the POM depth to 170 m (about two seismic wavelengths deeper than the preferred seismic interpretation), the resistivity of the third layer (i.e. the proposed POM) increased to 150 Ωm with a relatively large nonlinear uncertainty range of 90–405 Ωm computed using most-squares inversion. The RMS misfits between the observed and model-predicted data and all other resistivities and depths and their confidence limits remained essentially the same as for the preferred model. Even for the relatively extreme 170 m depth estimate, the resistivity of the third layer was clearly distinguishable from the resistivities of the overlying and underlying layers. Once the basement depth was firmly established, the depth to the top of the POM had a major influence on the estimated resistivity of the POM and its uncertainty range.

Geological and petrophysical interpretation

For the north-central Okavango Delta, the jointly inverted 4-layer ML models constrained by the seismically determined depths (Fig. 13) are interpreted from top to bottom as (Fig. 14b):

  1. Layer 1: an upper freshwater aquifer with 85–100 Ωm average resistivity, 20–200 Ωm pore fluid resistivity, ∼40 m average thickness, 10–32 per cent average interconnected porosity and ≤13  per cent average clay content that is affected by annual rainfall, inflow of the Okavango River and evapotranspiration.

  2. Layer 2: intercalated lacustrine clays and clayey sands with a low 15–20 Ωm resistivity, 105–115 m thickness that is the likely remnants of PLM (Podgorski et al.2013b). Interconnected porosity and clay content likely vary over ranges of 1–27 per cent and 15–60 per cent, respectively. The upper part of this layer may include deltaic sediments of the present delta that are saturated with brines from downward diffusion and density-driven flow (McCarthy & Metcalfe 1990; Bauer et al.2006; Zimmermann et al.2006; Milzow et al.2009).

  3. Layer 3: a freshwater aquifer with a moderate 40–50 Ωm resistivity, 10–200 Ωm pore fluid resistivity, ∼60 m thickness, 10–33 per cent average interconnected porosity and ≤15  per cent average clay content that is interpreted to be the remains of the POM (Podgorski et al.2013b). Diffusion or flow of salt water into this aquifer appears to be blocked by clay beds within the PLM sediments.

  4. Layer 4 (basal half-space): bedrock with resistivities of 3200–4000 Ωm.

The geological interpretation of the upper layering of our Wheat Fields geophysical models is comparable to the sedimentary sequence (Fig. 14a) recorded down to 92 m depth in borehole bh (location shown in Fig. 1) in the southwestern region of the delta. In this borehole, clay units are observed at depths >84 m, which is within the depth range of PLM at the Wheat Fields site. The very fine to fine grained sands immediately overlying this clay unit may correspond to the salt-water saturated sands of the present delta that we interpret as the top part of the conductive second layer at our field site. It should be borne in mind though that the borehole is 110 km from our field site (Fig. 1). Several other boreholes drilled in the same area as bh indicate that depositional conditions are spatially variable.

Our findings corroborate the evolutionary model for PLM and POM presented by Podgorski et al. (2013b). This model assumes that crustal flexuring and subsidence resulted in the redirection of major rivers and formation of PLM within the Kalahari Basin. Originally, the POM and an early southeastern part of PLM may have co-existed. Continued subsidence of the Kalahari Basin and displacements along the Kunyere and Thamalakane faults likely caused PLM to expand further north. As a consequence, the POM was buried underneath lacustrine sediments of PLM. Later tectonism may have re-directed major river systems resulting in the reduction of water flowing into PLM. Eventually, this lead to the desiccation of PLM and formation of the present Okavango Delta. As emphasized by Podgorski et al. (2013b), the timing of these events is uncertain.

Though the current environmental policies of the Botswanan government limit exploration, groundwater stored in POM may be a freshwater resource for the Botswanan population in the distant future.

CONCLUSIONS

Our work on the resistivity structure of the sediments and basement beneath the Okavango Delta began with the processing and inversion of HTEM data (Podgorski et al.2013a,b). Models derived from these data contained three units with a conductive layer sandwiched between overlying and underlying resistive layers. In seeking to constrain the interpretation of these models, we first acquired ERT, ground-based TEM, seismic reflection and tomographic seismic refraction data at various locations within the delta. Independent and joint inversions of the ERT and ground-based TEM data yielded 3-layer resistivity models that were essentially the same as the models based on the HTEM data. The depth to basement (i.e. the base of the sedimentary layers) was well defined in the seismic reflection images and tomographic seismic velocity models. Along the western margin of the delta, the top of the highly resistive layer (i.e. the half-space) in the airborne and ground-based electromagnetic models coincided with the seismically determined basement. In contrast, the top of the resistive layer was much shallower than the seismically determined basement depth in the north-central region of the delta. Based on these results, we inferred the existence of a resistive sedimentary layer (i.e. the POM) between the conductive layer and resistive basement in this region.

To examine the possibility of a deep resistive sedimentary layer overlying basement, we acquired deeper penetrating AMT and CSAMT data, the latter of which involved signals generated by both inductive and galvanic sources. Individual and joint inversions with the ground-based TEM data demonstrated that the AMT and CSAMT data could be explained by 3-layer ML models that were essentially the same as the 3-layer ML models derived from the TEM data alone, such that the depth to the basal resistive layer in the 3-layer ML models was again much shallower than the well-defined seismically-determined basement depth. The addition of the AMT and CSAMT did not resolve the equivalence problem for the POM.

Our next steps involved (i) imposing the seismically-determined basement depth on 3-layer joint inversions of the AMT, CSAMT and TEM data, (ii) jointly inverting for four layers without constraining the depths to any boundaries and (iii) jointly inverting for four layers with the seismically determined depths to the POM and basement as explicit constraints. Requiring a match between the depth to the basal resistive layer in the 3-layer joint inversions and the seismically determined basement depth resulted in 3-layer ML models with unacceptably high RMS misfits with the observed data. Although the upper part of the basal resistive layer in the 3-layer ML models was not basement, estimates of the resistivity and depth of the basal resistive layer were undoubtedly affected by the electromagnetic response of the basement. Allowing an additional layer in the joint inversions without any depth constraints produced 4-layer ML models with depths to the basal layer that were in the same range as the seismically determined basement depth, but with significant trade-offs between the resistivities and depths of the second and third layers. By requiring the depths to the lower two boundaries of the resistivity models to closely match the seismically determined depths to the POM and basement, the joint inversions yielded 4-layer ML models with well-determined resistivities and depths for all layers and RMS misfits that were as good as the unconstrained 3- and 4-layer ML models. The 4-layer ML models, which explain all ground-based electromagnetic and seismic data and are compatible with the airborne electromagnetic data, provide the first robust estimates for the intermediate resistivity of the proposed POM and the high resistivity of the basement beneath the Okavango Delta.

Our results support a simple layered model of the north-central Okavango Delta in which the top resistive layer of open water and freshwater-saturated deltaic sediments is underlain sequentially by a conductive layer of clay and saline-water-saturated sands (PLM), a freshwater-saturated gravel/sand layer with intermediate resistivities (POM) and a highly resistive basement. Furthermore, our resistivity models support the interpretation of the HTEM resistivity models that PLM extended across the entire area of the present Okavango Delta, such that the paleo lake once covered a total area of at least 90 000 km2.

For sedimentary environments like the Okavango Delta, petrophysical modelling of seismic P-wave velocities and electrical resistivities suffers from significant trade-offs between porosity and clay content. To obtain narrower ranges of porosity and clay content for the upper sedimentary units, it would be desirable to collect nuclear magnetic resonance and induced polarization data. Improved models of hydraulic permeability are needed to facilitate meaningful hydrological modelling of groundwater flow. On a local scale, such improvements could be obtained by drilling a borehole down to basement and logging distributions of porosity and grain size diameter. Owing to the intrinsic loss of resolution with depth, surface geophysical measurements alone may not be able to provide sufficiently accurate estimates of physical properties in the deeper parts of the sedimentary section.

This work was funded by the Swiss National Science Foundation under grant number 126981. Wilderness Safaris, Maun, provided invaluable logistical support. We would like to thank our colleagues Bashali Charles Jaba and Baphaleng Bantshang from the Department of Geological Survey of Botswana and F. Reiser and C. Bärlocher of ETH Zürich for helping in the field campaign. Niels Bøie Christensen, Aarhus University, kindly provided the fast Hankel transform routines employed in the Fourier sine transform to compute TEM responses. We thank journal editor Mark Everett and three anonymous reviewers for their constructive comments.

REFERENCES

Albouy
Y.
Andrieux
P.
Rakotondrasoa
G.
Ritz
M.
Descloitres
M.
Join
J.-L.
Rasolomanana
E.
Mapping coastal aquifers by joint inversion of DC and TEM soundings—three case histories
Ground Water
 
2001
39
1
87
97
Asten
M.W.
Full transmitter waveform transient electromagnetic modeling and inversion for soundings over coal measures
Geophysics
 
1987
52
3
279
288
Auken
E.
Christiansen
A.V.
Layered and laterally constrained 2D inversion of resistivity data
Geophysics
 
2004
69
3
752
761
Auken
E.
Christiansen
A.V.
Jacobsen
L.H.
Sørensen
K.I.
A resolution study of buried valleys using laterally constrained inversion of TEM data
J. Appl. Geophys.
 
2008
65
1
10
20
Bailey
R.C.
Inversion of the geomagnetic induction problem
Proc. R. Soc. A.
 
1970
315
185
194
Bastani
M.
EnviroMT—a new controlled source/radio magnetotelluric system
PhD thesis
 
2001
Uppsala University
Bauer
P.
Supper
R.
Zimmermann
S.
Kinzelbach
W.
Geoelectrical imaging of groundwater salinization in the Okavango Delta, Botswana
J. Appl. Geophys.
 
2006
60
2
126
141
Bauer-Gottwein
P.
Langer
T.
Prommer
H.
Wolski
P.
Kinzelbach
W.
Okavango Delta Islands: interaction between density-driven flow and geochemical reactions under evapo-concentration
J. Hydrol.
 
2007
335
3–4
389
405
Bauer-Gottwein
P.
Gondwe
B.N.
Christiansen
L.
Herckenrath
D.
Kgotlhang
L.
Zimmermann
S.
Hydrogeophysical exploration of three-dimensional salinity anomalies with the time-domain electromagnetic method (TDEM)
J. Hydrol.
 
2010
380
318
329
Berdichevsky
M.N.
Dmitriev
V.I.
Models and Methods of Magnetotellurics
 
2008
Springer
Boerner
D.E.
Wright
J.A.
Thurlow
J.G.
Reed
L.E.
Tensor CSAMT studies at the Buchans Mine in central Newfoundland
Geophysics
 
1993
58
1
12
19
Büker
F.
Green
A.G.
Horstmeyer
H.
Shallow seismic reflection study of a glaciated valley
Geophysics
 
1998
63
4
1395
1407
Büker
F.
Green
A.G.
Horstmeyer
H.
3-D high-resolution reflection seismic imaging of unconsolidated glacial and glaciolacustrine sediments: processing and interpretation
Geophysics
 
2000
65
1
18
34
Burrough
S.L.
Thomas
D.S.G.
Bailey
R.M.
Mega-Lake in the Kalahari: a late Pleistocene record of the Palaeolake Makgadikgadi system
Quat. Sci. Rev.
 
2009
28
15–16
1392
1411
Chave
A.D.
Chave
A.D.
Jones
A.G.
Estimation of the magnetotelluric response function, in
The Magnetotelluric Method: Theory and Practice
 
2012
Cambridge Univ. Press
165
218
chap. 5
Christensen
N.B.
Optimized fast Hankel transform filters
Geophys. Prospect.
 
1990
38
5
545
568
Christiansen
A.V.
Auken
E.
Viezzoli
A.
Quantification of modeling errors in airborne TEM caused by inaccurate system description
Geophysics
 
2011
76
1
F43
F52
Constable
S.C.
Parker
R.L.
Constable
C.G.
Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data
Geophysics
 
1987
52
3
289
300
Dvorkin
J.
Nur
A.
Elasticity of high-porosity sandstones: theory for two North Sea data sets
Geophysics
 
1996
61
5
1363
1370
Effersø
F.
Auken
E.
Sørensen
K.I.
Inversion of band-limited TEM responses
Geophys. Prospect.
 
1999
47
4
551
564
Fitterman
D.V.
Anderson
W.L.
Effect of transmitter turn-off time on transient soundings
Geoexploration
 
1987
24
2
131
146
Garcia
X.
Jones
A.G.
Atmospheric sources for audio-magnetotelluric (AMT) sounding
Geophysics
 
2002
67
2
448
458
Garcia
X.
Boerner
D.
Pedersen
L.B.
Electric and magnetic galvanic distortion decomposition of tensor CSAMT data. Application to data from the Buchans Mine (Newfoundland, Canada)
Geophys. J. Int.
 
2003
154
3
957
969
Geonics Ltd.
Tx Turn-on and Primary Field Measurement in GEONICS TEM Systems, unpublished technical manual
 
1998
Geonics Ltd.
PROTEM System Response
Unpublished Technical manual
 
2011
Gieske
A.
Vegetation driven groundwater recharge below the Okavango Delta (Botswana) as a solute sink mechanism: an indicative model
Botsw. J. Earth Sci.
 
1996
3
33
37
Gumbricht
T.
McCarthy
T.S.
Merry
C.L.
The topography of the Okavango Delta, Botswana, and its tectonic and sedimentological implications
S. Afr. J. Geol.
 
2001
104
3
243
264
Jackson
D.D.
Most Squares Inversion
J. Geophys. Res.
 
1976
81
5
1027
1030
Jones
A.G.
Imaging the continental upper mantle using electromagnetic methods
Lithos
 
1999
48
1–4
57
80
Jones
A.G.
Electromagnetic interrogation of the anisotropic Earth: looking into the Earth with polarized spectacles
Phys. Earth Planet. Inter.
 
2006
158
2–4
281
291
Jones
A.G.
Chave
A.D.
Jones
A.G.
Distortion of magnetotelluric data: its identification and removal
The Magnetotelluric Method: Theory and Practice
 
2012
Cambridge Univ. Press
219
302
chap. 6
Kalscheuer
T.
Pedersen
L.B.
A non-linear truncated SVD variance and resolution analysis of two-dimensional magnetotelluric models
Geophys. J. Int.
 
2007
169
2
435
447
Kalscheuer
T.
Commer
M.
Hördt
A.
Helwig
S.L.
Tezkan
B.
Electromagnetic Evidence for an Ancient Avalanche Caldera Rim on the South Flank of Mount Merapi, Indonesia
J. Volc. Geotherm. Res.
 
2007
162
1–2
81
97
Kalscheuer
T.
García
M.
Meqbel
N.
Pedersen
L.B.
Non-linear model error and resolution properties from two-dimensional single and joint inversions of direct current resistivity and radiomagnetotelluric data
Geophys. J. Int.
 
2010
182
3
1174
1188
Kalscheuer
T.
Hübert
J.
Kuvshinov
A.
Lochbühler
T.
Pedersen
L.B.
A hybrid regularization scheme for the inversion of magnetotelluric data from natural and controlled sources to layer and distortion parameters
Geophysics
 
2012
77
4
E301
E315
Kinabo
B.D.
Hogan
J.P.
Atekwana
E.A.
Abdelsalam
M.G.
Modisi
M.P.
Fault growth and propagation during incipient continental rifting: insights from a combined aeromagnetic and Shuttle Radar Topography Mission digital elevation model investigation of the Okavango Rift Zone, northwest Botswana
Tectonics
 
2008
27
3
TC3013
doi:10.1029/2007TC002154
Li
X.
Pedersen
L.B.
Controlled source tensor magnetotellurics
Geophysics
 
1991
56
9
1456
1461
Lines
L.R.
Treitel
S.
Tutorial: A review of least-squares inversion and its application to geophysical problems
Geophys. Prospect.
 
1984
32
2
159
186
Marion
D.
Nur
A.
Yin
H.
Han
D.
Compressional velocity and porosity in sand-clay mixtures
Geophysics
 
1992
57
4
554
563
Mavko
G.
Mukerji
T.
Dvorkin
J.
The Rock Physics Handbook: Tools for Seismic Analysis in Porous Media
 
2009
2nd edn
Cambridge Univ. Press
McCarthy
T.S.
Groundwater in the wetlands of the Okavango Delta, Botswana, and its contribution to the structure and function of the ecosystem
J. Hydrol.
 
2006
320
3–4
264
282
McCarthy
T.S.
The Okavango Delta and its place in the geomorphological evolution of southern Africa
S. Afr. J. Geol.
 
2013
116
1
1
54
McCarthy
T.S.
Metcalfe
J.
Chemical sedimentation in the semi-arid environment of the Okavango Delta, Botswana
Chem. Geol.
 
1990
89
1–2
157
178
McCarthy
T.S.
McIver
J.R.
Cairncross
B.
Ellery
W.N.
Ellery
K.
The inorganic chemistry of peat from the Maunachira channel-swamp system, Okavango Delta, Botswana
Geochim. Cosmochim. Acta
 
1989
53
5
1077
1089
McCarthy
T.S.
McIver
J.R.
Verhagen
B.T.
Groundwater evolution, chemical sedimentation and carbonate brine formation on an island in the Okavango Delta swamp, Botswana
Appl. Geochem.
 
1991
6
6
577
595
McCarthy
T.S.
Green
R.W.
Franey
N.J.
The influence of neo-tectonics on water dispersal in the northeastern regions of the Okavango swamps, Botswana
J. Afr. Earth Sci.
 
1993
17
1
23
32
Meier
P.
et al
Hydrogeophysical investigations in the western and north-central Okavango Delta (Botswana) based on helicopter and ground-based transient electromagnetic data and electrical resistance tomography
Geophysics
 
2014
79
5
B201
B211
Meju
M.A.
Biased estimation: a simple framework for inversion and uncertainty analysis with prior information
Geophys. J. Int.
 
1994
119
2
521
528
Meju
M.A.
Hutton
V.R.S.
Iterative most-squares inversion: application to magnetotelluric data
Geophys. J. Int.
 
1992
108
3
758
766
Milzow
C.
Kgotlhang
L.
Bauer-Gottwein
P.
Meier
P.
Kinzelbach
W.
Regional review: the hydrology of the Okavango Delta, Botswana – processes, data and modelling
Hydrogeol. J.
 
2009
17
6
1297
1328
MMEWR (Ministry of Minerals Energy and Water Resources)
Maun groundwater development project, Phase - 2, Resources Assessment and Wellfield Development (TB 10/3/5/2000–2001), Final Report, Volume 6 - Hydrogeological report, Part 4 : Gomoti Exploration Area, Tech. rep., Department of Water Affairs
2004
Gaborone, Botswana
Republic of Botswana
Modisi
M.
Atekwana
E.
Kampunzu
A.
Ngwisanyi
T.
Rift kinematics during the incipient stages of continental extension: evidence from the nascent Okavango rift basin, northwest Botswana
Geology
 
2000
28
10
939
942
Moore
A.E.
Larkin
P.A.
Drainage evolution in south-central Africa since the breakup of Gondwana
S. Afr. J. Geol.
 
2001
104
1
47
68
Moore
A.E.
Cotterill
F.P.D.
Eckardt
F.
The evolution and ages of Makgadikgadi palaeo-lakes: Consilient evidence from Kalahari drainage evolution
S. Afr. J. Geol.
 
2012
115
385
413
Nabighian
M.N.
Macnae
J.C.
Nabighian
M.N.
Time domain electromagnetic prospecting methods, in
Electromagnetic Methods in Applied Geophysics
 
1991
427
520
chap. 6
SEG
Nitsche
F.O.
Green
A.G.
Horstmeyer
H.
Büker
F.
Late Quaternary depositional history of the Reuss delta, Switzerland: constraints from high-resolution seismic reflection and georadar surveys
J. Quat. Sci.
 
2002
17
2
131
143
Patchett
J.G.
An investigation of shale conductivity
The Log Analyst
 
1975
16
06
3
20
Piessens
R.
De Doncker-Kapenga
E.
Überhuber
C.W.
QUADPACK: A Subroutine Package for Automatic Integration, Vol 1 of Series in Computational Mathematics
 
1983
Springer Verlag
Podgorski
J.E.
Auken
E.
Schamper
C.
Christiansen
A.V.
Kalscheuer
T.
Green
A.G.
Processing and inversion of commercial helicopter time-domain electromagnetic data for environmental assessments and geological and hydrological mapping
Geophysics
 
2013a
78
4
E149
E159
Podgorski
J.E.
Green
A.G.
Kgotlhang
L.
Kinzelbach
W.K.H.
Kalscheuer
T.
Auken
E.
Ngwisanyi
T.
Paleo megalake and paleo megafan in southern Africa
Geology
 
2013b
41
11
1155
1158
Podgorski
J.E.
et al
Integrated interpretation of helicopter and ground-based geophysical data recorded within the Okavango Delta, Botswana
J. Appl. Geophys.
 
2015
114
52
67
Qian
W.
On small-scale near-surface distortion in controlled source tensor electromagnetics
Geophys. Prospect.
 
1994
42
5
501
520
Raiche
A.P.
Jupp
D.L.B.
Rutter
H.
Vozoff
K.
The joint use of coincident loop transient electromagnetic and Schlumberger sounding to resolve layered structures
Geophysics
 
1985
50
10
1618
1627
Reeves
C.V.
Rifting in the Kalahari
Nature
 
1972
237
95
96
Reiser
F.
et al
Constraining helicopter electromagnetic models of the Okavango Delta with seismic-refraction and seismic-reflection data
Geophysics
 
2014
79
3
B123
B134
Revil
A.
Cathles
L.M.
Losh
S.
III
Nunn
J.A.
Electrical conductivity in shaly sands with geophysical applications
J. Geophys. Res.
 
1998
103
B10
23 925
23 936
Ringrose
S.
et al
Diagenesis in Okavango fan and adjacent dune deposits with implications for the record of palaeo-environmental change in Makgadikgadi-Okavango-Zambezi basin, northern Botswana
Geomorphology
 
2008
101
4
544
557
Routh
P.S.
Oldenburg
D.W.
Inversion of controlled source audio-frequency magnetotellurics data for a horizontally layered earth
Geophysics
 
1999
64
6
1689
1697
Scholz
C.H.
Koczynski
T.A.
Hutchins
D.G.
Evidence for Incipient Rifting in Southern Africa
Geophys. J. R. astr. Soc.
 
1976
44
135
144
Schön
J. H.
Physical Properties of Rocks: A Workbook, Vol. 8 of Handbook of Petroleum Exploration and Production
 
2011
Elsevier
Shaw
P.A.
Thomas
D.S.G.
Geomorphology, sedimentation, and tectonics in the Kalahari Rift
Isr. J. Earth Sci.
 
1992
41
87
94
Simpson
F.
Bahr
K.
Practical Magnetotellurics
 
2005
Cambridge Univ. Press
Spies
B.R.
Depth of investigation in electromagnetic sounding methods
Geophysics
 
1989
54
7
872
888
Sternberg
B.K.
The variability of naturally occurring magnetic field levels: 10 Hz to 8 kHz
Geophysics
 
2010
75
6
F187
F197
Thomas
D.S.G.
Shaw
P.A.
The Kalahari Environment
 
1991
Cambridge Univ. Press
Vozoff
K.
Nabighian
M.N.
The magnetotelluric method
Electromagnetic Methods in Applied Geophysics
 
1991
SEG
641
711
chap. 8
Vozoff
K.
Jupp
D.L.B.
Joint Inversion of Geophysical Data
Geophys. J. R. astr. Soc.
 
1975
42
3
977
991
Wannamaker
P.E.
Tensor CSAMT survey over the Sulphur Springs thermal area, Valles Caldera, New Mexico, U.S.A., Part II: Implications for CSAMT methodology
Geophysics
 
1997
62
2
466
476
Wannamaker
P.E.
Hohmann
G.W.
Ward
S.H.
Magnetotelluric responses of three-dimensional bodies in layered earths
Geophysics
 
1984
49
9
1517
1533
Ward
S.H.
Hohmann
G.W.
Nabighian
M.N.
Electromagnetic Theory for Geophysical Applications, in
Electromagnetic Methods in Applied Geophysics
 
1987
1
SEG
131
311
Theory, chap. 4
Zimmermann
S.
Bauer
P.
Held
R.
Kinzelbach
W.
Walther
J.H.
Salt transport on islands in the Okavango Delta: numerical investigations
Adv. Water Resour.
 
2006
29
1
11
29
Zonge
K.L.
Hughes
L.J.
Nabighian
M.N.
Controlled source audio-frequency magnetotellurics, in
Electromagnetic Methods in Applied Geophysics
 
1991
SEG
713
809
chap. 9

APPENDIX A: COMPUTATION OF TEM FORWARD RESPONSES AND SENSITIVITIES

In our approach to modelling TEM data, the electromagnetic fields and their sensitivities are initially computed in the frequency domain and then transformed to the time domain using a fast Hankel transform filter (Christensen 1990). In the frequency domain, we have based the computation of electromagnetic fields and sensitivities due to current excitation in a large horizontal loop on the representation of an extended cable source by a set of horizontal electric dipoles distributed along the cable (e.g. Ward & Hohmann 1987). To some extent similar to Kalscheuer et al. (2012), we have removed the terms arising from galvanic coupling of the individual dipoles such that only the inductively coupled terms of the dipole fields are retained. For a large horizontal transmitter loop and arbitrary transmitter-receiver geometry, fields and sensitivities are computed using adaptive Gauss–Kronrod quadrature (Piessens et al.1983) of the inductively coupled parts along all four sides of the loop. For a central-loop configuration, symmetry permits us to integrate the Hz field and its sensitivities only along one half of one loop side and then to multiply the result by eight. As for the CSAMT case, fields and sensitivities of individual dipoles are computed using fast Hankel transforms involving Wait's recursion formulae to account for layered media (Ward & Hohmann 1987). Since only the inductively coupled parts are considered, only the recursion formula of the TE mode is used. We compute sensitivities with regard to layer resistivities or thicknesses semi-analytically using direct differentiation in Wait's recursion formula.

Since individual electromagnetic field components are inverted, instrumental effects on the measured data are more severe than in methods that use transfer functions of fields (such as the AMT and CSAMT methods). The periodicity and finite duration of the turn-on and turn-off processes of the transmitter signal have substantial effects on the measured transients (Asten 1987; Fitterman & Anderson 1987). A smaller period (i.e. inverse of the repetition rate) of the current form results in the transient decaying more rapidly at late times, whereas a longer duration of the turn-off ramp diminishes the induced voltage at early times. Furthermore, the receiver coil and receiver filter impulse responses reduce the high-frequency components of transients, predominantly affecting the induced voltages at early times (Effersø et al.1999). A comprehensive description of the different instrumental effects that need to be accounted for in the system response is given by Christiansen et al. (2011). We incorporate all instrumental effects by convolving the TEM fields and their sensitivities with the respective system responses in the time domain.

APPENDIX B: IMPLEMENTATION OF DEPTH CONSTRAINTS

In our hybrid inversion scheme with additional depth constraints, the layered model section is parametrized as $${\bf m}_{L}=(\log \rho _{1},\ldots ,\log \rho _{N_{L}},\, \log t_{1},\ldots ,\log t_{N_{L}-1})^{T}$$ for NL layers. The additional depth constraint term $$Q_m^z$$ (eq. 13) is quadratic in the secondary parameter (layer depth) but non-quadratic in the primary parameter (log t). Hence, for the constrained interfaces, Taylor expansion of the modelled depths around the primary model parameters of the previous iteration k needs to be applied (e.g. Auken & Christiansen 2004), that is,  

(B1)
\begin{eqnarray} \bf {z}_{k+1} & = & \bf {z}_{k} + \bf {J}^z_k\left(\bf {m}_{k+1}-\bf {m}_{k}\right), \end{eqnarray}
where  
\begin{eqnarray*} \bf {J}^z_k & = & \left(\frac{\partial {}z_i}{\partial {}m_j}\right)_ {{i=1,\ldots ,N_{ic}\atop j=1,\ldots ,N_m}}, \nonumber \\ \frac{\partial {}z_i}{\partial {}m_j} & = & \left\lbrace \begin{array}{l@{\quad}l} \frac{\partial {}z_i}{\partial {}t_j} \frac{\partial {}t_j}{\partial {}m_j} = \delta _j^{1,\ldots ,l(i)-1} t_j & {\rm for }\, m_j = \log {}t_j,\nonumber \\ 0 & {\rm for\, all\, other\, types\, of}\\&\quad {\rm model\, parameters.} \end{array}\right. \end{eqnarray*}
Substituting eq. (B1) into eq. (13), computing the gradient $$\nabla {}_{\bf {m}_{k+1}} Q_m^z$$ with respect to mk + 1, adding $$\gamma \nabla {}_{\bf {m}_{k+1}} Q_m^z$$ to the gradient of eq. (14) in Kalscheuer et al. (2012) and equating the sum to zero yields the inversion model of the (k + 1)th iteration  
(B2)
\begin{eqnarray} \bf {m}_{k+1}(\lambda ,\gamma ) & = &\big [\left(\bf {W_{d}}\bf {J}_k\right)^T\bf {W}_d\bf {J}_k \!+\! \gamma \left(\bf {W}^{z}_{m}\bf {J}^z_k\right)^T\bf {W}^{z}_{m}\bf {J}^z_k \!+\! \lambda \bf {W}_{m}^T\bf {W}_m\big ]^{-1} \nonumber \\ & & \left(\left(\bf {W}_d\bf {J}_k\right)^T\bf {W}_d\hat{\bf {d}}_k+\gamma \left(\bf {W}^{z}_{m}\bf {J}^z_k\right)^T\bf {W}^{z}_{m}\hat{\bf {z}}_k\right) + \bf {m}_r,\nonumber \\ \end{eqnarray}
where  
\begin{eqnarray*} \hat{\bf {z}}_k & = & \bf {z}_r-\bf {z}_k + \bf {J}^z_k\left(\bf {m}_{k}-\bf {m}_{r}\right), \nonumber \\ \hat{\bf {d}}_k & = & \bf {d} - \bf {F}\left[\bf {m}_k\right] + \bf {J}_k\left(\bf {m}_k-\bf {m}_r\right). \end{eqnarray*}
In our notation, d is a vector of Nd electromagnetic field data, the inverse standard deviations of the field data are stored in a diagonal weighting matrix Wd, F[mk] is a vector with electromagnetic forward responses for the model of the kth iteration, Jk is the Jacobian matrix computed with respect to the model mk of the kth iteration, mr is a reference model for layer and distortion parameters (with entries of layer parameters equalling those of mk and entries for the distortion parameters set to zero to enforce minimum length solutions), and the matrix  
(B3)
\begin{equation} \bf {W}_m = \left( \begin{array}{r@{\quad}r}\bf {I}_{N_L} & \bf {0} \\ \bf {0} & \beta \bf {I}_{N_{PQ}}\\ \end{array} \right) \end{equation}
contains the Marquardt–Levenberg regularization operators $$\bf {I}_{N_L}$$ for layer parameters and minimum solution length regularization operators $$\bf {I}_{N_{PQ}}$$ for distortion parameters. The matrices $$\bf {I}_{N_L}$$ and $$\bf {I}_{N_{PQ}}$$ are appropriately sized identity matrices that are weighted by a manually chosen trade-off parameter β. The Lagrange multiplier λ is determined in each iteration with a line search, whereas γ is assigned a fixed value. As the implementation of prior depth information in eq. (13) mimics that of an additional data set, the weighted Jacobian matrices and data difference vectors of electromagnetic data can be conveniently augmented with the respective quantities of the a priori depths while applying the inversion.

As a measure of fit to the electromagnetic field data, the normalized misfit $${\rm RMS} = \sqrt{\frac{1}{N_d}\left\Vert \bf {W}_d\left(\bf {d}-\bf {F}\left[\bf {m}_{k+1}\right]\right)\right\Vert ^{2}}$$ is computed for every model.