Large-scale 3D inversion of semi-airborne electromagnetic data — Topography and induced polarization effects in a graphite exploration scenario

Electromagnetic (EM) methods are sensitive to the electrical conductivity of the subsurface and, therefore, play an essential role in the characterization of mineral deposits. Controlled-source EM (CSEM) methods in the time domain, e.g., Strack (1992), Haroon (2016), and Cai et al. (2022), or frequency domain, e.g., Gehrmann et al. (2019), Malovichko et al. (2019), and Védrine et al. (2023), are a widely used tool to investigate subsurface resistivities down to several kilometers. Using grounded high-power transmitters down to several kilometers depth, high-quality data can also be obtained in urbanized areas (Streich et al., 2013; Mörbe et al., 2020).

In the past few decades, semi-airborne EM (semi-AEM) systems have been developed to combine the advantage of high signal strength and, therefore, high penetration depth using grounded transmitters with fast areal data coverage of helicopter- and drone-towed receiver systems (Seigel, 1971; Elliott, 1998; Mogi et al., 2009; Wu et al., 2019). Within the project deep EM sounding for mineral exploration (DESMEX), methods for the efficient exploration of deposits at depths down to 1 km are developed. Due to advancements in sensors systems (Stolz et al., 2022; Thiede et al., 2022) and processing routines (Becken et al., 2020; Mörbe et al., 2020), high-quality data in a broad frequency range are obtained. The feasibility of the DESMEX system is demonstrated in earlier studies by Smirnova et al. (2019, 2020), Becken et al. (2020, 2022), Steuer et al. (2020), and Kotowski et al. (2022).

Due to elongated transmitters and the resulting large footprint of the semi-AEM technique, a multidimensional evaluation of the data set is, in most cases, inevitable (Cai et al., 2022). Although the advancements in high-performance computing and 3D CSEM inversion algorithms allow improved imaging of real-world scenarios, a 3D inversion of multiple transmitters with several thousands of receiver stations over a large frequency range remains computationally costly. The 3D inversion codes in the frequency domain comprise finite-difference or finite-element algorithms using a secondary field (e.g., Grayver et al., 2013; Cai et al., 2021) or total field approach (e.g., Schwarzbach and Haber, 2013; Wang et al., 2018). Smirnova et al. (2019) use the finite-difference secondary field approach presented by Grayver et al. (2013) for 3D inversion of semi-AEM data, which does not include topography. If strong topography is present in the survey area, the resulting model can be biased and the interpretation might be incorrect. The finite-element open-source toolbox customizable electromagnetic modeling (custEM) (Rochlitz et al., 2019) uses a total E-field approach and considers the topography as well as the accurate source geometry in the forward operator. Most recently, Rochlitz et al. (2023) present an inversion framework by adding explicit sensitivity calculations to custEM and combining the forward operator with the open-source inversion framework Python library for inversion and modeling in geophysics (pyGIMLi) (Rücker et al., 2017).

The presence of highly polarizable material such as graphite (Pelton et al., 1978) can cause strong induced polarization (IP) effects that superpose EM induction. IP is widely used in the characterization of mineral deposits as an additional indicator besides electric conductivity, especially for disseminated mineral resources (Revil et al., 2022). However, IP can significantly bias the EM inversion results and lead to false models when neglected in the interpretation. There are several studies investigating the effect of IP on EM data in the time and frequency domains for ground and (semi)airborne EM data (Nabighian and Macnae, 1991; Macnae and Hine, 2016; Mörbe et al., 2020). Interpretation schemes commonly use 1D approaches, either by deconvolution of IP effects (Kratzer and Macnae, 2012) or by highly constrained 1D inversion (Seidel and Tezkan, 2017). For multidimensional subsurface scenarios, several 3D complex resistivity forward-modeling algorithms and various modeling studies exist (Qi et al., 2019; Porté et al., 2023). Recent interpretation schemes include the decoupling of EM and IP effects in the time domain, e.g., as proposed by Kang and Oldenburg (2018) for grounded source data. A subtraction of IP effects from EM data or vice versa is often not feasible for complex 3D resistivity models because the IP effect on the EM transfer function is nonlinear and can affect a broad frequency range. Cox et al. (2023) propose a full inversion of time-domain AEM data, including IP effects. However, a multidimensional inversion of EM data superposed by IP remains challenging due to severe equivalence issues. The fitting of a complex resistivity will lead to a strongly underdetermined and highly nonunique problem. Adding strong constraints to the inversion is often not practicable because the distribution of IP parameters, especially for a 3D case, is not known in advance.

