Probing the Southern African Lithosphere With Magnetotellurics—Part I: Model Construction

The Southern African Magnetotelluric Experiment (SAMTEX) involved the collection of data at over 700 sites in Archean to Proterozoic southern Africa, spanning features including the Kalahari Craton, Bushveld Complex, and voluminous kimberlites. Here, we present the first 3D inversions of the full SAMTEX data set. In this paper, we focus on assessing the robustness of the 3D models by comparing two different inversion codes, jif3D and ModEM, and two different subsets of the data, one containing all acceptable data and the other containing a smaller selection of undistorted, high‐quality data. Results show that the main conductive and resistive features are imaged by all inversions, including deep resistive features in the central Kaapvaal Craton and southern Congo Craton and a lithospheric‐scale conductor beneath the Bushveld Complex. Despite this, differences exist between the jif3D and ModEM inverse models that derive mainly from the differences in regularization between the models, with jif3D producing models that are very smooth laterally and with depth, while ModEM produces models with more discrete conductive and resistive features. Analysis of the differences between these two inversions can provide a good indication of the model resolution. More minor differences are apparent between models run with different subsets of data, with the models containing all acceptable data featuring higher wavelength conductivity variations than those run with fewer stations but also demonstrating poorer data fit.


The SAMTEX Data Set
MT is a passive electromagnetic technique to infer the resistivity of the subsurface from measurements at the surface (Chave & Jones, 2012).Together with seismic tomography and potential field methods (e.g., gravity, magnetics), MT is one of the foremost geophysical techniques to image the structure of the lithosphere-asthenosphere system.It has been used to investigate continental lithospheric structures in many regions around the world (e.g., Gatzemeier & Moorkamp, 2005;Jones, 1999;Rao et al., 2014;Selway, 2018;Wannamaker et al., 2017) and has been a component of large national programs in the United States (Kelbert, 2019), Australia (Kirkby et al., 2020), and China (Dong et al., 2013).Based on simultaneous measurements of naturally occurring variations of the horizontal components of the electric field E and the magnetic field B, we can estimate the frequency-dependent, complex-valued MT impedance Z, viz.

𝐄𝐄 = 𝐙𝐙𝐙𝐙.
(1) This estimation process is based on robust statistical methods and thus gives formal estimates of data uncertainties (e.g., Chave & Thomson, 2003) although these are typically increased before use in an inversion algorithm (e.g., Miensopust, 2017, and discussion below).
The MT response estimates used for our inversions were acquired during the SAMTEX between 2003 and 2008 (Jones et al., 2009).They comprise >700 BBMT and LMT measurements across South Africa, Botswana, and Namibia (see Figure 1 for locations).The primary goal of SAMTEX was to image the lithospheric architecture of the cratons and mobile belts in the region, and thus measurements were taken at intervals of roughly 20 km and at periods 0.01-2,000 s for the BBMT sites, and 60 km at 10-10,000 s for the LMT sites.Due to long recording 10.1029/2021JB023117 3 of 26 times and favorable noise conditions in many parts of the study area, the data quality at many sites is excellent throughout the period range.To date there have been a number of publications performing modeling and inversions for different subregions and profiles (Evans et al., 2011;Finn et al., 2015;Hamilton et al., 2006;Khoza et al., 2013;Miensopust et al., 2011;Moorkamp et al., 2019;Muller et al., 2009).In addition, maps of resistivity directly derived from the data at selected periods have been used to investigate the structure and composition of the lithosphere-asthenosphere system (Jones et al., 2012(Jones et al., , 2013) ) and multiobservable petrological-geophysical models have been created based on subsets of the data (Fullea et al., 2011).Still, to date, no three-dimensional resistivity models of the lithosphere based on all the available data have been published.
We show representative data from six stations across the array in Figure 2. Their locations are shown as yellow stars in Figure 1.We concentrate on the period range of 1-10,000 s as lithospheric-scale structures are the focus of our modeling efforts.The sites show good data quality throughout the plotted period range, although at the longest periods (>1,000 s), the scatter and error estimates increase at some of the sites (e.g., KAP45).The off-diagonal apparent resistivity sounding curves highlight a moderately conductive shallow subsurface (∼10 Ωm at short periods) at the northern sites RAK011A and ZIM117.Further south, at sites KIM428 and KAL014, the short period apparent resistivities are generally comparable but slightly higher, while the southernmost sites, KAP045 and KAP019, show significantly higher short period apparent resistivities of ≥500 Ωm.Beneath this surface layer, the sounding curves of the four northern sites suggest that resistivity increases with depth to a maximum of ∼500 Ωm at ∼500 s, before either remaining constant or decreasing slightly at longer periods.In contrast, at the two southern sites (KAP45 and KAP19), most of the apparent resistivity curves have higher initial values and decrease with period.

