Exploration of a Copper Ore Deposit in Elbistan/Turkey Using 2D Inversion of the Time-Domain Induced Polarization Data by Using Unstructured Mesh

We present the results of a direct current (DC) resistivity and time-domain induced polarization (TDIP) survey exploring a copper ore deposit in Elbistan/Turkey. The ore deposit is elongated below a valley and is of disseminated form with sulfide content. DC and IP data were acquired using the pole-dipole array on eight parallel profiles crossing the valley perpendicularly. The length of each profile was 300 m with an inter-profile distance of about 50 m. The data were interpreted by a newly developed 2D DC/TDIP inversion algorithm. The finite element algorithm uses a local smoothness constrained regularization on unstructured meshes. The finite element forward solution, as well as the inverse problem, is solved by an iterative preconditioned conjugate solver. The depth of investigation (DOI) was determined from cumulative sensitivities of the 2D inversion algorithm results. Because of the dissemination of the ore, the 2D inversion of the DC data was ambiguous: However, due to the sulfide content, a strong chargeability anomaly associated with the ore body was detected. We show that chargeability anomalies can generally be detected in the absence or presence of corresponding resistivity anomalies. This highly chargeable structure was confined in lateral direction. Although the lower boundary of the structure could not be resolved by the applied field set-up, a rough estimation of it could be derived at a depth of 90 m using synthetic modeling analyses. The 2D chargeability models are consistent with existing borehole information.