Here, we present a large-scale multitransmitter multicomponent 3D semi-AEM exploration study in a graphite mining district in southern Germany. The study has three key objectives: (1) derivation of a robust large-scale 3D subsurface model with resolution down to several hundred meters depth, (2) investigation of topography effects and testing the feasibility of comparably low-cost 2D interpretation schemes, and (3) subsequent investigation of IP effects on semi-AEM data.

First, we introduce the geologic setting, the survey design, and the applied processing routines. Subsequently, we present the derived multidimensional large-scale 3D inversion result and an elaborated interpretation. Because algorithms for imaging real-world scenarios such as custEM are still not routinely applied and require high-performance computing facilities, we provide possible pathways for alternative lower-cost interpretation schemes in an elaborate “Discussion” section. We first analyze the influence of topography and demonstrate the applicability of using low-cost 2D interpretation schemes. In the next step, we test model robustness by using an alternative 3D inversion code. The presented cross validation with different 2D and 3D inversion schemes also serves as a validation of model reliability. In the final step, we study the effects of IP affecting the field data and how they map into the derived 3D conductivity model.

Industrial graphite mining activities in the Kropfmühl area in the Bavarian Mountains, Germany, started in 1870 and are still ongoing. Although graphite is a critical raw material, nowadays, there exist only a few known economically feasible graphite deposits throughout Europe. The deposit in Kropfmühl is the only active graphite mining area in Germany with an average carbon content in the graphite ores ranging between 20% and 35%. The depth and the spatial extent of the deposit are only partly known from the drillings. We conducted large-scale semi-AEM field experiments over the eastern extent of the Kropfmühl graphite deposit and covered an area of $15km×8km$ (see Figure 1). Measurements were carried out on five flight days, with seven flights in total.

Geologically, the deposit is part of the Moldanubikum in the southern Bohemian Massif. The graphite enrichment occurs in the Variegated Group (Bunte Serie), which consists of highly metamorphic Precambrian paragneises (Krüger et al., 2017). The graphite-bearing garnet-(plagioclase)-biotite gneisses are interbedded with bright quartzitic gneisses, quartzitic insertions, amphibolite layers, calc-silicate marbles, and calc-silicate rocks (Krüger et al., 2017). Around the Kropfmühl deposit, the graphite is oriented from east to west (see Figure 1c). The seams were folded and exhibited extensive small-scale tectonic faulting. Graphite enrichment is often found along the 5°–15° west-dipping fault axis of the striped paragneises. Within the Variegated Series, graphite can occur locally in up to 30 layers with thicknesses in the range of a few millimeters to centimeters interbedded in the biotit gneisses or up to decimeter- to meter-thick seams in a marble-rich series. After regional metamorphosis and faulting, Variscian granite intrusions disrupted the Variegated Group locally (Teuscher, 1981).

Borelog information from the Variegated Group close to the Kropfmühl deposit down to a depth of approximately 200 m shows a highly resistive (⁠$ρ>500Ωm$⁠), nonpolarizable material of several meters thickness interbedded with a highly conductive, highly polarizable material with charge abilities up to 500 mV/V (Schulz, 1983). The conducted direct current (DC)-IP measurements in the survey area (Schiebel et al., 2023) indicate shallow chargeable material in the central part of the survey area.

On the ground, seven horizontal electrical dipole transmitters with lengths between 1.5 and 5 km were deployed (see Figure 1a). Each of the seven flight areas covered one ground-based transmitter with offsets up to 4 km to the transmitter dipole with a total area coverage of $8km×4km$ per flight area. The distance between the flight lines was 250 m, with tie lines approximately every kilometer. For higher data coverage and to obtain reliable data below the transmitters, flight areas were overlapping. The direction of transmitters 1–6 was chosen parallel to the regional strike direction of graphite in the survey area. To have dense data coverage in the central part, we set up an additional north–south-directed transmitter in the east (Tx 7). As source signal, a rectangular current function with base frequencies of 5.95 and 11.9 Hz and current strengths between 6 and 16 A was used. The high-power fast-switching transmitters, GGT-30 and GGT-10 from Zonge International, were used. Both are powered by a 400-cycle alternator connected to a ZMG generator. The internal switch box of the GGT-30 regularizes the pulse width of the transmitted signal. The continuous current was recorded with a current clamp with an edge frequency of 10 kHz and a Metronix ADU-07e logger with the same settings as for the semi-AEM acquisition. Magnetic field sensors for data recording are installed in a helicopter-towed bird approximately 40 m below the aircraft. Typical sensor heights above the ground were between 40 and 100 m. A sensor triple of obliquely arranged induction coils (Metronix, MFS-11) is installed in the bird for continuous magnetic field recording. As a data logger, an ADU-07e (Metronix) with a sampling frequency of 16,384 Hz was used. We synchronized the current and semi-AEM data recordings via the global navigation satellite system times. In a preprocessing step, the flight attitude, measured by an inertial navigation system (iNAT-RQT-4003), was used for motional noise compensation in airborne magnetic field recordings. For more technical details about the helicopter-towed system and the preprocessing, we refer to Becken et al. (2020).