Inversion Algorithms
We invert the observed MT impedances with two different inversion algorithms (Kelbert et al., 2014;Moorkamp et al., 2011) and two different strategies for each code.This helps us to address model uncertainties related to algorithm-specific choices, such as regularization, error floor, model discretization, and precision of the forward modeling engine.As a result of these choices, each code will fit different aspects of the data, including those affected by noise, to different degrees.Both algorithms used here are well established and have been compared to each other and other inversion code based on synthetic data (Miensopust et al., 2013).We therefore only give a brief summary of each algorithm and focus on comparing the differences and the potential impact on the results.
For two inversions, we use the MT inversion module of the joint inversion framework jif3D (Moorkamp et al., 2011).The numerical basis of the forward modeling engine and the gradient calculation is presented in Avdeev and Avdeeva (2009) and Avdeeva et al. (2015).It utilizes an integral-equation-based forward engine x3d (Avdeev et al., 1997) and includes a correction for galvanic distortion at each site (Avdeeva et al., 2015;Moorkamp et al., 2020).Galvanic distortion of MT impedances is typically caused by charge accumulation at small structures (compared to the induction length scale), and it can mathematically be described as a site-specific, frequency-independent multiplication of the impedances with a real-valued matrix C (Chave & Jones, 2012).In jif3D distortion correction is achieved by estimating the elements of C as part of the inversion and multiplying the synthetic impedance Z synth at each site with the corresponding distortion matrix when calculating the data misfit.Details on this methodology can be found in Avdeeva et al. (2015) and similar approaches since then have been implemented by Kordy et al. (2016) and Usui et al. ( 2017), e.g.The inversion algorithm has been used on a range of MT data sets, both commercial and academic, including imaging the fault structure for an intraplate event in Botswana (Moorkamp et al., 2019), detecting hydrothermal fluids in the central Andes (Pearce et al., 2020) and investigating the influence of a plume on oceanic crust (Franz et al., 2021).Due to the integral-equation-based forward modeling algorithm, the horizontal cell sizes need to be constant in both orthogonal horizontal directions.In the vertical direction, cell sizes can vary and are typically fine near the surface, increasing by a constant factor with depth to match the decreasing resolution of MT data.The discretized region is embedded in a layered half-space that is kept constant throughout the inversion.Some care must be taken to avoid the strong influence of this background conductivity structure on the inversion results.In order to enforce positive conductivity values during optimization and restrict model conductivities to realistic values, conductivities in each model cell are transformed using the generalized model For each site, we plot apparent resistivity and phase of the four impedance elements.Off-diagonal (xy and yx) apparent resistivities are plotted with a consistent y axis to highlight differences in average resistivity between different regions, while apparent resistivity for the diagonal components is plotted with a different scale for each site for better readability.
parameter scheme described in Moorkamp et al. (2011).This allows us to use an unconstrained optimization algorithm based on a limited-memory quasi-Newton method (L-BFGS; Avdeeva & Avdeev, 2006) to minimize the objective function.Within jif3D, we regularize the inversion through a first-order approximation of the spatial gradient of the generalized model parameters.Applying the regularization to the generalized parameters instead of the resistivity values in each cell has the advantage of equalizing the vast range of Earth conductivities (∼10 −1 -10 6 Ωm) to a range between approximately −2 and 2, and ensuring that the regularization operates similarly in all parts of the model.As the regularization is purely smoothness based, it has the potential disadvantage that structures may horizontally or vertically smear into regions of low resolution, such as those with poor site coverage in the heterogeneous SAMTEX array.However, this could be considered as a form of natural interpolation.
The other two inversions were performed using ModEM, a well-established and freely available 3D MT inversion code (Egbert & Kelbert, 2012;Kelbert et al., 2014).It is widely used in the academic community and has seen applications on data sets around the world (e.g., Dong et al., 2020;Kelbert & Egbert, 2012;Meqbel et al., 2014;Robertson et al., 2020;Ye et al., 2020;Zhang et al., 2016).Its forward engine is based on a finite-difference formulation (e.g., Egbert & Kelbert, 2012;Mackie et al., 1994) and its modular structure allows for inversion of different combinations of electromagnetic data (e.g., Campanya et al., 2016).Compared to jif3D, the gridding requirements are less strict with variable-sized rectilinear cells in all three coordinate directions.Furthermore, no background layered half-space is prescribed; instead, the grid must be extensive enough that secondary electromagnetic fields are insignificant at the model boundaries.Thus, a typical strategy for designing inversion grids in ModEM is to use an inner core with constant horizontal cell size and padding cells of increasing size around it.In ModEM, the natural logarithm of conductivity is used as a model parameter which enforces positivity of conductivity and has an equalizing effect similar to the generalized model parameters in jif3D.ModEM does not allow the range of permitted conductivities to be directly limited, but the regularization limits the difference from a prior model, often the starting model, and lateral variations of conductivity simultaneously (Egbert & Kelbert, 2012).Compared to a pure smoothing-based regularization, this combined approach should reduce smearing but can result in artificial changes in conductivity if the prior model is not representative of the average conductivity in the region.In this case, poorly resolved regions of the model will be kept at prior conductivity values, whereas well-resolved regions will exhibit a different conductivity.
It would be interesting to perform additional inversions with algorithms that use a significantly different optimization approach, e.g., data-space optimization (Siripunvaraporn & Egbert, 2009;Siripunvaraporn et al., 2005) or the increasingly popular finite-element method for calculating the synthetic data (e.g., Grayver, 2015;Jahandari & Farquharson, 2017;Kordy et al., 2016) in order to further assess the influence of these factors.This is, however, beyond the scope of this discussion.