Introduction
Mineral exploration is an important task that serves to investigate the characteristics, such as size and shape, of a mineral deposit prior to the exploitation. The application of geophysical methods serves the non-destructive characterization of the subsurface. If applied, covering a large area, it is a powerful tool to define outlines of a potential mineral deposit and suggests promising locations for boreholes or excavations. In mineral exploration, the direct current (DC) resistivity method is commonly combined with the induced polarization (IP) method (Pelton et al., 1978). The advantage of using IP in addition to pure DC is that the resistivity signature caused by disseminated material is hard to detect whereas the chargeability signature resolved by the IP method might be strong and rather independent from the geometry of the resistivity structure. In this case, the DC method gives an overview of the general geology while the IP method provides information about a potential ore deposit (Spitzer & Chouteau, 2003).
The IP method is strongly associated with the DC method. It makes use of polarization effects in a lowfrequency range that are analogous to a capacitor like reversible storage of energy . In frequency-domain (FD), mostly a phase shift between transmitter signal and measured voltage is investigated, which is dependent on the transmitter frequency. A special variation of FDIP is the spectral induced polarization (SIP) where several frequencies are employed. Time-domain (TD) applications make use of the voltage decay after transmitter turn-off. Besides information about the subsurface resistivity distribution, the IP method offers the determination of further properties. The most important property is the chargeability which is related to the ability of the subsurface to store energy. Several IP models exist that describe the instant voltage drop and the voltage decay curve after transmitter turn-off.
In the model according to Seigel (1959), the effect of the presence of chargeable material results in an increase of effective resistivity. Other models are e.g. the Cole-Cole model (Pelton, 1977) or the model by Dias (1968), just to name a few. They describe a complex frequency-dependent resistivity. Further parameters such as relaxation time and frequency exponent are introduced allowing to draw conclusions about the texture of the polarizable material. They are represented by an equivalent circuit including a capacitor. Dias (2000) provides an overview of the most commonly used IP models. Effective 2D inversion algorithms for time domain IP data were developed to take into account the whole IP transient (Fiandaca et al., 2013;Hönig & Tezkan, 2007). Besides mineral exploration (Pelton et al., 1978;Seigel et al., 1997;Vanhala & Peltoniemi, 1992), the IP method is also employed for environmental issues, e.g. characterization of contamination of the subsurface Yuval & Oldenburg, 1996), hydrogeological aspects (Hördt et al., 2007;Kemna et al., 2002), investigation of porous rocks and sediments in field surveys (Kemna, 2000;Titov et al., 2004), or laboratory experiments (Bairlein et al., 2014;Zimmermann et al., 2008). A review about recent developments of the DC and the IP method is given in Loke et al. (2013), Revil et al. (2012) and Kemna et al. (2012).
The basis of each inversion algorithm is the implemented forward calculation which serves to create synthetic data for a given subsurface model. In TDIP, commonly the model according to Seigel (1959) is employed, e.g. in the algorithms by Oldenburg and Li (1994), Loke and Barker (1996) or Karaoulis et al. (2013). Then, the synthetic IP datum ''apparent chargeability'' is derived from DC forward calculation.
In general, most geophysical forward calculations make use of the finite difference (FD) approach (Dey & Morrison, 1979;Mufti, 1976) or the finite element (FE) method (Coggon, 1971). Besides the work of Coggon (1971), further basics of using the FE method for DC problems are presented by Pridmore et al. (1981), Queralt et al. (1991) and Li and Spitzer (2002). The forward problem is to be solved under consideration of boundary conditions that restrict the potential at the boundaries of the computational domain, such as the Neumann boundary condition or mixed boundary conditions (Dey & Morrison, 1979). In FE and FD methods, the computational area is discretized into a grid consisting of mesh elements and nodes. The introduction of boundary elements (Okabe, 1981;Queralt et al., 1991) at the outer mesh boundaries facilitates the incorporation of boundary conditions. Every DC forward problem is eventually reduced to a sparse linear system of equations that is to be solved with respect to the electrical potential.
An overview over efficient solvers is given in Spitzer and Wurmstich (1999). While the FD approach is usually applied on structured rectangular grids, the FE method is suited for the use of unstructured triangular meshes in 2D (Ö zyıldırım et al., 2020) or tetrahedral grids in 3D . The advantage of unstructured meshes over structured grids is that local refinement is achieved where accuracy is needed without increasing the number of mesh elements significantly.
Oldenburg and Li (1994) proposed three TDIP inversion approaches. All of them are based on the IP perturbation model according to Seigel (1959), and they recover a resistivity and a chargeability model after two inversion steps. The first step is a DC data inversion recovering a resistivity model. The second step serves to recover a chargeability model. In the first method, the IP data is linearized, and a linear inverse problem is solved under the assumption that chargeability is small. In the second approach, a chargeability model is approximated after two resistivity inversions of perturbed resistivity models. The third method involves a non-linear chargeability inversion.
Several years ago, a copper ore deposit was discovered in Elbistan/Turkey (e.g. Ö zyildirim et al., 2017;2020). The ore deposit is assumed to be elongated below the valley. A DC/IP survey was planned and afterwards carried out to derive the lateral and, if possible, vertical borders of the ore body. The DC/ TDIP data were observed on eight profiles perpendicular to the strike direction of the valley. The data should be interpreted with 2D inversion algorithms.
For this study, we developed a new 2D DC/TDIP inversion algorithm (Adrian, 2017) and used it for the 2D inversion of the observed DC/TDIP data from the survey area in Elbistan/Turkey. The algorithm applies a local smoothness constraint regularization approach for 2D inversion of DC and TDIP data including error weighting. The basis of the inversion algorithm is a FE DC forward calculation using unstructured triangular meshes. The inversion procedures are also conducted on an unstructured mesh enabling the incorporation of surface topography. The IP chargeability is determined by a slightly modified version of the third method that was proposed by Oldenburg and Li (1994). It is a non-linear IP inversion subsequent to a non-linear DC inversion. In this approach, there is no need for the chargeability to be small. Therefore, strong and weak chargeability anomalies can be reproduced by the inversion. For more detailed information about the algorithm the reader is referred to Adrian (2017), where modeling studies with synthetic data are also included which test the applicability of the algorithm, and the results are compared with the existing conventional DC/IP algorithms introduced by Günther (2004),  and Loke (2002Loke ( , 2016.
In the following, we introduce the essential features of the applied inversion algorithm and present the inversion results of the DC/TDIP field data. The resulting resistivity models show high resistive structures below low resistive valley sediments, but from the resistivity models alone the presence of an ore body cannot be derived. The chargeability models contain, however, strong IP anomalies with a chargeability of partly more than 100 mV/V, which indicates the presence of the expected disseminated copper ore. From the inversion models, the upper boundary of this structure is clearly seen, but the lower boundary is not revealed by the inversion.
In order to investigate whether the introduction of a lower boundary of the structure followed by low chargeability is realistic, we present a forward modeling study of synthetic data based on an alteration of the chargeability inversion models. It is shown that the datafit can be improved by a lower boundary and a chargeability decrease.

Geological Background of the Survey Area
The survey area, close to the town of Elbistan, is located in the province of Kahramanmaras in southeastern Turkey. It lies at an altitude of 1170 m to 1200 m above mean sea level in the East Anatolian Fault Zone (Fig. 1) including the Tauride thrust belt where several ophiolitic structures occur (Karaoglan et al., 2013). In general, ophiolitic structures appear in connection with ore deposits (Press & Siever, 2003).
During the collision of the Arabian platform and Anatolia, an oceanic arc of the Neotethys Ocean was partly subsided below the continental crust or was accreted. The accreted material forms the South-East Anatolian Ophiolites which are cropping out along the thrust belt (Kozlu et al., 2014). Since the Göksun Ophiolite ( Fig. 1), which lies close to the survey area, belongs to the same tectonic system, it is expected that similar rocks and ore paragenesis are present as they are found at the Isbendere Ophiolite. Near the Isbendere Ophiolite, several sulfide ores were detected in veins and in disseminated form including copper and iron mineralization such as chalcopyrite (CuFeS 2 ), bornite (Cu 5 FeS 4 ), chalcosine (Cu 2 S), covelline (CuS), pyrite (FeS 2 ), sphalerit (ZnS), and marcasite (FeS 2 ) (Yildirim et al., 2012). The ophiolite is mainly non-metamorph and contains ultramaficmafic cumulates, isotropic gabbro, a sheeted dyke complex, plagiogranite and volcanic units. Extrusive mafic rocks are exposed e.g. in the Nurhak-Elbistan area. They are comprised of basalt, basaltic andesite, andesite, dacite and rhyolite. Furthermore, volcanic sediments and intercalcations of limestone and mudstone occur.
Generally, the ores are typically found in the oceanic crust where they are formed in hydrothermal deposits at oceanic arcs (Press & Siever, 2003). Where the oceanic crust is subducted beneath the continental crust, ore deposits are found in accreted fragments of the oceanic crust (ophiolites), or they originate from melting and upward intrusion of the subsiding material (Press & Siever, 2003). The deposited ores are then found along veins or are disseminated (finely distributed) within the fractured, intruded material (Press & Siever, 2003;Yildirim et al., 2012).
The survey area covers a copper ore deposit with sulphide being mainly of disseminated form which was discovered several years ago. According to information based on boreholes, the hostrock consists mainly of dacite and copper ore is expected at a depth of more than 8 m below the surface. The ore deposit is assumed to be elongated beneath a valley structure which is oriented approximately from north to south and is flanked by gentle hills (Figs. 2, 3). The valley and the hills are covered by agricultural fields with a dry surface. On the hills, several outcrops of mineralized rocks were observed.

Figure 1
Major tectonic features and distribution of the Neotethyan ophiolites (green) in the eastern Mediterranean region; after Dilek and Flower (2003), Karaoglan et al. (2013). The Survey area is indicated by a red circle, the Göksun and the Isbendere Ophiolite by a red cross, respectively. b, c Sketch of the development of SSZ-type ophiolites, granitoid intrusions and volcanic sediments, related to the arc, that were subducted or accreted during the late Cretaceous in the northern part of the southern Neotethys (Karaoglan et al., 2013). The sketched subduction is oriented from southeast (SE) to northwest (NW)

Data Acquisition and Measured Data
DC and TDIP data were acquired on eight parallel profiles perpendicular to the valley and the assumed strike direction is as shown in Fig. 2. The total length of each profile was 300 m while the inter-profile distance was around 50 m. A Pole-Dipole configuration was applied with 8 potential electrodes (P1 to P8). The whole setup was moved along the profile. The procedure was done twice, once for a dipole length of a = 10 m and once for a = 20 m. Since eight potential electrodes were used, the acquired data contains seven n-levels for each dipole length with n-spacing of n = 1 to n = 7 (Fig. 4). The largest array spread was 160 m using a dipole length of a = 20 m and an n-spacing of n = 7. Data acquisition was performed using the equipment of Zonge Engineering (Zonge, 2003), namely the transmitter GGT-3 and the receiver GDP32-II with 8 channels. The transmitting signal has a 50% duty cycle and a cycle frequency of f = 0.125 Hz. The measured apparent chargeability is an integrated chargeability given in Newmont Standard (Swift, 1973). The device provides the measured apparent chargeability data in mV/V as an approximation of the Newmont Standard according to (Zonge, 2003): Þdt, with T = 1/f, t 1 ¼ 451:38ms, t 2 ¼ 1097:32ms, the measured primary voltage V p , secondary voltage V s (t), and the internal Zonge constant s zonge ¼ 1:53 (Zonge, 2003). s swift ¼ 1:87 is the swift contant (Swift Jr, 1973).
Data errors were estimated as a sum of a geometry error, voltage error and an additional noise floor of 3%. The geometry error DG G is the inaccuracy of the geometry factor G used for the calculation of the apparent resistivity and it is based on an electrode displacement which we assume to be Dx±0.15 m. The voltage error DU U is based on the inaccuracy of the voltage measurement itself which we estimated to be DU = ± 0.5 mV. Both, geometry and voltage error are dependent on the electrode spread and amount to 2-6% and 0-3%, respectively. Exemplary for all profiles, the measured data of Profile 3 (Figs. 2, 3) are shown in Fig. 5a,b. The ordinate indicates the n-level in the form a/n. The first n-level refers to a = 10 m and n = 1, which is denoted as 10/1, and the last n-level refers to a = 20 m and n = 7, which is denoted as 20/7. The other 7 profiles show similar pseudo-sections. The measured apparent resistivity values (Fig. 5a) are mainly between 10 and 300 Xm. The minimum values are found in the first n-levels in the center of the profile and in the last n-levels at the end of the profile. The pseudo-section of the measured apparent resistivities give no clear indication of an ore body. The low resistivities in the center of the valley indicate mainly the near surface valley sediments. On the other hand, the measured apparent chargeabilities (Fig. 5b) are low close to the surface. With growing pseudo-depth, the chargeability increases building two maxima in the beginning and in the center of the profile with values of around 100 mV/V, which can be associated with the ore body (Tezkan & Adrian, 2017).

2D Inversion of the Observed DC/IP Data
The data interpretation was conducted using the 2D DC/TDIP inversion algorithm of Adrian (2017). A brief introduction of this algorithm will be given in the following. A detailed description of the algorithm, including the forward calculation and the comparison of the algorithm using analytical solutions, can be found in (Adrian, 2017). The algorithm applies a smoothness constraint regularization approach with explicit local regularization. For an inversion mesh with N p elements, the model Þ T with M data-points and the vector containing data errors is ¼ e 1 ; 1 ; . . .; e M ð Þ T . These vectors contain either resistivity or chargeability parameters. The complete functional W for iteration k reads The left norm in Eq. (1) refers to the data functional W d . It contains the error weighting matrix W d with diagonal entries 1/e m , the Jacobian or sensitivity matrix, and the model update vector Dp. J is the Jacobian or sensitivity matrix with entries j i;j ¼ À Á and the nth model parameter p k n . The data residual vector Dd = f(p)-d meas contains the difference between measured data d meas and the synthetic data according to the model p. Since Gaussian noise is expected, a L2-norm is applied. The second norm in Eq. (1) is the model functional (or stabilizing functional) W p (Tikhonov & Arsenin, 1977). It regulates the roughness of the model by applying a smoothness constrained matrix C on the model update Dp. Since the smoothness constrained is applied on the model update instead of the model itself, sharp boundaries would be allowed when they are included in the initial model. The smoothness constrained matrix C is a discrete approximation of the partial differential operators . That means it regulates changes of model parameters Dp between neighboring elements. For an element i with neighbouring elements n j (j = 1,2,3), the roughness, which is the inverse of the smoothness, is The simplest case is to set the diagonal entries C i,j = 3, the entries C i,nj = -1, and the remaining equal to zero. Like this, the contribution of all element boundaries is considered equally. The constraint matrix of this approach is denoted as C a . To consider different contributions of element vertices according to their varying length and varying element size, the length of the mutual sides L nj and the circumference of the element L all are inserted instead. Then, the entries are C i,nj = L nj , C i,I = -L all and 0 else. The corresponding constraint matrix is denoted as C b .
The functional in Eq.
(1) has to me minimised with respect to Dp which leads to In every iteration k, Eq. (3) is solved using a preconditioned conjugate gradient solver (Freund & Nachtigal, 1996). The result is the model update Dp k . The data-fit v is determined by All forward and inverse calculations are conducted on unstructured meshes facilitating the incorporation of topography.

Two-Step IP Inversion
As commonly applied in TDIP, the induced polarization is incorporated according to the model of Seigel (1959). In the presence of a non-zero chargeability m, the intrinsic conductivity r is reduced to an effective conductivity r 0 : Following this simple relation, synthetic apparent chargeability data m a can be calculated by double DC forward calculation f dc of the intrinsic and the effective conductivity as: This approach requires an additional forward calculation, but since the computational mesh is the same, the element matrices only need to be calculated once. Thus, the extra computational effort is rather small.
The used algorithm applies a two-step IP inversion approach as proposed by Oldenburg and Li (1994). It assumes that the resistivity model q' according to the effective conductivity r' = 1/q' is recovered by the resistivity inversion in the first step. The chargeability model m is obtained by a nonlinear chargeability inversion in the second step. In this second step, a resistivity model q = 1/r corresponding to the intrinsic conductivity, is also derived using Eq. (5). This resistivity model q is required for the forward calculation according to Eq. (6).
The Jacobian used in the chargeability inversion step can be obtained by simply scaling the DCsensitivity (Oldenburg & Li, 1994). The partial derivative of the data-point d i = m a,i with respect to the chargeability of the j-th element yields The potential U 0 is produced by the effective conductivity r' and U by the intrinsic conductivity r (cf. Equation 6). Under this consideration, J dc ij ¼ Different approaches of smoothness constraint matrices were introduced and compared in the 2D inversion algorithm of Adrian (2017). Reasonable results are obtained by the matrix denoted as C a that is the general form of the discrete approximation of first order derivatives. However, the employment of the matrix denoted as C b is recommended because it includes the consideration of the varying length of triangle sides and element areas in the applied unstructured mesh. Moreover, it is shown that the application of a stepwise cooling produces useful results. This procedure starts with the creation of a smooth model by using a rather large regularization parameter in the first inversion iterations. The regularization parameter is only decreased when the data-fit is not decreasing further. This procedure is continued until either the convergence criterion is met or a predefined minimum regularization parameter is reached. Moreover, it was shown that the application of an inexact line-search for an optimum model update step-length in the form of the two-point parabola technique serves to improve convergence and reduce the required number of iterations (Adrian, 2017).

Coverage and Depth of Investigation
In general, the coverage cov j or cumulated sensitivity of an element j is obtained by summation of the j-th column of the Jacobian. Obtaining the coverage from summation over W d J corresponds to the effect of error weighting (e.g. Martin, 2009). Normalizing the sum with the overall maximum coverage cov max ensures that values of different measurements are comparable because the maximum normalized coverage is always equal to one. Since the sensitivity is dependent on element size, it is reasonable to normalize the coverage additionally with the area of the element A j yielding The normalization ensures that the coverage of different elements is comparable and independent of discretization and element size. Visualizing the normalized coverage as colour-coded map reveals the sensitivity of the dataset regarding each element's parameter (resistivity or chargeability). Since the sensitivity is influenced by conductivity, the coverage also reacts to conductivity changes. We propose a method to derive a coverage threshold cov t that approximates a depth of investigation DOI app. Martin (2009) used the coverage to investigate the DOI of transient electromagnetic measurements and applied predefined coverage thresholds to separate well-resolved model domains from poorly resolved domains. Yogeshwar (2014) and Yogeshwar and Tezkan (2017) expanded this method by choosing a threshold based on the coverage of a homogeneous initial model.
For the DC case, we use the DOI according to a rule of thumb DOI t to derive the threshold cov t . This estimation is based on sensitivity studies of homogeneous models, and it is dependent on the array type and the array spreading (Barker, 1989). For the largest applied Pole-Dipole spread with a dipole length of a ¼ 20 m and a distance between current electrode and first receiving electrode of a Á n ¼ 20 m ð ÞÁ7 ¼ 140 m, the estimated depth according to (Barker, 1989) is Þ¼53 m: The coverage threshold cov t is derived from the normalized coverage of the initial model which is a homogeneous half-space. The coverage value that coincides with the DOI t in the center of the profile is defined as cov t . For the presented survey, the homogeneous initial model has a resistivity of q 0 = 200 Xm. In the center of the profile in the depth of DOI t = 53 m, the normalized coverage of the homogeneous half space is cov t = 10 -2.15 . For the inhomogeneous 2D inversion result (Fig. 6a,b), we use the cov t -isoline of the corresponding normalized coverage (Fig. 6c,d) as an indication of the approximate DOI app. It separates well resolved from poorly resolved model domains.

Inversion Results
All inversion results presented in this chapter are derived by the algorithm of Adrian (2017) using the following settings: smoothness matrix C b , stepwise reduction of the regularization parameter k, a Vol. 179, (2022) Exploration of a Copper Ore Deposit in Elbistan/Turkey Using 2D 2263 homogeneous initial half space model with q 0 = 200 Xm and m 0 = 10 mV/V, consideration of topography and surface distances. The resistivity and the chargeability inversion results of Profile 3 are presented in Fig. 6a and b respectively The main structures of the models are labelled as A-C. The corresponding normalised coverage is shown in Fig. 6c,d. Most prominent in the resistivity inversion model (Fig. 6a) is the low resistive structure (A) with less than 15 Xm and with a thickness of approx. 15 m close to the surface in the centre of the profile which is dipping towards deeper depths at profile distance 200 m. It is followed at depth by a high resistive structure (B2) of [ 100 Xm with a thickness of about 40 m. The near surface resistivity towards the beginning (B1) and the end of the profile (B3) is higher (100-300 Xm) than in the centre, as was expected from the measured data. The resistivity decreases below the highly resistive structures (B1-B3), which is also supported by a decrease in the measured apparent resistivity in the last n-levels.
The chargeability inversion model (Fig. 6b) shows low chargeability values close to the surface. There is a highly chargeable structure (C) below 20 m depth with a lateral extension of x = 20-200 m and more than 80 mV/V. The maximum chargeability values reach [ 100 mV/V and are limited to to x = 130-180 m. A lower boundary of the chargeable structure is not revealed by the inversion model.
From the surface mesh with similar element sizes and for the same electrode setup, the normalized coverage threshold of 10 -2.15 , corresponding to the approximate DOI app , was derived. According to the approximate DOI app , the depth of investigation of the resistivity model varies between a depth of 60 m below high resistive structures (B1, B2) and 80 m below low resistive structures to the left and right of the centre. The DOI app corresponding to the chargeability model has a similar form but, according to the sensitivity calculations, is shifted downwards by about 10 m.
The normalized differences displayed in Fig. 7 point out that the data are mostly fitted within the error bounds (within ±1). There are only two outliers visible according to the apparent chargeability. But, for the apparent resistivity, there are several poorly fitted data-points, especially at x = 230-250 m. This shows that the inversion model does not fully reflect the resistivity contrast which is seen in the measured apparent resistivity data in that part of the profile. Figure 8 shows the 2D inversion results of all eight profiles. The resistivity models (Fig. 8a) vary from profile to profile, but they have some features in common. The low resistive structure (\ 15 m) that lies close to the surface (cf. structure A) appears in all profiles more or less in the centre. It has a thickness of about 15 m. The dipping of this structure towards the bottom right is most prominent in Profile 3. In Profiles 1 and 2, too, there is a connection of the structure to rather low resistivity in deeper depths to the bottom on the right. Below this structure, the resistivity increases. The shape and size of the high resistive structure (cf. structure B2) varies from profile to profile as well as the amplitude of resistivity. All in all, the resistivity of the structure is up to 100-300 Xm. In most profiles, the structure continues from beneath the low resistive structure towards the beginning of the profiles, i.e. structures B1 and B2 connect. Towards greater depth, it is confined by low resistivities in Profiles 2-5. There is another high resistive structure (B3) of 100-300 Xm close to the surface at the end of Profiles 1-4. This is where the profiles ascend a hill. The resistivity in this part of the models decreases from Profile 4 to Profile 8 where the resistivity is only approx. 10 Xm. Profiles 7 and 8 do not ascend this hill and therefore different surface material is expected.
In all chargeability models (Fig. 8b), the chargeability is below 10 mV/V close to the surface, and a highly chargeable structure (cf. structure C) with amplitude of more than 100 mV/V appears at depths below approx. 10 m. The main part of this structure seems to be shifted along the x-axis from profile to profile. From that, the strike direction of structure C can be roughly approximated to be oriented from the beginning of Profile 1 to the end of Profile 8 (see sketched strike direction in Fig. 2). In Profile 8, the upper boundary of the structure is lower than in the other profiles and lies close to the approximate DOI. The presence of the ore deposit is not clearly revealed by the resistivity models alone which is certainly due to its disseminated form. However, the chargeability models show values of partly more than 100 mV/V below a depth of ca. 10 m below the surface. These values indicate the presence of sulfidic ore content. According to the inversion results, the lateral extension can be estimated, and the strike direction of the deposit is roughly oriented from the beginning of Profile 1 to the end of Profile 8. However, the final images do not show a lower boundary of this structure.
By the comparison of the resistivity and chargeability results with borehole information, the obtained physical parameters can be linked to the stratigraphy. Figure 9 shows the inversion models of Profile 1 together with a visualization of borehole information obtained from HKS4 (Fig. 2) which is located at profile distance 250 m. The colours of the borehole stratigraphy refer to resistivity values derived from the inversion results. The dacite layers can be associated with resistivity values of [ 200 Xm for the altered dacite, 50-100 Xm for the filic dacite and \ 50 Xm for the potassic dacite. These values are, however, in a smaller order of magnitude than the resistivity of dacite listed in (Telford et. al., 1990). A possible explanation is that the texture and composition of the present dacites leads to altered resistivities. According to the simplified geological map of the survey area (Robertson et. al. 2006), the surface layer in the survey area is covered by sediments of the Miocene-Pliocene. This fits to the low resistive surface layer in the centre of the valley. The copper ore content below a depth of 8 m exhibit increased chargeability values of partly [ 100 mV/ V. According to literature values, a content of 2-8% leads to a significant IP effect. The chargeability increases with a growing sulphide content. Furthermore, copper sulphides chalcocite, chalcopyrite and bornite, and copper produce a significant IP effect even for concentrations of 1%.

Lower Boundary of the Ore Deposit
The chargeability inversion did not reveal the lower boundary of the highly chargeable structure (C) beneath the investigated 8 profiles using the applied field set-up, at least not within the vertical Overview of the inversion results: a Resistivity models, b Chargeability models extension of the inversion mesh which is, for example, 200 m for profile 3 (Fig. 10). The corresponding calculated synthetic data is shown in Fig. 10b. The data fit for profile 3 (Fig. 7) suggests that the data is fit within the error bounds. In comparison to the original measured data (Fig. 5), the calculated data is mostly similar, which is supported by the data fit. However, the apparent chargeability of the data points at about x = 150 m in the last n-levels (see black ellipse in Fig. 11a) are higher than in the original data. This is emphasized by the relative differences between both data sets shown in Fig. 11a.
To investigate whether the introduction of a chargeability decrease in depth can improve the relative difference, the chargeability model was altered in the following way: Below a certain depth z D , we introduce a structure D for which the maximum chargeability is set to a certain limit m D (see Fig. 12 for an example). The resistivity model stays unchanged. Then, a forward calculation of the altered chargeability models is conducted, and the resulting synthetic data is displayed in terms of relative difference with respect to the measured data (cf. Fig. 11b-i). The left column refers to models in which structure D has a maximum chargeability of m D = 10 mV/V. This is similar to the values that appear close to the surface. The right column refers to models with a maximum chargeability of structure D of m D = 50 mV/V which is similar to the values that appear at depths below 50 m in the right part of the chargeability model. The boundary between structures C and D is set to z D = 70 m in sub-figures (b,c), z D = 80 m in (d,e), z D = 90 m in (f,g), and z D = 100 m in (h,i). For every altered model, synthetic data was calculated and compared to the measured data. For all tested models with a chargeability maximum of m D = 10 mV/V, the relative difference indicates that the calculated data is lower than the measured data. The relative difference is influenced from the 3rd n-level downwards. The relative difference according to the models with a chargeability limit of m D = 50 mV/V in the depth improves in the lower n-levels. The results of the boundary at a depth of z D = 90 and z D = 100 m are similar, but the boundary at a depth of 100 m influences the intermediate n-levels less. According to these model calculations we expect the boundary below structure C at a depth between 90 and 100 m which is slightly below the DOI app . Surely, this investigation is only an approximation, but it reveals that the data-fit of the lower n-levels is improved when a lower boundary with a transition to an intermediate chargeability value is introduced.

Discussion and Conclusion
We presented the results of a DC/TDIP survey on a copper ore deposit in Turkey. The resistivity model revealed a near surface low resistive structure in the centre of the profile, which implicates the presence of sediments in the trough of the valley. The high resistive structures are presumably associated with mineralized rock with dacite content. This is supported by borehole information and the observation of outcropping rocks on the hill at the end of the profile. Due to the state of the ore deposit, which is mainly disseminated, it was expected that the presence of the metallic particles does not produce a prominent resistivity anomaly. Thus, the deposit is difficult to detect with the DC resistivity method alone. However, the strong induced polarization effects support the hypothesis of a disseminated sulphide content. As expected, highly chargeable structures with more than 100 mV/V are recovered that are associated with the ore deposit. The expansion of the structures is not clearly distinguishable in the resistivity models.
For the modelling of the DC/TDIP field data, the 2D inversion algorithm developed by Adrian (2017) was applied. The algorithm consists of a novel b Figure 11 Profile 3: relative difference between measured data and the chargeability data calculated from models including structure D as shown in Fig. 12 Here, the data misfit is determined using linear parameters composition of finite element forward calculation and a smoothness constrained local regularization inversion approach including error weighting and a nonlinear IP chargeability inversion. All forward and inverse calculations are conducted on unstructured triangular meshes. The IP chargeability model confines the ore body in the lateral direction, but the lower boundary cannot be identified in the inversion models. By forward modelling of altered models, the location of the boundary is roughly estimated. According to this modelling study, the presence of a boundary at a depth slightly below 90 m improves the data-fit.

Acknowledgements
Part of this research was funded by the Federal Ministry of Education and Research (BMBF, project no: 01DL12028) of Germany and the Scientific and Technological Research Council of Turkey (TÜ Bİ-TAK, project no: 110Y343). We thank especially our colleagues at the Department of Geophysical Engineering/Ankara University, Turkey for the collaboration during the field survey and for providing the Zonge equipment.
Author Contributions JA: code development, field measurements, data processing, modelling; BT: supervision of the modelling, field survey, modelling; MEC: field survey, modelling.

Funding
Open Access funding enabled and organized by Projekt DEAL. Part of this research was funded by the Federal Ministry of Education and Research (BMBF, project no: 01DL12028) of Germany and the Scientific and Technological Research Council of Turkey (TÜ BİTAK, project no: 110Y343).

Data Availability
The authors confirm that all data and materials as well as software application support their published claims and comply with field standards. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code Availability
Not applicable.

Declarations
Conflict of Interest The authors have no conflicts of interest to declare that are relevant to the content of this article.
Ethical Approval Not applicable.
Consent to Participate Not applicable.

Consent for Publication
The content of the publication is agreed to by all authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.