Figure 2a shows the color-coded amplitude of the magnetic field components for an evaluation frequency of 1024 Hz, exemplarily for transmitter 2. Figure 2b shows the transfer functions for all components, exemplarily for the relatively central flight line 10 as an absolute value of the real and imaginary parts. Sign reversals occur along the local minima and are indicated in Figure 2b. The displayed errors consist of statistical errors from processing. In addition, an absolute error floor of 0.5 pT/A was set, and a relative error floor of 5% was added. In general, measured amplitudes decrease with increasing offsets to the transmitter. In the case of a homogeneous conductivity distribution, some of the typical features of the magnetic field for a dipole transmitter oriented along the x-axis consist of the following. (1) The $Bx$-component is zero along the centered axis perpendicular and inline to the transmitter. This behavior can be observed in the field data set in Figure 2a. The weak amplitudes in comparison to the strong horizontal component $By$ can also be observed along flight line 10, as shown in Figure 2b. (2) The strong horizontal magnetic field component $By$ decays slower (⁠$∝1/r3$⁠) than the $Bz$-component (⁠$∝1/r4$⁠) in the far field of the transmitter (Kaufman and Keller, 1983) and exhibits, therefore, higher amplitudes than the $Bz$-component at large offsets. The results of this decay pattern can also be observed in Figure 2a and 2b. Amplitudes of the $Bz$-component are lower in the far field of the transmitter (see Figure 2a). The vertical magnetic field component exhibits a stronger decay of amplitude over distance and, therefore, shows a higher error than the $By$-component for large offsets (see Figure 2b). The strong horizontal magnetic field component $By$ exhibits the slowest decay of amplitudes with offset and has, therefore, the highest signal-to-noise ratio. (3) The $Bz$-component inline and directly along the transmitter line is zero for a homogeneous or layered conductivity distribution. A drop in the amplitudes of $Bz$ inline to the transmitter can be observed in Figure 2a. In Figure 2b, absolute values of the real and imaginary parts of $Bz$ reach a local minimum when crossing the transmitter line. We do not expect a value of zero because the measurements were not taken directly at the transmitter line but were spatially averaged over 50 m.

In addition to the observation of some of the typical behaviors of the three magnetic field components, we can observe indications for a 3D conductivity distribution in the field data. One indication is the asymmetric behavior of the transfer functions in Figure 2b on both sides of the transmitter. For example, the strong drop of amplitudes for the real and imaginary parts of the vertical magnetic field $Bz$ around profile meter $y=−2000m$ followed by a plateau indicates the presence of a conductive anomaly. Fields are channeled in the anomalous structure, leading to a weakening of the $Bz$-component over the conductor.

For the 3D inversion, we used the recently developed open-source inversion tools based on custEM and pyGIMLi (Rochlitz et al., 2023). For automated mesh generation, custEM uses the software package Tetgen (Hang, 2015). For forward modeling, highly accurate second-order finite-element forward solutions are computed (Rochlitz et al., 2019). Here, we use the implemented total electric field formulation and calculation of explicit sensitivities in combination with a fast-converging Gauss-Newton minimization algorithm (Rücker et al., 2017).

For the final 3D inversion model, $Bx$⁠, $By$⁠, and $Bz$ complex-valued magnetic field-to-current transfer functions from all seven flight areas were inverted jointly. The topography was interpolated from a 25 m grid digital elevation model. The used transmitter-receiver geometry is shown in Appendix  A. To minimize the computational costs, only every fourth station was considered, leading to a station distance of approximately 200 m and an approximately even spatial distribution in the x- and y-directions. In total, the data from seven transmitters and 3093 stations are inverted jointly.