Data Selection
Within each inversion algorithm, we ran an "all data" inversion of the entire SAMTEX data set (with only clearly erroneous stations removed) and another "selected" inversion of a subset of the data.This approach was designed to test the impacts of heterogeneous station coverage and of noisy and distorted data on the results of each inversion algorithm.A similar strategy was pursued by Lee et al. (2020), e.g., to homogenize data coverage.
SAMTEX station coverage is highly heterogeneous compared to other large-scale initiatives such as USArray (Kelbert, 2019) or AusLAMP (e.g., Kirkby et al., 2020;Robertson et al., 2016;Thiel et al., 2020).Due to logistical constraints and the still prevalent two-dimensional inversion approaches at the time of planning the measurements, data were collected in relatively dense transects separated by significant gaps.In a 3D regional model, the crustal structure will therefore be strongly represented in the data near those profiles and completely absent in regions without coverage.In contrast, deeper features (50-200 km) will at least be partially sensed by the data even in regions without direct station coverage.For this reason, the focus of our inversions will largely be on the regional imaging of the mantle lithosphere-asthenosphere system.To image the mantle, dense sampling along the profiles could be either beneficial or detrimental.On the one hand, dense coverage should result in redundant information and thus reduce the influence of noise for deep imaging, but on the other hand, dense measurements can be highly affected by local structures that cannot be represented well in the regional model.This issue might be further exacerbated by the need to choose a global regularization parameter for the model, as localized structures in densely covered areas might require a small regularization parameter.However, small regularization parameters might be inappropriate for regions without dense coverage.
The noise levels of SAMTEX data are also heterogeneous (Figure 3).Some sites show significant noise across the whole period range, with either highly scattered or physically unrealistic data (e.g., ELG010A).These sites were excluded from all inversions.Of the remaining sites, some (e.g., BOT405) display smooth sounding curves with phases in quadrant but also demonstrate large offsets between the apparent resistivity curves, indicating local static distortion.Others (e.g., KAP047 and WIN011) show similar signs of static distortion and additionally display rapidly varying phases that extend out of the quadrant, which could indicate local noise or strong resistivity contrasts in the shallow subsurface.Even though jif3D can correct for distortion and strategies have 10.1029/2021JB023117 7 of 26 been devised for ModEM to mimic the effects of distortion (Meqbel et al., 2014), it is unclear to what degree the information from these sites is useful to constrain deep structures and whether fitting distorted sites prevents fitting other data.To investigate the impact of such sites, they were therefore retained in the "all data" inversion but excluded from the "selected" inversion.To further reduce station density and to assess the impact of heterogeneous station coverage, additional stations with low-quality long-period data (>500-10,000 s) were also removed.After the selection procedure, the resulting "selected" data set has a station spacing along the station transects of ∼30-80 km in most regions, compared to ∼20 km for the "all data" data set (compare blue and red dots in Figure 1).

Inversion Setup
To be able to accommodate the entire region in a single model with acceptable computational run times, the horizontal discretization for the core region was chosen to be 15 km in the northing and easting directions for all inversions.Including the padding cells for the runs with ModEM, the inversion domain comprises 132 × 133 × 53 cells with a vertical discretization of 50 m for the topmost cells increasing up 141 km at the bottom of the domain.
Information on ocean bathymetry was introduced from the ETOPO1 global topographic data set, and seawater was assigned a resistivity of 0.3 Ωm.On land, no topography was considered, and a starting resistivity of 100 Ωm was assigned to all cells.
Inversions for ModEM were run with error floors of 5% of  √  on all tensor elements.The starting λ was set to 10 and decreased by a factor of 5 when the inversion when root-mean-square (RMS) misfit difference is <0.002.An isotropic smoothing operator was constructed with the covariance matrix set to 0.4 in all directions.For the inversions with jif3D, we used the same error floor as for the inversions with ModEM.We removed the outer padding cells from the grid used for modeling in ModEM resulting in a mesh with 119 × 120 × 48 cells and chose a fixed background resistivity of 100 Ωm.
The inversions for jif3D were run with a similar approach to regularization as the inversions with ModEM.However, we used different values for the regularization parameter, staring with λ = 1,000 and reducing it to λ = 1 in the final iterations, since the influence of the regularization on the inversion is different between the two algorithms.The initial iterations did not include any distortion correction, but this was enabled after the first regularization change as this has been shown to yield stable results (Moorkamp et al., 2020).

