## 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.

About 6 km^{3} 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 km^{3} 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 km^{2} 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 km^{2} 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.

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 km^{2}, 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 (*v _{p}* ∼ 4500 m s

^{−1}) from the overlying sedimentary units with low seismic

*P*-wave velocities (

*v*∼ 1800 m s

_{p}^{−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.

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:

*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

*E*and

_{x}*E*(in V m

_{y}^{−1}) and the horizontal and vertical magnetic fluxes

*B*,

_{x}*B*and

_{y}*B*(in nT) are measured. Subsequently, the magnetic fluxes are converted to horizontal and vertical magnetic fields

_{z}*H*,

_{x}*H*and

_{y}*H*(in A m

_{z}^{−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):

*Z*= −

_{xy}*Z*and

_{yx}*Z*=

_{xx}*Z*= 0 V A

_{yy}^{−1}. To provide a more obvious link to the resistivity distribution of the subsurface, impedance tensor elements

*Z*can be transformed to apparent resistivities ρ

_{ij}_{a, ij}and phases φ

_{ij}through:

The vertical magnetic transfer function tensor **T** = [*T _{x}*

*T*] relates the vertical and horizontal components of the magnetic field as:

_{y}*T*=

_{x}*T*= 0, because electromagnetic fields entering the Earth at AMT frequencies are refracted towards the normal, independent of the angle of incidence.

_{y}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 **Z**_{0} and **T**_{0} 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:

**I**is the identity matrix,

*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 ∂*b _{z}*/∂

*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

*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 *N _{L}* 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

*N*≤

_{ic}*N*− 1 layer interfaces be close to

_{L}*N*

_{ic}*a priori*or reference values

**z**

_{r}, the additional term takes the form:

*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.

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 | 2 | 1.0 | 1.0 | 1.3 | 1.0 | 1.0 |

3 | 4 | 1.0 | 1.0 | 1.5 | 1.0 | 1.0 |

5 | 12 | 1.2 | 1.3 | 1.6 | 1.3 | 1.3 |

7 | 6 | 1.1 | 1.2 | 1.7 | 1.2 | 1.2 |

9 | 8 | 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 | 2 | 1.0 | 1.0 | 1.3 | 1.0 | 1.0 |

3 | 4 | 1.0 | 1.0 | 1.5 | 1.0 | 1.0 |

5 | 12 | 1.2 | 1.3 | 1.6 | 1.3 | 1.3 |

7 | 6 | 1.1 | 1.2 | 1.7 | 1.2 | 1.2 |

9 | 8 | 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 *Z _{xy}* and

*Z*that only involve sets of scalar distortion parameters

_{yx}*P*,

_{xx}*P*,

_{yy}*Q*and

_{xy}*Q*. Although

_{yx}*Q*and

_{xx}*Q*may deviate from zero when inverting the off-diagonal

_{yy}**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 *Z _{xy}* and

*Z*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

_{yx}*et al.*(2013b).

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).

### 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 *Z _{yx}* and

*Z*(i.e. impedance tensor elements related to the

_{yy}*E*field) are almost equal and about half that of

_{y}*Z*(Fig. 6b). Two observations suggest, that energy leaked from

_{xy}*Z*into

_{yx}*Z*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

_{yy}*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).

In inverting the CSAMT data, the strong source overprint in *Z _{yx}* and

*Z*results in

_{yy}*P*and

_{yx}*P*values of about 0.5 and −0.6, respectively (Fig. 6a). As expected (Qian 1994; Kalscheuer

_{yy}*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

*Z*and

_{xx}*T*, 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).

_{y}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).

### 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.

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.

### 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.

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 *t*_{3} becomes substantially larger than *t*_{2}.

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.

### 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.

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.

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 *t*_{2} can vary from 0.56 · *t*_{2} to 1.40 · *t*_{2}. 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.

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 |

t_{1} | 0.97 | 1.04 | 0.99 | 1.01 |

t_{2} | 0.56 | 1.40 | 0.97 | 1.03 |

t_{3} | 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 |

t_{1} | 0.97 | 1.04 | 0.99 | 1.01 |

t_{2} | 0.56 | 1.40 | 0.97 | 1.03 |

t_{3} | 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).

### 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:

Whereas electrical resistivity is only dependent on ϕ

_{i}, velocity is dependent on ϕ_{t}. Estimated porosities are consistent only if ϕ_{i}< ϕ_{t}.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.

(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 | |||||

v [ms_{p}^{− 1}] | v [ms_{p}^{− 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 | |||||

v [ms_{p}^{− 1}] | v [ms_{p}^{− 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 *v _{p}* = 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 (*v _{p}* ∼ 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 *v _{p}* = 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):

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.

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).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.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 km^{2}.

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

### 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 *H _{z}* 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 *N _{L}* 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,

**m**

_{k + 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

**d**is a vector of

*N*electromagnetic field data, the inverse standard deviations of the field data are stored in a diagonal weighting matrix

_{d}**W**

_{d},

**F**[

**m**

_{k}] is a vector with electromagnetic forward responses for the model of the

*k*th iteration,

**J**

_{k}is the Jacobian matrix computed with respect to the model

**m**

_{k}of the

*k*th iteration,

**m**

_{r}is a reference model for layer and distortion parameters (with entries of layer parameters equalling those of

**m**

_{k}and entries for the distortion parameters set to zero to enforce minimum length solutions), and the matrix

*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.