To reduce artifacts due to very strong primary fields in the vicinity of transmitters (Nazari et al., 2023), we excluded receiver stations with an offset smaller than 300 m from the corresponding transmitter line. Due to the overlapping flight areas, the data gap under the transmitter can be closed by data from adjacent transmitters. The generated inversion mesh consists of 107,860 model cells.

In total, data at 15 frequencies in a range between 5.95 Hz and 1024 Hz were used, distributed to five log-spaced frequencies per decade. The input data contain real and imaginary parts of the calculated transfer functions from all three magnetic field components. Overall, 138,649 data points were included in the inversion. The statistical errors achieved from processing were used as an error model. In addition, an absolute error floor of 0.5 pT/A was set, and a relative error floor of 5% was added. Data points with a relative error larger than 80% were excluded from the inversion.

As the starting model for the inversion, a homogeneous half-space of $500Ωm$ was used. In Table 1, the input and inversion parameters for the final 3D inversion model are summarized. As a measure of the data fit, we use the global error-weighted root-mean-square value (called $χ$⁠). Note that an $χ$ value of one means that the data are explained by the model within its error bounds. With a final global $χ$ of 3.3, the data are underfitted. Several error models were applied for the inversion, including increased absolute error floors, which resulted in a lower final $χ$ value. However, these do not improve the convergence behavior during the inversion and resulted in a similar 3D model. Due to the vast extent of graphite occurrence, there are several reasons that can lead to a nonoptimal data fit. These include the impact of the IP effects and the presence of anisotropy.

The regularization parameter $λ$ serves as a trade-off between the data fit and the model smoothness. We ran the inversion using several $λ$ values. To achieve an acceptable fit and to account for the expected sharp conductivity contrasts between the resistive host rocks and the conductive graphite structures, we chose a relatively low regularization parameter $λ=1$ for the final 3D inversion run.

We show the results of all three components and one frequency along one flight line. In addition, we show the error-weighted relative difference color coded for all frequencies and the complete data set corresponding to transmitter 6 in Figure 3. Note that the data density differs between the field components. The amplitude of the $Bx$-component parallel to the x-directed transmitter is typically smaller, and data points with a high relative error level were excluded prior to the inversion. The strong horizontal component $By$ typically decreases more slowly for large offsets than $Bz$ and exhibits, therefore, the highest data density above the noise level. The color-coded misfit plots in Figure 3a indicate the areas with over- or underestimated field amplitudes. Typically, the misfit tends to be lower for larger transmitter offsets due to the assigned relatively higher error floor.

Figures 4 and 5 show the final 3D conductivity model. Overall, the resistivities range between the fixed resistivity limits of 1 and 10.000 Ωm and exhibit sharp contrasts. Figure 4a and 4b shows the conductivity model for different x-y depth slices and y-z slices. Note that the topography surface level in the survey area ranges between 400 and 750 m above sea level (m.a.s.l). The area covered by the receiver stations is outlined in black. Areas in the northeast and southwest are not covered by data and are not considered in the interpretation. An overview of the coverage of the final inversion model is shown in Appendix  A.

In the shallow subsurface, conductivity structures arise locally in the complete survey area. A general trend of especially deep-rooted conductivities (at 200 m.a.s.l.) in the east–west direction can be observed. The outlines of the dominant east–west anomaly C are marked with the dashed black line in Figure 4a.

Model cross sections running north–south at x = −1500, 0, and 1500 m are shown in Figure 4b. Anomalous high-conductivity structures can be observed along all three cross sections and down to approximately −1000 m.a.s.l.

In Figure 5, we compare the obtained 3D model with the occurrence of graphite in borelogs and apparent resistivity isolines from the airborne helicopter electromagnetic (HEM) survey (Rehli et al., 1979). The 900 Hz HEM data provide rather shallow information within the upper 100 m, depending on the conductivity. Isolines from the HEM survey marking resistivities $<100Ωm$ are highlighted in black. For comparison, we show a rather shallow surface slice at 450 m.a.s.l. of the semi-AEM 3D inversion model. Apparent resistivities drop down to $1Ωm$ at approximately y = −2000 m and agree well with the marked main anomaly C in the semi-AEM 3D model. The drilled occurrence of graphite is in good agreement with the dominant conductive structures.

With a final global $χ$ of 3.3, the data are not fitted optimally, and therefore, small-scale anomalous features should be interpreted carefully. Furthermore, graphite usually occurs in thin bands with thicknesses in the decimeters to several meters range. These bands exhibit high conductivities, and it is not possible to resolve each small graphite band individually; rather, the bulk resistivity of the zones with interbedded graphite layers is resolved.