Data Fit
For the selected data sets, the inversion algorithms reach a final RMS of 1.7 (jif3D) and 2.3 (ModEM) after 200 and 146 inversion iterations, respectively.For the inversion of the full data sets, the corresponding RMS values are 4.4 and 5.0, respectively.We show the final RMS misfit at each site for all frequencies in Figure 4.When we only invert the selected data (bottom row) using both jif3D and ModEM, we achieve a relatively homogeneous Figure 5. Map of distortion estimates for the inversions with jif3D in the central region of the array.We show the estimates for the inversion with all data as circles and for the selected data inversion as squares.These have been displaced north from the original locations for better visibility.Colors mark the deviation of the distortion matrix elements from the identity matrix.We highlight sites BOT405 and KAP047 shown in Figure 3 with black stars.10.1029/2021JB023117 9 of 26 RMS between 1.5 and 2.5 at the majority of sites, and only a few sites exceed RMS values of 4.5.While there are some differences in how well sites are fit, the overall pattern is comparable and some sites are fit better in one inversion or the other.In contrast, when inverting the maximum amount of data, the distribution of RMS becomes much more heterogeneous.Many sites are still in the 1.5-2.5 range, but some sites exceed RMS values of 10.This effect appears to be more pronounced for ModEM than jif3D and some sites that were fit well in the "selected" inversion are now fitted significantly worse.These observations confirm that some data that were excluded are, in fact, problematic for the inversion.At least for some of these sites, the distortion correction used Comparison between observed data (symbols) and predicted data (lines) for the inversion run with jif3D (blue and red lines) and ModEM (black and green lines) and all data for four selected sites marked in Figure 1. by jif3D helps to achieve a better fit.The question remains, though, to which degree this impacts the final models.⎞ ⎟ ⎟ ⎠ for the two jif3D inversions at the central area of the array around the sites BOT405 and KAP047 (Figure 3) that were previously identified as distorted (a version with all stations can be found in the model and data archive https:// doi.org/10.5281/zenodo.5870584).For the selected data inversion, C is close to the identity matrix at virtually all sites indicating little to no galvanic distortion.This demonstrates that the data selection process successfully removed stations with significant distortion and that the inversion algorithm does not introduce artificial distortion, e.g., to achieve a low data misfit with a smooth model.When inverting the complete data set, some but not all of the additional sites show significant distortion and sites BOT405 and KAP047 are among the most distorted (Figure 5).
In theory, if galvanic distortion is caused by structures that are small compared to the typical induction scale length at short periods (Chave & Jones, 2012) and the distortion correction only represents this structure-related distortion, the estimates of C at neighboring sites should show little correlation.Although this is the case in some regions, we also see clusters of sites with very similar distortion estimates, e.g., south of site BOT405.Here, the estimate of C xx is consistently smaller than unity and C yy larger than unity at most sites.There are two possible explanations for this phenomenon: (a) It is possible that these sites were all installed in similar geological conditions, e.g., when looking for softer ground in an environment dominated by outcropping bedrock.(b) More likely, the distortion estimates capture variability in structures that can, in principle, be resolved by MT measurements but cannot be represented by the chosen horizontal discretization of 15 km, i.e., they account for the so-called model discrepancy (Kennedy & O'Hagan, 2001).These stations are located at the northern end of the Kaapvaal Craton crossing into the Magondi Mobile Belt, and thus it is likely that significant deformation is recorded in the crust.
We compare the data fit for the distorted site BOT405 and the exemplary sites RAK011A, KIM428, and KAP045 for the ModEM "all data" inversion and the jif3D "all data" inversion (Figure 6).The difference between observed and predicted data for site BOT405 clearly shows how the distortion correction in jif3D helps to achieve a better fit to the off-diagonal apparent resistivity curves.Whereas the model response from ModEM converges to a common apparent resistivity value at short periods, jif3D reproduces the constant offset between the two curves.Interestingly, although the off-diagonal phases are fit differently by both inversions, there is no clearly superior fit by either of the two models.At sites RAK011A and KIM428, both models produce virtually identical responses for the off-diagonal apparent resistivities and phases and match the observed data well.At site KIM428, the models reproduce all variations of the curves, while at site RAK011A, the overall shape is reproduced well by the models, but the phase anomaly in the xy-component at periods between 50 and 100 s is not fully reproduced by either model.At both stations, the diagonal elements are significantly smaller than the off-diagonal elements and are matched better by the jif3D inversion than the ModEM inversion.It is our experience that distortion correction helps to match diagonal elements better even when these are small (Moorkamp et al., 2020).The data at site KAP045 are matched differently by the two inversions, and again the difference is more pronounced in the phase than in apparent resistivity.The response from ModEM reproduces the short period phases well but shows small but consistent differences in the overall shape.In comparison, jif3D appears to reproduce aspects of the general shape of the yx-component better but does match the phases exactly in any period range and shows some discrepancies for long periods of the xy-component.
The observed differences in model fit highlight that different inversion algorithms reproduce different aspects of the observed data that go beyond the changes expected from simply modifying the regularization in a single inversion algorithm.This contrasting behavior illustrates the value of inverting data with multiple inversion algorithms.It also shows that distortion correction can help fit certain aspects of the data, as demonstrated by the misfit maps, but this does not necessarily imply that all aspects of the data are matched more closely.To investi- gate the degree to which distortion correction helps in these inversions, we plot the misfit at each site for ModEM versus the corresponding misfit with jif3D in Figure 7. Here, color indicates the amount of distortion, i.e., the magnitude of deviation from the identity matrix ‖C − I‖.This plot confirms that at sites with little estimated distortion, both codes achieve similarly low data misfits.All sites with large data misfit also show significant distortion and at the majority of these sites the distortion correction in jif3D helps to achieve a lower misfit.However, there are some sites where ModEM achieves a lower misfit despite significant distortion.It would be interesting to further investigate the characteristics of the data and models at the sites showing extreme behavior.However, this is beyond the scope of this manuscript, but would make for a valuable future contribution.All in all both inversions with all data match the observations at the majority of sites well.We therefore expect both models to provide reasonable representations of electrical resistivity in the vicinity of the measurement sites.

Resistivity Models
We show horizontal cross-sections through the derived inversion models between 50-km and 200-km depth in 50-km intervals  as well as vertical slices in the east-west direction at latitude 22° south (Figure 12) and along the Kaapvaal (Figure 13) and Kimberley (Figure 14) profiles (see also Figure 1 for location of these profiles).The different inversions show very similar large-scale structures, e.g., a generally resistive (≥500 Ωm) central region below 50-km depth, which is significantly more resistive than the starting model (100 Ωm).Embedded in this resistive lithosphere are several large conductors, typically associated with boundaries of different geological units.Even though the large-scale picture is similar for all models, there are significant differences in the detailed resistivity structures and values between the inversion results.We will therefore start with a description of the main features based on the "all data" jif3D model and discuss how these are expressed in the other models.In the next section, we use the differences and similarities between the models to appraise the robustness and resolution of inversion results.
We observe the maximum resistivity (≥5,000 Ωm) around the south-eastern part of the array (labeled the Kaapvaal Resistor (KR) on the horizontal slices) and the north-western part of the array, north of the Damara Conductive Belt (DCB).In both cases, the maximum resistivity is located at depths between 50 and 100 km and appears to decrease at 150-km depth and below.These observations are compatible with the thick, dry lithospheric mantle associated with the roots of the Kaapvaal Craton and the Congo Craton, respectively (e.g., Evans et al., 2011;Jones et al., 2013;Khoza et al., 2013).In the central part of the array, around latitude 24° south, is a roughly east-west striking band of reduced resistivity (∼100 Ωm) in the deeper slices (150 km and below) which becomes more resistive in the shallower parts of the model.We term the central structure in this band at ∼24° east the Molopo Farms Conductor (MFC).It can be identified as a zone of decreased resistivity (<20 Ωm) on the 150-km and 200-km depth slices from the two ModEM inversions.The jif3D-based inversions only show a weak signature at 150 km but show a structure with similarly low resistivity in a similar location at 200-km depth.The conductor associated with Bushveld Complex (Bushveld Conductor, BC) appears on all different modeling schemes north of the KR (Figure 13).Its resistivity (<10 Ωm) and horizontal location are imaged similarly in the maximum data inversions.However, it appears to be slightly wider and extend to shal-Figure 9. Horizontal slices through the inversion models at a depth of 100 km.For an explanation of abbreviations, see Figure 8. lower depth in the model produced by jif3D.Overall, its spatial extent is consistent between models and it is consistently positioned beneath the surface expression of the Bushveld intrusive complex, suggesting that it is a robustly modeled feature.
Further south, all inversions indicate a low resistivity zone (20 Ωm) at a depth of 50 km near the south-western terminus of the Kaapvaal Craton, which we term the Southern Kaapvaal Conductor (SKC).The "selected data" ModEM inversion shows this low resistivity extending to depths ≥150 km, but this is less clearly visible in the inversion with all data.Both jif3D inversions show decreased resistivities of the same magnitude as ModEM with similar differences between the selected and maximum data sets and there appears to be some degree of downward smearing.At 200-km depth, the expression of the SKC is stronger in the jif3D models than in the ModEM models.We will discuss this in more detail below.
In the north-western part of the array, the signature of the DCB, previously identified by Khoza et al. (2013), is apparent at a depth of 50 km in all inversions.It is an east-west striking band of decreased resistivity (∼10 Ωm), interpreted to be associated with the collision between the Congo Craton and the adjacent mobile belts.At depths of >100 km, the inversions with all data also contain an approximately north-south striking, low resistivity feature.The inversions with selected data also show slightly decreased resistivity in the same region, but it appears that some information on this feature is contained in the sites excluded in the selection process.Furthermore, its expression is stronger in the jif3D-based inversion than in the results obtained from ModEM.The jif3D model indicates two parallel north-south striking features while the ModEM model only contains one.It has to be noted that the western of the two conductive bands in the jif3D model coincides with a region of no data coverage.We therefore only consider the eastern band, which is expressed very similar in the ModEM model, as robust feature and suspect that the western band is, at least to some degree, an effect of regularization.
In addition to these four features discussed above, the model contains a variety of other structures.We do not go into further detail on all these features here but in the second part of this study (Özaydın et al., 2021a, 2021b) we investigate the relationships between the geoelectric lithospheric architecture, composition, tectonic, and magmatic history of southern Africa in detail.
The vertical slices through the model shown in Figures 12-14 confirm the inferences made by comparing the horizontal slices, showing similar low and high resistivity features.However, the exact locations, shapes, and resistivity values vary between the different inversions.In all cases, the "all data" inversions show stronger resistivity contrasts and more localized features than the "selected data" inversions for the same inversion algorithm, particularly in the upper 50 km.Below this depth, the differences between using all data and selected data are less pronounced, but persist.For example, the ModEM "all data" slice along the KAP line (Figure 13) shows a low resistivity zone at a depth of 200 km toward the northern (right) end of the profile.This feature is not clear in the "selected data" inversion, which shows resistivity values comparable with the starting model, possibly suggesting that there is little resolution in this region.Comparing the results from the two different inversion algorithms, ModEM appears to favor more concentrated features at depth while in jif3D the features are generally more distributed with less sharp edges.2011) along the KAP line (see Figure 1) which crosses the central Kaapvaal Craton.On the large scale, the 2D results show strong similarity with the new 3D models presented here.The SKC, KR, and BC are image in the same locations as in the 3D results and with roughly comparable geometry.When comparing details, e.g., the exact shapes or resistivity values, we can identify significant differences.For example, the BC is imaged as a near vertical feature in the 3D inversions, but appears to extend further south-west in the 2D results.Similarly, the dip of the SKC is outside the range of variation indicated by the four 3D models.In contrast, the KR exhibits a similar depth extent and shape as the 3D inversions.

Model Appraisal
To provide a quantitative view on the differences between the models, we plot model difference matrices at a depth of 50 km in Figure 15 and 150 km in Figure 16.In both cases, we show a horizontal slice through each model at the respective depth on the diagonal.Plots above the diagonal show the difference in logarithmic resistivity for the different model combinations, while plots below the diagonal show the corresponding resistivity difference histograms.The model histograms show significant differences in resistivity between all model combinations of up to 2 orders of magnitude (2 in logarithmic units), even though for the vast majority of model cells, the difference is <±1 order of magnitude.The histograms appear to be slightly wider at 50-km depth than at 150-km depth.Most histograms are centered around a difference of zero, suggesting that there is no significant overall bias in the resistivities retrieved in each inversion, except the histogram for the two inversions with selected data, which is centered around ≈0.2, indicating that the model produced by jif3D is consistently more resistive.
The histograms clearly illustrate that the largest resistivity differences are produced by using different algorithms to invert the same data set, while smaller differences are produced by inverting different subsets of data with the same inversion algorithm.At both depths and for both inversion algorithms, the histograms comparing the "selected" and "all data" inversions show highly symmetric shapes and a concentrated peak at zero, while the other histograms are generally broader and exhibit more structure.
The spatial difference plots in Figures 15 and 16 add more detail to the global resistivity differences displayed in the histograms.Spatial comparisons between the jif3D and ModEM models using both "selected data" and "all data" data sets at both depths consistently show that the jif3D models have higher average resistivities over much of the model space than the ModEM models, except for in the south-western part of the array where jif3D produces a consistently less resistive model.The south-western region is the part of the model most poorly constrained by station coverage.These differences sum to a resistivity difference histogram that centers on zero.In all difference plots, we also see a correlation between locations of large-scale tectonic boundaries and changes in sign of the resistivity difference, particularly along the northern margin of the Namaqua-Natal Belt and the  1).From top to bottom, the inversion runs are: jif3D with selected data, ModEM with selected data, jif3D with all data, ModEM with all data.We also show the 2D inversion result presented in Evans et al. (2011) in the bottom panel.margins of the Damara and Ghanzi-Chobe belts.While the details vary, this phenomenon is observed in all combinations of models to varying degrees.This indicates that the differences in the models are not merely due to fitting aspects of the data differently or the influence of noisy measurements, but each inversion images the Earth in a different way.
Comparisons between the two "selected data" inversions and the two "all data" inversions demonstrate that data selection has a significant impact on model differences.The spatial difference plots for the jif3D and ModEM "selected data" inversions reveal broad zones of consistent resistivity differences, while those for the two "all data" inversions show much more inhomogeneous, spatially varying resistivity differences.This results from the stronger influence of regularization in the "selected data" inversions, leading to overall smoother models.When adding data, the wavelength of the patterns decreases, and we see more fine-scaled differences, together with relatively sharp changes between positive and negative differences.Some of the largest differences are located in regions without site coverage, e.g., southeast of the KAP line or in the gaps between measurement lines in the northern part of the array.
The most likely candidate for causing many of these differences is the different regularization schemes.This factor is most clearly seen in the comparison between models produced with the "selected" data set as the influence of regularization is strongest there.Where the Earth is more resistive than the starting model in both inversions, jif3D consistently estimates higher resistivities than ModEM.Conversely, where the inversions indicate lower resistivities than the starting model, jif3D underestimates resistivity compared to ModEM on the larger scale.Both observations can be explained by the fact that ModEM minimizes the difference to the reference model and smoothness simultaneously, while jif3D only aims at recovering a smooth model.This behavior can explain the observed correlation between major tectonic boundaries and changes in sign of the resistivity difference.A resistive geological region is likely to be modeled with a higher resistivity in jif3D than ModEM, and an adjacent conductive geological region is likely to be modeled with a lower resistivity in jif3D than ModEM.The model difference plot, therefore, highlights the boundary between these two regions.Without additional information, we cannot say which of the inversions is more representative of the true resistivity within the Earth.However, we tried to reduce the effect of regularization in ModEM by running an additional inversion with a starting and reference model based on laterally smoothed apparent resistivities.To construct this model, for each measurement site, we construct a circle with a 4-degree radius centered on the site and take the median apparent resistivity value of all sites within the circle at periods longer than 100 s.The resulting resistivity value is assigned to all cells below this site.We then perform a linear interpolation of logarithmic resistivity between these values to determine the resistivity in each model cell.The resulting model (Figure 17) shows The resulting median inversion model (Figures 18 and 19) fits the selected data to an RMS comparable with the inversion run from a homogeneous half-space.Compared to the homogeneous inversion run, the average resistivities at 50-km depth (Figure 18) and 150-km depth (Figure 19) are higher, particularly in regions that are not 10.1029/2021JB023117 20 of 26 directly covered by sites.Conductive anomalies show a very similar pattern to the previous inversions, although the shape and location differ slightly in some cases, including some of the individual conductors that form the DCB at 50-km depth (compare Figures 8 and 18) or the SKC at 150-km depth.These changes are not significant enough to imply a different geological interpretation of these structures.
The spatial difference plot and difference histogram comparing the median inversion and the inversion of the same data with jif3D reveal some interesting changes compared to the inversion with a homogeneous starting model.Visually, the spatial difference plot for the median model contains a lot less long-wavelength structure and is dominated by more small scale differences.This contrast is particularly visible at 50-km depth where jif3D previously produced consistently higher resistivities in the central model region.However, the spatial difference plot for the median model displays a much more variable pattern where the sign of the conductivity difference changes within smaller distances.At 150-km depth, the effect is less pronounced yet still observable.
The difference histogram comparing the median inversion and the jif3D inversion at 50-km depth has a maximum very close to zero, while the histogram comparing the homogeneous inversion and jif3D is offset to slightly positive values with a maximum at ≈0.2.While still not fully symmetric, the maximum and minimum (most negative) differences now show a similar magnitude.At 150-km depth, the impact of the starting model on the histogram is even more pronounced and the median model histogram is more significantly offset to negative values associated with the higher average resistivity of the median model.While the average resistivities inverted from the homogeneous ModEM are downward biased compared to jif3D, the median model ModEM resistivities are upward biased.