The joint inversion of multicomponent magnetic field-to-current transfer functions from seven overlapping flight areas with high data density results in a consistent 3D subsurface model that is in good agreement with a priori geologic and geophysical information (see Figure 5).

In practice, often, high-performance computational resources are not available for 3D interpretation using modern inversion schemes, such as realistic transmitter geometries, topography, and flexible finite-element meshing. In other scenarios, data are collected solely along single far-distant flight lines, or a fast focused profile inversion is required during field measurements. For these cases, a 2D inversion can be a feasible tool for data interpretation. In the following, we provide possible pathways for alternative lower-cost interpretation schemes for interpreting large-scale semi-AEM data: (1) we discuss the effects of topography, (2) we evaluate the feasibility of using a 2D inversion scheme for semi-AEM data, and (3) subsequently, we test model robustness by using an alternative 3D inversion code without topography. Because secondary effects such as IP can superpose the recorded semi-AEM signal, especially for mineral prospects such as graphite, we investigate how these effects map into the derived 3D conductivity model discussed in Figure 5.

For real-world scenarios with rugged topography, the induced EM fields may be significantly distorted. A model can be heavily biased if the inversion algorithm only considers 2D or flat topography. Often, the shallow subsurface relating to high-frequency data is not imaged correctly.

For models with high-conductivity contrasts, the overall impact of topography might be less severe because the minimization algorithm is mainly driven by the sensitivity toward heterogeneous conductivity distributions. Therefore, we investigate, in the following, the impact of topography on the field data set.

We inverted the $Bz$ transfer functions from one selected flight line using the 2D implementation of custEM consisting of a 3D forward model and a region-based model reduction scheme. In Figure 6a and 6b, we compare the inversion model, including the topography along the flight line, with the inversion model assuming a flat surface. For the flat model, we placed the receiver points at the corresponding height above ground and projected the transmitters onto the flat surface. Differences in the resulting models are visible mainly as the increased resistivity magnitudes around the high-conductive anomaly at the profile meter 7000. The overall spatial extension of the resistivity anomalies (C1–C5) and the overall resistivity amplitudes remain mostly unaltered. Both models exhibit a similar misfit. We conclude that for field data scenarios with moderate topography and strong conductivity contrasts, the topography has a minor impact. Nevertheless, the effects can be severe for strong topography or less pronounced conductivity differences. Therefore, we recommend using an inversion algorithm wherein the topography can be included.

One major drawback of performing a full 3D inversion is the required computational resources. Even if the geometry of the transmitter is in three dimensions, an inversion for 2D conductivity can be achieved, either through transforming the 3D forward problem in the wavenumber domain (Key, 2016) or by model reduction of the inverse problem (e.g., Ronczka et al., 2020). The former approach is implemented within modeling with adaptively refined elements for 2D electromagnetics (MARE2DEM), which we refer to as the 2.5D inversion approach. Within custEM, we follow the latter strategy using a 3D triangular prism model wherein all the prisms with the same triangle are treated as one unknown, implemented in pyGIMLi (Rücker et al., 2017). We refer to this as the 2D inversion approach.

Table 2 summarizes the computational requirements using custEM to produce the final 3D model for spatial flight data (Figure 5), a 3D profile data inversion, and a 2.5D profile data inversion. The 3D inversion of the full field data set (approximately 70,000 complex data points) results in huge computational resource requirements (approximately 1 TB of memory and 7 days of central processing unit [CPU] time). Often, such great resources are not available for 3D interpretation using modern inversion schemes. Moreover, a 3D inversion is often underdetermined, and the data density might not be sufficient in the off-profile direction or not spatially evenly distributed. For these cases, a 2D inversion or a 3D inversion for single flight lines can be a feasible tool for data interpretation. Hence, we evaluate the feasibility of using such subset approaches for semi-AEM data.

In Figure 6b and 6c, we compare the 2D implementation of custEM with 2D topography, assuming noncrooked transmitters with the 3D inversion of the same single flight-line data set. The starting model was set to $500Ωm$⁠, and only the vertical transfer function was inverted. The 3D inversion of the single flight-line data exhibits an improved data fit (⁠$χ=2.8$ for the 2D inversion and $χ=1.5$ for the profile 3D inversion). Computational resources for the two inversion runs, as well as for the full 3D inversion, are provided in Table 2. The 3D inversion of the single flight-line data set assumes a 3D conductivity distribution of the subsurface and a full 3D topography, which is not accounted for in the 2D inversion. The resulting inversion models for both approaches depicted in Figure 6b and 6c exhibit a well-comparable conductivity distribution. However, we notice that the lateral location and depth of some conductive structures, e.g., around the southern transmitter, deviates between both inversion models. Furthermore, the 2D inversion exhibits higher resistivity contrasts and an overall higher resistivity range.

We conclude that a profile-wise 2D inversion for the given data set can reproduce the general structures and is, therefore, a suitable tool for the inversion of 2D distributed data or a first, fast, computational resource-saving interpretation. However, for flight lines which cross the transmitter in the central part, a joint inversion of all three magnetic components can only be performed by a 3D inversion because the $Bx$-component equals zero along the centered axis perpendicular to the transmitter for a 2D case. Furthermore, in a 2D inversion routine, transmitters must be simplified as straight transmitter lines, which is not a suitable approach for very crooked transmitters. Therefore, for a 3D conductivity distribution with arbitrary transmitter geometry, we recommend a full 3D inversion for final data interpretation.

To validate the presented 3D subsurface image in Figure 4, we compare the obtained model with results using different multidimensional inversion schemes. We consider a structure as robust if it prevails for different inversion algorithms.

For a single flight-line inversion, we compare the 3D custEM implementation with the 2D custEM implementation and MARE2DEM (Key, 2016). Note that, for all comparative studies, we focused on the vertical magnetic field component due to computational reasons. For both 2D inversion runs, a 2D topography along the flight line and noncrooked elongated transmitters are assumed. The comparison of the 2D custEM and the 2.5D MARE2DEM inversion results shown in Figures 6b and 6d exhibits mostly small-scale differences. Differences appear mainly around profile meters −3000 to −4000. Conductivity structures in the custEM models reach to the surface, whereas they reach deeper in the MARE2DEM model. All other main features are reconstructed well, and the derived models are in good agreement. The most obvious differences relate to the degree of smoothing for the MARE2DEM model, which is due to the regularization schemes. In custEM, several inversions with varying regularization parameters were performed, wherein a fixed regularization parameter of 10 gave a good trade-off between data fit and model smoothness. MARE2DEM uses an optimization line search scheme for the determination of the regularization parameter in each iteration. Because we encountered improved convergence behavior for the MARE2DEM inversion using amplitude and phase as data input, we performed a data transformation accordingly. For both algorithms, amplitudes and phases were inverted with a 5% amplitude error level and 2.8° absolute error on the phases. The starting model was set to $100Ωm$ for all inversion runs.

To further validate the robustness of the 3D inversion, we compare the 3D custEM results with the 3D finite-difference algorithm 3DInv (Grayver et al., 2013) using the complete 3D data set. For the latter, only a flat surface can be included.

In Figure 7, the results from both algorithms are compared as a slice plane view at the surface. Both depth slices show well-comparable conductivity structures. Increased differences are visible along the northernmost transmitter between −2000 and 0 m along the x-axis. In general, conductivity structures obtained with custEM seem to be more continuous than the 3DInv model result, mainly due to the rather coarse mesh used for the latter. Differences in the vicinity of the transmitter are a result of the flattened transmitter geometry in the 3DInv model. Note that small-scale differences can originate from different smoothness constraint strengths, the flat surface approximation, and the differences in mesh discretization. The discretization of the transmitters was validated in advance for a flat homogeneous half-space model for both codes and exhibited differences well below 1%, except directly over the source singularity.

Due to the good agreement of all presented inversion models in Figures 6 and 7, our study serves as cross validation between the different inversion algorithms. Although, our results cannot be generalized for all subsurface scenarios, they provide possible pathways for alternative lower-cost interpretation schemes. We showed that for the presented field scenario, the topography effects are only minor. Nevertheless, we recommend using a full 3D topographic model if the expected conductivity contrasts are low.

Depending on its mineral texture, graphite exhibits strong IP properties, and these secondary effects can superpose the inductive EM response (e.g., Mörbe, 2020; Porté et al., 2023). IP effects will result in a frequency-dependent resistivity $ρrock(ω)$⁠. The measurement of the IP parameters is often used to deliver complimentary information, especially in mineral exploration prospects. However, IP effects are difficult to handle if they superpose the EM response, and they can bias the EM inversion if not incorporated into the modeling routine (Mörbe et al., 2020).