Discussion and Conclusions
The main goal of this paper is to present a new 3D conductivity model for southern Africa and use the different inversion methodologies to understand uncertainties in the results better.An additional result is that the detailed comparisons of the models also reveal some technical aspects of inversions and regularization that are of interest to both algorithm developers and practitioners and thus warrant some discussion before describing some of the geological interpretations implied by these models.
It is our impression that most of the differences between the results from the two inversion algorithms stem from the different regularization philosophies.The purely smoothness-based approach followed by jif3D spreads out structures to their maximum possible extent, most clearly visible in Figure 13.Selecting such a smoothing operator has the disadvantage that conductive anomalies can be smeared out, and their boundaries can be challenging to identify.In contrast, the mixed regularization approach pursued by ModEM typically produces more localized structures.On the flip side, the regularization toward a reference model appears to bias the large-scale resistivity toward this model, particularly in regions of low sensitivity.This observation is mirrored by the systematic study of Robertson et al. (2020).Taken together, we conclude that for large arrays with heterogeneous coverage such as this, jif3D produces models with more representative large-scale resistivity values, while ModEM produces more focused and localized anomalies.To some degree, a more representative large-scale resistivity can be obtained with ModEM with a median-based starting model, as shown by the comparison at 50-km depth.Still, the shift in bias at 150 km shows that possibly a more detailed starting model with varying resistivity with depth is necessary to obtain good average resistivities over large areas.Alternatively, one could design a regularization scheme where the balance between smoothness and damping toward a reference model can be finely adjusted.While it seems that such an approach could combine the advantages of both regularization approaches, it is questionable how an optimal balance could be objectively found and how practical such a scheme would be for routine application.
The effect of inverting for only selected, high-quality data, or the maximum amount of data is similar regardless of the inversion algorithm.In both cases, the inclusion of more data increases the misfit of the final models as potentially problematic sites are introduced.To some degree, this effect is reduced by the distortion correction employed in jif3D which can deal with problems associated with galvanic distortion and achieve a better fit at many sites.At the same time, the models with more data exhibit stronger resistivity contrasts and additional structures, e.g., the north-south striking conductor in the north-western part of the study area that extends from the Congo Craton into the Damara Belt.Given the similarity of these features for both inversion algorithms and the acceptable misfit for the inversion with jif3D, we conclude that these are not artifacts caused by noisy data, but that these features are due to information about the resistivity of the Earth contained in the measurements included in the "all data" models.Still, the inversions with selected data contain the same general features as the inversions with all data.Based on the similarity with other models and the data fit, our two preferred models are the ones produced by jif3D with all data and the ModEM inversion with a median starting/reference model and selected data.
The most prominent features of our two preferred models are: (a) A resistive core of the Kaapvaal Craton as indicated by the KR.region of high resistivity (>1,000 Ωm) extends to depths of 150 km (ModEM) to 200 km (jif3D) and indicates a dry lithospheric mantle in line with previous 2D interpretations (Evans et al., 2011) and experimental electrical conductivity of common mantle minerals (e.g., Karato & Wang, 2012;Özaydın & Selway, 2020).(b) Other high resistivity regions at depths of 100 km and greater include the Congo Craton in the north-west (Kamanjab Inlier) and northern Botswana in the north-eastern part of the array suggesting the presence of lithosphere with a broadly similar composition in these regions.(c) These resistors are intersected by several deep-seated conductors that are present to varying extents in all inversion models.These include the MFC, the BC, and the north-south striking feature below the Congo Craton and Damara Belt.In the latter case, a possible interpretation is a shallower lithosphere-asthenosphere boundary compared to the thick cratonic roots of the Kaapvaal and Congo Cratons (Celli et al., 2020).In contrast, the MFC and the BC are likely expressions of emplacement of metasomatic material during episodes of magmatism (Beukes et al., 2019).
A more detailed interpretation of the resistivity structures recovered by these inversions requires careful consideration of the geological history of the region and the relationships between resistivity, composition, and temperature.These considerations are beyond the scope of this study and are presented in a companion paper (Özaydin et al., 2021a, 2021b).
In summary, we have constructed a set of 3D models for southern Africa based on two subsets of the SAMTEX data set and utilizing two independent inversion algorithms.Despite some differences in the shape of structures and the recovered resistivities, the models show strong similarities.Previous efforts using these data either used the whole data set but did not perform inversions or were concentrated on regional subsets of the data.Thus, the models presented here are the first large-scale resistivity models of the region and can serve as a resource for further investigations and integration with other observations such as gravity and seismology.

Figure 1 .
Figure 1.Map of the study area.We show the Southern African Magnetotelluric Experiment (SAMTEX) measurement sites considered in the full data inversions as blue dots and the sites considered in the inversions with selected data as red dots.Yellow stars indicate exemplary stations for different regions shown in Figure 2 and red stars poor quality data excluded from some of the inversions shown in Figure 3.The red, green, and yellow lines mark the locations of vertical model profiles along the Kimberley, Kaapvaal, and 22° south lines, respectively.Black lines mark the boundaries of tectonic provinces based on McCourt et al. (2013).

Figure 2 .
Figure2.Six exemplary sites representing different regions within the inversion domain.For each site, we plot apparent resistivity and phase of the four impedance elements.Off-diagonal (xy and yx) apparent resistivities are plotted with a consistent y axis to highlight differences in average resistivity between different regions, while apparent resistivity for the diagonal components is plotted with a different scale for each site for better readability.

Figure 3 .
Figure 3.Examples of data excluded from the inversions with selected data.Site ELG10A shows overall problematic data and has been excluded from all inversions while the other sites show potentially problematic features as discussed in the text but have been retained for the inversions with maximum data.

Figure 4 .
Figure 4. Map of misfit for the four inversion runs.We show the error normalized root-mean-square (RMS) across all frequencies at each site.