The resistivity and IP borelog data located a few kilometers west of the survey area show a clear correlation between high conductivity and high chargeability (Schulz, 1983). Thin (⁠$<5m$⁠) conductive layers with chargeabilities as large as the resolution limit of the logging tool (500 mV/V) are interbedded in a resistive background that correlates with low-chargeability values. Only poorly resolved chargeability borehole logging data are available from within the survey area. IP values for further modeling are based on the early studies by Pelton et al. (1978), who give an overview of Cole-Cole IP parameters from measurements over different graphite deposits. These indicate chargeabilities between $m=700$ and 1000 mV/V and a relaxation exponent $c$ close to 0.25. The observed high-frequency dependence of the electrical properties of graphite is in accordance with the results of Katona et al. (2021), who investigate graphite from a former mining area in Lower Austria. In general, $τ$ is greatly material dependent, with values ranging between $10−5$ and $104s$⁠. The investigated graphite deposits by Pelton et al. (1978) exhibit very high values varying from several tenths of seconds up to several thousand seconds.

To understand the possible effect of the time constant on the EM data for the chosen field setup, we performed a synthetic 1D modeling study. For modeling, we used the 1D algorithm tipforward1d (Hoheisel et al., 2004), which calculates the EM field responses considering the frequency-dependent resistivities based on the Cole-Cole model. The effect of the varying $τ$ is displayed in Figure 8, with the underlying 1D model shown in Figure 8b. A 50 m thin conductive and chargeable layer is interbedded in a resistive nonpolarized half-space. We varied the time constant $τ$ between 0.05 and 500 s in comparison to a nonpolarized 1D model. The resulting transfer function, exemplarily shown for a station with an offset of 2000 m broadside to the 1 km extended transmitter dipole, is shown in Figure 8a. The IP leads to a drop in amplitude for the real part in a frequency range between approximately 10 and 200 Hz. The imaginary part is affected in a wider frequency range between 1 and approximately 600 Hz. We find that a larger time constant results in a larger shift of the transfer function along the frequency axes toward lower frequencies, resulting in a larger relative difference of the IP disturbed transfer functions compared with a nonpolarizable model response. However, due to the 3D conductivity and 3D IP distribution, a 1D model study is not sufficient to estimate the impact of the IP effects on the semi-AEM data set.

The EM inverse problem is highly ill-posed and, therefore, suffers from nonuniqueness. This holds for single parameter problems, but it is especially pronounced when the IP parameters are included (Seidel and Tezkan, 2017; Porté et al., 2023). A joint inverse problem for the electrical resistivity together with the Cole-Cole parameters for the derived large-scale 3D subsurface model shown in Figure 4 is highly nonunique and ill-posed. To estimate the possible 3D model reconstruction bias produced by the IP effects in the field data, we carry out a forward-modeling study considering the following steps:

1. Conductive anomalies with $ρ≤5Ωm$ shown in Figure 9a are assumed as graphite-rich structures.

2. In accordance with the borehole-logging results and the literature values given by Pelton et al. (1978), the Cole-Cole parameters (⁠$m=0.8$⁠, $c=0.25$⁠, $τ=500s$⁠) are assigned to all mesh cells with $ρ≤5Ωm$⁠. As demonstrated in the synthetic modeling study in Figure 8, larger time constants have a stronger effect on the EM response. Therefore, we performed the modeling study with a high time constant.

3. The 3D forward calculation of the field data model in Figure 4 using the exact semi-AEM survey geometry and IP parameters assigned to each cell is performed.

4. The 3D resistivity inversion of the synthetic IP-affected semi-AEM data set using the same inversion parameters and error model as for the field data set.

Figure 9b and 9c compares a depth slice at 450 m.a.s.l. for the field data model and the synthetic IP imposed inversion result. We conclude some key aspects as follows:

1. All large-scale anomalies, especially in the central modeling region, remain predominantly unaltered.

2. The magnitude of some features is altered slightly, e.g., the magnitude of the conductivity of the main structure.

In Figure 9d, the transfer functions of transmitter 2 for the observed field data and the corresponding model response are displayed against the IP distorted data set and the corresponding model response for one selected frequency along a profile crossing the main anomaly. In Figure 9e, we show both sets of transfer functions for one selected station on top of the main anomaly. A visible deviation of the IP distorted data set over the complete frequency range can be observed. The sign reversal of the IP distorted data is shifted toward lower frequencies when compared with the original field data set.

For further analysis, the distribution of the magnetic fields (⁠$Hz$⁠) is visualized in Appendix  A for the IP- and non-IP-affected 3D model in Figure 9b and 9c, showing that the anomalous field distribution is predominantly related to the conductivity distribution and only weakly related to the imposed IP effects.