Figure 6 .
Figure 6.Comparison between observed data (symbols) and predicted data (lines) for the inversion run with jif3D (blue and red lines) and ModEM (black and green lines) and all data for four selected sites marked in Figure1.

Figure 7 .
Figure 7. Misfit at each site for the two inversions with maximum data.We plot the misfit for the inversion with ModEM versus the misfit for jif3D.Color indicates the total magnitude of distortion and the black line indicates equal misfit in both inversions.

Figure 10 .
Figure10.Horizontal slices through the inversion models at a depth of 150 km.For an explanation of abbreviations, see Figure8.

Figure 11 .
Figure 11.Horizontal slices through the inversion models at a depth of 200 km.For an explanation of abbreviations, see Figure 8.

Figure 13
Figure13also shows a comparison with the 2D model ofEvans et al. (2011) along the KAP line (see Figure1) which crosses the central Kaapvaal Craton.On the large scale, the 2D results show strong similarity with the new 3D models presented here.The SKC, KR, and BC are image in the same locations as in the 3D results and with roughly comparable geometry.When comparing details, e.g., the exact shapes or resistivity values, we can identify significant differences.For example, the BC is imaged as a near vertical feature in the 3D inversions, but appears to extend further south-west in the 2D results.Similarly, the dip of the SKC is outside the range of variation indicated by the four 3D models.In contrast, the KR exhibits a similar depth extent and shape as the 3D inversions.