Conclusively, the introduced 3D modeling study clearly demonstrates that IP effects affect the observed field data for the introduced modeling parameters. However, they do not alter the main reconstructed conductivity distribution of the presented large-scale Kropfmühl 3D model using the introduced error model. Note that this study cannot be generalized to arbitrary subsurface conductivity scenarios but serves as an orientation on how to address IP-affected data in 3D CSEM imaging.

We present a robust and well-validated large-scale 3D inversion model derived from semi-AEM exploration in the Kropfmühl graphite mining district. The recently developed open-source finite-element 3D modeling and inversion algorithm custEM was successfully applied to the large-scale semi-AEM data. In total, approximately 138,000 data points within a frequency range of 5 Hz to 1 kHz were jointly inverted. The wide frequency range, high 3D data density from seven flight areas, and complementary information from all three jointly inverted magnetic field components result in a consistent and comparatively detailed 3D inversion model. Connected high-conductive 3D subsurface structures could be imaged within several hundred meters of depth.

Reconstructed anomalous high-conductive structures show sharp contrasts compared with the surrounding high-resistive host rock. Surface anomalies correlate well with high-frequency HEM data and the known local occurrence of graphite indicated by the geologic background information. The features mainly orient from the east to west direction, where the main central conductivity anomaly reaches down to approximately 700 m depth. Although the conductivity distribution is 3D, we demonstrate that the data set can be stably interpreted profile wise using a 2D algorithm. This is especially significant if computational resources are limited or data are only collected along several profiles, as is often the case for ground-based surveys. For an elaborate data interpretation, we recommend a 3D interpretation of the data, including all three components of the magnetic field-to-current transfer functions.

Neglecting topography and IP effects in the inversion process might lead to false models and, consequently, an incorrect interpretation. Different interpretation scenarios for large-scale semi-AEM data, however, demonstrate that the effect of topography is secondary compared with the response characteristics driven by conductive anomalies. However, this conclusion holds only for the investigated survey area with moderate topography, shallow conductors, and strong conductivity contrasts. The robustness of the 3D semi-AEM model was tested and validated by running two different 3D inversion algorithms. Although differences occur, the overall structures are qualitatively well comparable. Subsequently, we included the IP effects in a realistic 3D inversion using synthetic data. Our study demonstrates that the IP effect does not alter the overall conductivity structure and, hence, the derived 3D model is reliable. Obviously, the inversion model is driven by the strong conductivity contrasts rather than the distortion produced by the IP.

We would like to thank the field team of the Kropfmühl Survey 2020 and the DESMEX working group. Furthermore, we would like to thank N. Krüger (Graphit Kropfmühl GmbH) for discussions about the regional geology and insight into the Kropfmühl graphite deposit. The presented work was financially supported by the Federal Ministry of Education and Research (BMBF) of Germany as part of the DESMEX-II and DESMEX-real projects under grant nos. 033R130CN and 033R385E. We thank the editor and two reviewers, R. Peng and A. Hördt, for their suggestions, which helped us to further improve the manuscript significantly.

Data associated with this research are available and can be obtained by contacting the corresponding author.

The coverage of the final 3D inversion model, displayed as the sum of the absolute values of all the entries of the Jacobian matrix of the final 3D inversion model, is provided in Figure A-1. As a threshold, we used an error-weighted and normalized coverage value of 1e − 6. Coverages below that value are clipped. Central areas with higher data density exhibit high coverages. Coverages are the highest in the central part of the survey area due to a high data overlap and at the grounding points of the source. The overall coverage decreases with increasing depth.

Figure A-2 shows the comparison of fields at the surface considering EM induction with and without additional IP effects. As a base model, the final inversion model of the field data is used (see Figure 4). Exemplarily displayed is the vertical magnetic field component $Hz$⁠. We show a rectangular slice plane view of the simulated data at the earth’s surface for an evaluation frequency of 10 Hz and for transmitter Tx 2. The overall field distribution is rather heterogeneous but well comparable for both models. To obtain realistic insight into the impact of the IP on the field distribution, we weighted the data set by the error floor model, which was applied to the field data set (0.5 pT/A absolute error and 5% relative error). Due to the application of an absolute error floor, the impact of the IP in the imaginary part is down weighted. Error-weighted relative differences are mainly below 10%. The anomalous field distribution is driven by the conductivity distribution and is only weakly related to the imposed IP effects. However, toward the southern part where regional IP-affected conductive structures are present, the error-weighted relative differences are increased.

Biographies and photographs of the authors are not available.