Figure 12 .
Figure 12.East-west slice through the four inversion models at 22° southern latitude (yellow line in Figure 1).

Figure 13 .
Figure13.Vertical model slices through the four inversion models along the KAP line (green line in Figure1).From top to bottom, the inversion runs are: jif3D with selected data, ModEM with selected data, jif3D with all data, ModEM with all data.We also show the 2D inversion result presented inEvans et al. (2011) in the bottom panel.

Figure 14 .
Figure14.Vertical model slices through the four inversion models along the Kimberley line (red line in Figure1).From top to bottom, the inversion runs are: jif3D with selected data, ModEM with selected data, jif3D with all data, ModEM with all data.

Figure 15 .
Figure 15.Difference matrix for the four inversion runs at a depth of 50 km.We plot the resistivity slice for each model on the diagonal.Plots above the diagonal show the difference in logarithmic resistivity between pairs of models as labeled above each column, plots below the diagonal show the corresponding histogram.

Figure 16 .
Figure 16.Difference matrix for the four inversion runs at a depth of 150 km.We plot the resistivity slice for each model on the diagonal.Plots above the diagonal show the difference in logarithmic resistivity between pairs of models as labeled above each column, plots below the diagonal show the corresponding histogram.

Figure 17 .
Figure 17.Starting model derived from median apparent resistivity for the inversion with ModEM.Note the reduced range of resistivities in the color bar compared to the other model plots.

Figure 18 .
Figure 18.Difference matrix for the ModEM inversion runs with a homogeneous starting model and a median apparent resistivity starting model at a depth of 50 km.

Figure 19 .
Figure 19.Difference matrix for the ModEM inversion runs with a homogeneous starting model and a median apparent resistivity starting model at a depth of 150 km.