Prerpints.org logo
Preprint
Article

Application of Audio-Magnetotelluric and Dual-Frequency Induced Polarization Joint Inversion in a Pb-Zn Deposit Exploration in Inner Mongolia, China

Altmetrics

Downloads

140

Views

70

Comments

0

Submitted:

02 June 2024

Posted:

03 June 2024

You are already at the latest version

Alerts
Abstract
Models of subsurface structures are important for successful deposit exploration but are challenged by the need to integrate data from different geophysical methods. In the present study, we evaluated a method of joint inversion in which audio-magnetotelluric (AMT) and dual frequency induced polarization (DFIP) data sets are inverted simultaneously to produce a consistent 2D resistivity model to show a clear image of subsurface structures. To achieve the objectives, we conducted AMT and DFIP surveys first on the same survey line within the Dongjun lead-zinc deposit in inner Mongolia by developing 31 AMT survey sites with an interstation separation of 40 m on a 1440 m survey track and operated in fifty-three frequencies in the range 1–10,400 Hz to record the resistivity distribution of subsurface to depths exceeding 800 m. The same survey setup up was applied to the DFIP method by utilizing the SQ-3C model with a pole-dipole array configuration and conducted measurements at frequencies of 4 Hz and 4/13 Hz. We first separately inverted AMT and DFIP measurements. The two-dimensional (2D) model obtained from AMT data revealed distinct low resistivity anomalies in the middle of the 2D inversion model. In contrast, the DFIP inversion model showed a high resistive body in the same region. In response to the discrepancies observed in the separate 2D inversion models, we implemented a joint inversion for both the AMT and DFIP datasets. The joint inversion resistivity model shows surficial conducting bodies and a high conductive body along the profile with relatively high Percent Frequency Effect (PEF) indicating high chargeability. The final joint inversion resistivity model clearly images the large silica alteration zone and the Pb-Zn mineralization. This study demonstrates the feasibility of a joint inversion methodology and highlights the value of integrating geophysical methods through joint inversion for enhanced characterization and exploration of lead-zinc ores.
Keywords: 
Subject: Environmental and Earth Sciences  -   Geophysics and Geology

1. Introduction

In recent years, there has been a notable surge in the advancement of geophysical exploration technology, which has significantly contributed to the field of mineral resource exploration [1,2,3]. The electromagnetic (EM) method and the induced polarization (IP) method are considered the most significant technologies developed for mineral exploration [4]. The development of electromagnetic (EM) geophysical methods in the 1950s has played a crucial role in mapping the lateral and vertical variations in subsurface resistivity. These methods, such as natural source audio magnetotelluric (AMT) and controlled source (CSAMT), have found wide-range of applications in metallic mineral exploration [5,6,7,8,9], groundwater studies [10], and geothermal system investigations [11,12,13,14]. Induced polarization (IP) also has a longstanding history in geophysics, initially employed in the field of mining geophysics to delineate and localize ore bodies [15,16]. In addition to the existence of advanced geophysical exploration techniques, it is important to note that each method possesses distinct advantages, limitations, and a variety of potential outcomes when applied for mineral exploration [17]. However, when the data from different geophysical techniques are inverted separately, the results can produce inconsistent inversion models; this can lead to difficulties in interpretation and create confusion about the real subsurface structures. The observed discrepancies in the inversion models indicate the potential benefits of employing a joint inversion approach to seek a single, consistent inversion model and reduce the ambiguities in the final model, which satisfies the physical properties of the two geophysical methods.
In a joint inversion methodology, multiple sets of data are simultaneously inverted [7,18,19,20,21]. The existing literature provides numerous examples of methodologies and practical instances of joint inversion techniques. Recently, several geophysical approaches have been merged for joint inversion in order to overcome the shortcomings of individual geophysical prospecting techniques. In this study, we carried out joint inversion of a geophysical dataset comprising natural source AMT and controlled source dual-frequency domain induced polarization DFIP data. Initially, individual data sets of AMT and DFIP were inverted separately to produce 2D inversion model to estimate resistivity structures sensitive to each specific survey. This process involved analyzing the unique characteristics that each method illuminated in the subsurface structures. The resulting models from these independent inversions were then carefully examined to identify both the similarities and differences in the resistivity profiles they revealed. Subsequently, a joint inversion of the AMT and DFIP data sets was conducted, allowing for simultaneous processing and integration of the information from both techniques. Compared to analysis of single data sets, this approach permits a more comprehensive and accurate representation of the subsurface resistivity structures, leveraging the strengths and mitigating the limitations of each individual method. We have three primary research objectives of this study 1) to develop a joint inversion method that incorporates multiple spatially overlapping geophysical AMT and DFIP data sets. 2) To prove that a joint inversion method increases the accuracy of the resulting 2D resistivity model 3) by applying joint inversion technique to produce the clear subsurface resistive and conductive structures with high resolution. Here we apply and test these methods in a case study of the Dongjun Pb-Zn-Ag deposit in the Central Asian Orogenic Belt. By conducting a comprehensive analysis of the inversion results, including sensitivity analyses and comparisons with existing geological information, we seek to validate the effectiveness of the joint inversion approach of DFIP and AMT and its potential for practical implementation in mineral exploration projects.

2. Regional Geology

The geographical area known as the Great Xing'an Range is characterized by its location across the Siberian, North China, and Pacific tectonic plates (Figure 1a)[22]. It is widely recognized as a complex assemblage of multiple microcontinental blocks, specifically referred to as the Erguna, Xing'an, Songnen, and Jiamusi blocks (Figure 1b)[23]. These blocks are delineated by distinct fault lines: the Tayuan-Xiguitu fault separating the Erguna and Xing’an blocks, the Hegenshan–Heihe fault between the Xing’an and Songnen blocks, and the Mudanjiang fault demarcating the boundary between the Songnen and Jiamusi blocks [22]. The Phanerozoic period witnessed significant tectonic changes, including the closure of the Paleo-Asian Ocean, the Mongol-Okhotsk Ocean, and the subduction of the Pacific Ocean. These processes led to the formation of the overall structural framework, which occurred through the northwest to southeast amalgamation of microcontinental blocks [23,24]. Although there is ongoing debate regarding the specific ages and processes involved in the amalgamation, geological researchers generally agree that the Xing'an block was incorporated into the Erguna block along the Tayuan-Xiguitu fault during the Early Paleozoic era. Additionally, it is widely accepted that the Songnen block fused with the composite block along the Hegenshan-Heihe fault during the late Paleozoic era [22,24]. The Paleo-Asian Ocean underwent its ultimate closure during the late Permian to the Early Triassic period along the Xilamulun-Changchun fault subsequent to which there was a period of regional extension [25,26]. The magmatic events during the Early-Middle Jurassic period were associated with the closure of the Mongol-Okhotsk Ocean towards the Xing'an Mongolia Orogenic Belt [26,27]. Additionally, the combined effects of the Mongol-Okhotsk Ocean closure and the subduction of the Paleo Pacific Oceanic plate could have played a role in the extensive magmatism and associated mineralization [28,29]. Furthermore, the Jiamusi block is widely regarded as an exotic block that underwent tectonic amalgamation with the Asian continent along the Mudanjiang fault in the Jurassic epoch [30].
The lithologies observed in the northern region of the Great Xing'an Range consist of metamorphic rocks belonging to the Paleoproterozoic Xinghuadukou Group. This group encompasses a Precambrian crystalline basement, metamorphic rocks from the Neoproterozoic Jiageda Group, and a Paleozoic cover sequence comprising clastic and carbonate rocks from the Cambrian, Ordovician, Silurian, Devonian, Carboniferous, and Permian periods. Additionally, the region exhibits volcaniclastic rocks and coal-bearing seams from the Jurassic and Cretaceous periods [31]. The emplacement of intrusive rocks was predominantly observed during the late Paleozoic and Mesozoic periods, with a minor occurrence during the early Paleozoic era. The primary occurrence of magmatic activity during the early Paleozoic era was concentrated in the Mohe, Tahe, and Duobaoshan regions of Nenjiang County [26]. Basic-ultrabasic rocks were primarily generated during the late Paleozoic era, predominantly emerging at the interfaces of geological blocks. During the late Paleozoic and Mesozoic periods, there was the formation of substantial intermediate felsic intrusive rocks [32]. The emplacement of igneous bodies occurred within a shallow crustal environment, manifesting as midhypabyssal, hypabyssal, and ultra-hypabyssal intrusions composed of felsic and intermediate materials. The primary rock formations consist of dacite porphyry, granite porphyry, quartz porphyry, and quartz monzonite porphyry.

3. Ore Deposit Geology

The Dongjun Pb-Zn-Ag deposit is located in the Hulun Buir area, 20 km northeast of the city of Erguna in the northern part of the Great Xing’an Range, within the eastern segment of the Central Asian Orogenic Belt (50° 21′ 30′′–50°23′N, 120°17′–120°23′E), which is in the centre of the Erguna Block and to the northwest of the Tayuan Xiguitu Fault.
The Dongjun deposit comprises several distinct strata, including the Tamulangou formation, Manketouebo formation, and Quaternary sediments (Figure 2). The Tamulangou formation consists of a sequence of volcanic rocks with intermediate to elemental compositions, such as andesite, basaltic andesite, andesitic tuff, sedimentary tuff, volcanic breccia, and small dacite. On the contrary, the Manketouebo formation primarily comprises intermediate to felsic volcanic and volcaniclastic rocks, including rhyolite, volcanic breccia, and volcanic agglomerate. The Quaternary Holocene is distributed throughout the river valley. It mostly comprises gray-black humic soil, fine sand, medium sand, pebbles, gravel, and other alluvial materials; humic silt and other swamp deposits; and residual slope (Figure 2). The region shows several fractures and faults that follow a pattern of NEE, NNW, and N-trending orientations. These geological features are believed to be genetically related to the Genhe Fault. The distribution of the orebodies in the Dongjun deposit is controlled by NNW and N-trending faults, which serve as secondary structures of the Genhe Fault [31].The mineralization occurrence in the Dongjun Pb-Zn-Ag deposit is closely linked to late Yanshanian intrusions in the area. The intrusive rocks found in the Dongjun deposit during the Yanshanian period primarily consist of granite porphyry. Granite porphyry and adjacent andesitic tuff, andesite, and sedimentary tuff characterize the Tamulangou formation. These host rocks hold significant importance in the formation.

3.1. Characteristics of orebodies

The primary host rocks for the orebodies consist predominantly of andesite, andesitic tuff, and sedimentary tuff belonging to the Tamulangou formation. Additionally, the contact zone between the granite porphyry and the Tamulangou formation also serves as a significant location for the occurrence of these orebodies. Crypto-explosive breccia bodies are observed in various zones within the Dongjun Pb-Zn-Ag deposit. The crypto-explosive breccia bodies show a prevalent morphology characterized by either a spherical or dumbbell shape. These bodies display a distinct zoning pattern, with the central region consisting of hydrothermal crypto-explosive breccia surrounded by a broader envelope of fragmented breccias. The cementing material in the crypto-explosive breccias consists of silica, andesitic tuff powder, or sulfides. The breccias primarily comprise andesitic tuff, andesite, and fragments of sedimentary tuff in the upper portion of the breccia formations. At deeper levels, the andesite fragments demonstrate significant silicic, potassic, and polymetallic mineralization. The clasts within the crypto-explosive breccia exhibit a notable variation in size and shape. Specifically, the clasts located at the periphery of the breccia bodies are characterized by their large size and angular morphology. In contrast, the clasts found in the breccia bodies' central region are relatively smaller and exhibit sub-angular or sub-rounded shapes. The breccia bodies exhibit intrusive relationships with the andesite, andesitic tuff, and sedimentary tuff units within the Tamulangou formation. Several granite porphyry dikes were observed within the underlying crypto-explosive breccia bodies, indicating a potential genetic association between the granite porphyry, the magmatic crypto-explosion, and the Pb-Zn-Ag mineralization [31].

3.2. Mineralization and alteration

Most host rocks in the Tamulangou formation, such as andesite, andesitic tuff, sedimentary tuff, and granite porphyry, have undergone varying degrees of alteration. The wall-rock alteration can be classified into three zones based on their proximity to the granite porphyry: the potassic-silicic-sericitic alteration zone, the phyllic alteration zone, and the propylitic alteration zone. The potassic-silicic-sericitic alteration zone can be distinguished because it has silicic, alkalic, sericitic, and carbonate alterations. The phyllic alteration zone is characterized by sericitic, silicic, and carbonate alterations. Substantial chlorite, epidote, and carbonate minerals characterize the propylitic alteration zone. The magnitude of alteration typically diminishes as the distance from the underlying granite porphyry increases. The Tamulangou formation exhibits the most extensive mineralization with high levels of silicic alteration [34].
Preprints 108150 g003
Figure 3. Geological map of the study area along with AMT and DFIP survey lines.
Figure 3. Geological map of the study area along with AMT and DFIP survey lines.
Preprints 108150 g003

4. Methodology and Data Acquisition

The careful selection of appropriate techniques is of utmost importance prior to conducting an extensive geophysical survey plan in a prospective area [35]. To find out if there is mineralization and its geophysical properties, as well as how it related to different rock formations, geological structures, and mineralization background, two different geophysical techniques were used to model the subsurface resistivity and predict the area of mineralization.

4.1 Audio Magnetotelluric (AMT) Method

The magnetotelluric technique for mineral exploration was initially proposed by Andrey Nikolayevich Tikhonov in 1950 [36] and subsequently refined by Louis Cagniard in 1953 [37], and underwent further advancements through the contributions of Cantwell Thomas in 1960 [38,39] and Keeva Vozoff [40]. The method probes into the sub surface’s electrical structure, employing either an artificial (controlled) or natural source electromagnetic (EM) field as the primary source of the field [41]. In the case of artificial or controlled magnetotellurics, electromagnetic signals are generated through dedicated EM transmitters. A notable illustration of artificial source magnetotellurics is the radio magnetotelluric (RMT) method, which relies on civilian and military radio transmitters operating within the frequency spectrum of 10 kHz to 1 MHz [42].
Conversely, in natural magnetotelluric exploration, the electromagnetic fields originate from global lightning phenomena, known as sferics (generating short-period signals), and solar wind activities in the ionosphere, which produce long-period signals. Natural magnetotelluric fields can be classified into audiomagnetotellurics (AMT) with frequencies within the range of f = 1–10,000 Hz and broadband magnetotellurics (BBMT) spanning the frequency spectrum of f = 0.001–300 Hz [43]. Fundamentally, the Earth is conceptualized as a horizontal medium, with the magnetotelluric fields representing plane electromagnetic waves projected vertically onto the Earth's surface[13]. Upon striking the Earth's surface, a substantial portion of these waves undergoes reflection, while only a minor fraction penetrates into the subsurface. This phenomenon is driven by electromagnetic induction, specifically the fluctuating magnetic field, which induces telluric currents to propagate into the subsurface, with the magnitude of these currents being contingent on the electrical conductivity properties [44]. The skin depth (δ), representing the depth within the subsurface where electromagnetic wave amplitude diminishes to 1/e of its value at the surface, is mathematically expressed as follows [45].
δ = ρ π f μ   ≈   500 ρ f
The skin depth is determined by the resistivity (ρ) measured in ohm meters (Ω.m), the frequency (f) expressed in hertz (Hz), and the magnetic permeability (μ) in henry per meter (H/m). It is notable that the primary factors governing the skin depth are the conductivity (reciprocal of resistivity) of geological formations and the operational frequency employed. Geological formations exhibiting enhanced conductivity in the subsurface are commonly associated with materials such as graphite or carbon films, interconnected metallic minerals, aqueous fluids and partial melt [46].
The above provided expression illustrates that the skin depth exhibits an inverse relationship with frequency. Consequently, lower frequencies possess the ability to penetrate to greater depths, whereas higher frequencies are confined to shallower regions. At the Earth's surface, orthogonal electromagnetic field components are observed, and these components provide insights into the frequency response, reflecting the distribution of electrical properties within the subsurface medium [47]. The variation in the magnetotelluric field component over time is transformed into a frequency spectrum, enabling the computation of magnetotelluric frequency domain responses, such as apparent resistivity and impedance phase. The calculation of apparent resistivity can be expressed as follows:
ρ = 1 5 f E x H y 2
In this equation, f represents the frequency in hertz (Hz), ρ signifies the resistivity in ohm meters (Ω.m), Ex denotes the electric field x-component, and Hy represents the magnetic field y-component.

4.2 Dual Frequency induced Polarization Method

The Dual-Frequency Induced Polarisation (DFIP) method, innovated by Chinese academician Jishan He [49], represents a significant advancement in geophysical exploration techniques, particularly in the domain of mineral exploration [48]. The (DFIP) method is an advanced geophysical technique operating within the induced polarisation (IP) frequency domain. This system integrates both a transmitter and a receiver, essential for its functionality. The core mechanism relies on the utilization of both high frequency and low-frequency electrical currents. The transmitter, a pivotal component of the DFIP system, is tasked with synthesizing and energizing the electromagnetic field source. This generated field is then strategically deployed into the subsurface geological layers. When this field interacts with subterranean rock formations and ore deposits, it induced polarization effects within these materials [49]. The DFIP system's receivers are specifically designed to detect these induced polarisation responses. The distinctive feature of this methodology is its proficiency in detecting the differential frequency characteristics presented by diverse rock and ore types. The principal measurements obtained by this system include the high-frequency potential difference (∆VH), the low-frequency potential difference (∆VL), from which we calculate the resistivity (Ohm-m) and the percent frequency effect (Fs) (Table 1). These parameters are critical in analyzing the subsurface geological structures and identifying potential mineral deposits [49].

4.3 Data Acquisition

In this research, we carried out an AMT survey at the Dongjun deposit area, deploying 31 AMT stations along a 1440-meter profile line aligned perpendicular to the geological strike (figure 3). The stations were spaced 40 meters apart. We utilized the GSEM-W10 system, developed by Giant Sequoia Artificial Intelligence Technology Co., Ltd., for collecting time-variant field or time series data. The AMT data encompassed a frequency range from 1 Hz to 10,400 Hz, covering 53 frequencies. At each station, data collection lasted for 35 minutes, capturing two horizontal electric field components (Ex and Ey) and two orthogonal magnetic field components (Hx and Hy). The orientations for the X- and Y-directions were north and east, respectively. We measured magnetic field variations using induction coil magnetometers (ICM) and electric field variations with two pairs of non-polarizable lead-chloride electrodes (Pb-PbCl2). To lower contact resistivity, each electric field measurement point was pre-saturated with water. We also assessed each station to identify potential sources of interference, such as roads, high-voltage power lines, and communication cables, which were prevalent along the survey lines. The AMT time series field data collected were processed using the GSEM-pros software. This involved converting the data to the frequency domain and calculating the cross-power spectra. The cross-power spectra calculations were crucial for estimating the impedance tensor, which varies with frequency. This impedance tensor is key for understanding the dimensionality and strike direction of subsurface structures in our study area.
In the same survey area, the Dual-Frequency Induced Polarisation (DFIP) method using a pole-dipole array configuration was conducted (figure 3). The receiver array, comprising non-polarizable potential electrodes (MN), was arranged at 40-meter intervals along a 720-meter length. Data acquisition were performed at 31 stations using dual frequencies of 4 Hz and a combination of 4/13 Hz, utilizing the SQ-3C model. The measurement process initiated with an initial current electrode (AB) spacing of 80 meters. The supply dipole's center was positioned at the midpoint of the survey line. Subsequently, the distance of the power-supply dipole within the survey line was progressively increased by 80 meters for each measurement, following a linear sequence from 80 meters up to 800 meters (in increments of 80, 160, 240, up to 800 meters). Beyond the survey line, the distance of the power-supply dipole (AB) was extended linearly in steps, starting from 960 meters and continuing until it reached the maximum distance of 1440 meters (specifically at intervals of 960, 1200, and 1440 meters).

5. Inversion

5.1. AMT Data

In this research, we conducted a two-dimensional inversion of AMT data using the Occam inversion algorithm, using ZONDMT2D software package [50]. The Static shift, an artifact stemming from near surface in homogeneities, has the potential to introduce significant complexities and misleading interpretation of apparent resistivity curves. Such complexities could lead to inaccuracies in delineating geoelectric structures. Therefore, we corrected for static shift by using manual adjustments of curve levels while referencing adjacent curves. These adjustments were executed using the ZONDMT2D software.
For AMT inversion, several critical parameters and procedures were selected to ensure the precision of our results. We initiated the inversion process by selecting height value of 5 and applying an incremental factor of 1.05, which systematically adjusted model parameters in subsequent iterations. An initial half-space resistivity value of 250 Ωm, smoothing factor of 1, depth smoothing of 1, a smoothness ratio of 0.5 was selected to optimize the inversion process based on conditions specified in reference material [50]. Additionally, we imposed common model limits, setting a minimum resistivity of 10 Ωm and a maximum resistivity of 10,000 Ωm to bound the range of resistivity values within the model. The inversion process was executed iteratively, initially spanning 20 iterations and subsequently undergoing 10 more iterations with parameter adjustments aimed at improving the root mean square (RMS) error. The RMS error, quantifying the difference between the modeled and measured datasets, was assessed as a percentage, with a stopping criterion of RMS errors below 10% considered acceptable for our two-dimensional inverse modeling.

5.2 DFIP Data

We initiated the inversion process by selecting height value of 5 and incremental factor of 1.05, which systematically adjusted model parameters in subsequent iterations. An initial half-space resistivity value of 300 Ωm was established as a reference point for subsurface resistivity. The parameters including a smoothing factor of 0.01, depth smoothing of 1, a smoothness ratio of 1, and a focused threshold of 0.05, each selected to optimize the inversion process based on conditions specified in reference [51].
The inversion results showed a root mean square (RMS) error for the resistivity model and PFE inversion model below 10%, reflecting a high level of accuracy in our results.

5.3 Joint inversion Data

We used a joint inversion method to combine data sets from DFIP and AMT surveys, using ZONDRES2D software[51]. We selected settings similar to those we used for DFIP. During the joint inversion, an important step was merging both sets of data and adjusting how much weight the MT data carried in the software. We fine-tuned this weight within a range from 0.25 to 1, aiming for a closer match between the phases of the data to reduce the misfit. Our joint inversion analysis produced a root mean square (RMS) error below 10% for the inversion Model.

5.4. Inversion results

5.4.1. AMT result

The resistivity model derived from AMT data over the Dongjun lead-zinc deposit reveals distinct features characterized by low electrical resistivity. The inversion model (Figure 4) identifies three zones with moderate to high conductivity (resistivity < 600 Ω.m), labeled as C1, and C2, and two zones with high resistivity (resistivity > 3500 Ω.m), designated as R1 and R2. Notably, the C1 zone is predominantly situated in the surface layers, exhibiting a horizontal extension towards the NE and could be attributed to the quaternary alluvial deposits of clays and silt. In contrast, the R1 and R2 zones are located in the SW and NE sections of the inversion model and be attributed to the basaltic rocks. These resistive zones are notably interspersed by the conductive zone C2. In the middle portion of the inversion model, the conductive bodies C2 is prominent and vertically continuous and could be related to the lead-zinc ore.

5.4.2. DFIP result

The resistivity model derived from DFIP survey conducted over the same area shows numerous unique characteristics. The resistivity model (figure 5a) showed one high conductive body (resistivity < 150 Ω.m) labeled as C1 and two moderates to high resistive bodies (resistivity < 3000 Ω.m) as R1 and R2. The C1 conductive layer is present in the surficial layer, extends in NE direction, and could be attributed to the quaternary alluvial deposits of clays and silt. The moderate resistive body R1 is present in the SW section of the inversion model and could be related to sandy soil gravel or conglomerates, while the high resistive body R2 is present in the middle of the inversion model that could be attributed to basaltic rocks like Rhyolite and Breccia. In contrast with the AMT inversion results, the DFIP shows the presence of a high resistivity body R2 in the central part of the survey, where the AMT result identified it as a low-resistivity zone. However, the Percent frequency effect (PFE) of the same anomaly shows values greater than 12 percent indicating IP anomaly as shown in Figure 5b, which could be attributed to disseminated lead-zinc ore within a basaltic rock host.

5.4.3. Joint inversion result

The resistivity model derived from Joint inversion, as depicted in (Figure 6a) showed a range of subsurface anomalies that differ from the individual models. These range from highly resistive to highly conductive structures, encompassing a highly conductive layer (C1), a moderately resistive body (R1), a structure that has high resistivity (R2) at the center of the profile, which overlies a conductive body (C2). The (C1), high conductive is present in the surficial portion of the inversion model have resistivity value less 150 Ωm, which could be attributed to alluvial deposits of clays and silt, while (R1) is located in SW section and its resistivity range is less than 2000 Ωm may be related to sandy gravel or conglomerates. In the upper-middle segment of the inversion model, a prominently high resistivity structure, designated as (R2), exhibiting a resistivity exceeding 6000 Ωm at the center of the profile is attributed to volcanic rocks such as Rhyolite and Breccia. Situated beneath this high resistivity zone is a moderately to highly conductive feature, labeled as (C2), characterized by resistivity values below 700 Ωm, and displaying a notable vertical continuation could be attributed to disseminated lead-zinc deposits within silicate host. Additionally, the PFE value associated with this region is greater than 14 percent (figure 6b), which confirms the presence of IP anomaly related lead-zinc mineralization.

6. Discussion

The AMT inversion model (Figure 4) showed several areas of high resistivity and high conductivity along the Profile from SW to NE sections. In the SW section of the profile, the strata from the shallow to greater depths are observed to be highly resistivity range between 800-4000 Ωm. The high resistive structure (R1) extends along 400 m (from -500 to -100) along the horizontal distance and located at a depth of 400 m (from 700 to 300) from the surface. After studying geology of the area and resistivity signature, we could interpret that the (R1) structure could be attributed to the basaltic rocks like rhyolite, breccia tuff, andesitic basalt, volcanic basalt and breccia, belonging to Manketouebo and Tamulangou formations.
The central section of the profile has three distinct anomalous features, which area several units of high conductivity structures (C1) in shallow depths and moderate to high conductive body (C2) in the middle of the profile and extending downward. The high conductivity structure (C1) has resistivity <150 Ωm and extends 500 m horizontally (from 0 to 500) and 100 m vertically (from 700 to 600) with a northeast orientation. Based on their resistivity signature and the geology of the area we interpret the (C1) layer to be quaternary deposits of silt, and alluvial clay and that the enhanced conductivity is due to fault (F1). The significantly low electrical resistivity of the conductor C2, (<600 Ωm) is restricted by 400 m of horizontal distance (from -100 to 300) and is evident from the minimum depth of 300 m from the surface (from 600to 300), which extends into deeper strata. We deduce that this deep-seated low resistivity anomaly is due to mineralization of lead and zinc base metals. This is supported by the PFE model and geological cross section which coincide and identify the location of the highest conductors in the study area. The central section consists of fault: F1, (Figure 4) which are likely to exhibit low resistivity. According to Unsworth (1999), faults can increase permeability parallel to the faults line and inhibit the movement of materials perpendicular to the fault plane. In addition, faults can trap fluids in gouge- and fault- breccia. In this work, conductivity enhancing effects of fault F1 are evident in the AMT resistivity model.
The NE section of the profile consists of high resistivity structure (R2) with a resistivity range from 800-4000 Ωm and strikes southeast and is 150 m (point 350 to -500) in length and 300m (point 600 to 300) in height. We interpret that the highly resistive structures (R1 and R2) are composed of rhyolite, breccia tuff, andesitic basalt, volcanic basalt and breccia of Manketouebo and Tamulangou formations based on geology and resistivity signature. The resistivity signature of the structures presents in the of SW and NE section are similar and little enhancement in resistivity is attributed to silicification effect.
In the DFIP inversion model showed in (Figure 5a), a resistive body (R1) and conductive body (C1) are observed at a shallow depth. The moderate to high resistive body (R1) have resistivity >3500 Ωm, constrained within a horizontal distance of 50 m (from -310 to -270) and located at a depth of 130 m from the surface. Based on the resistivity signature and geological context, this resistive body (R1) is interpreted as basaltic rocks such as andesitic basalt, volcanic basalt, and volcanic breccia from the Tamulangou (J2tm) formation. The high conductive body (C1) with resistivity < 150 Ωm, extending over a horizontal distance of 650 m (from -250 to 400) in the SW-NE direction (from 700 to 400). This conducting body shows the presence of Quaternary alluvial clay and silt deposits, based on the geology and resistivity signature. However, the main contradiction between the resistivity models of DFIP and AMT is in the middle of the profile. In the middle of the DFIP inversion model, a high resistive body (R2) with resistivity > 7000 Ωm is observed, along a horizontal distance of 450 m (from -150 to 300) and extending with depth (from 650m to 350m) which is different from the AMT inversion model results as well as the electrical measurement results of rock and ores specimens from the study area. The electrical measurement results of rock cores and outcrops conducted in the mining area (Table 1) and the ore bearing specimens generally exhibit high Percent frequency effect value and low resistivity anomalies. In reference to [31] the Dongjun lead-zinc deposit is composed of rhyolite, basalt, andesitic basalt, volcanic tuff, and exhibits alteration characterized by high levels of silicification. Considering the geology, literature review, and the high resistivity signature, it can be inferred that the high resistive structure (R2) observed in the DFIP inversion model, unlike the low resistive structures (C2) observed in AMT, is attributed to disseminated mineralization within a silica-rich host, which creates a resistive barrier that increase the resistivity in DFIP but increase the PFE values. The significant IP anomaly zone is located in the middle of the profile (Figure 5b), with a maximum PFE value of 12 percent. This suggests that the IP anomaly zones are associated with deep lead-zinc mineralization.
The resistivity model derived from the joint inversion of AMT and DFIP data sets improve the subsurface structure and resolve structures that are not shown by DFIP. The joint inversion resistivity model in (Figure 6a) shows surficial conducting layer (C1), highly resistive body (R1) and moderate to high conductive structure (C2) along the profile. A conductive zone (C1) which have (<150 Ωm) is evident in the surficial depth of 50m from surface with a horizontal extension of 600m (from -100-500) and extend in the NE direction. We interpret that the conductive layer (C1) to be due to Quaternary alluvial clay and mudstone deposit based on their resistivity signature and in reference to the geology of the area. The (R1) resistive body is embedded at the upper layers and extend at about 60m depth from surface and horizontal extension of 100m (from no -400 to -300), that could be interpreted as breccia, volcanic basalt or andesitic basalt of Tamulangou (J2tm) formation. The high resistive zone (R2) which have (>6000Ωm) are located just below the surficial conductive body (C1). The high resistive body (R2) is located along 450m (from -280 to 220) along the horizontal distance and is located at a depth of 250m from surface which extends with depth. In reference to the geology of the area, resistivity signature and literature review we interpret that the high resistive zone (R2) reflect basaltic rocks such as, rhyolite, andesitic basalt, breccia tuff of the of Manketouebo and Tamulangou formations. Below the resistive alteration zone there is a large area of moderate to high conductivity (C2) are recovered in joint inversion resistivity model which have resistivity range (<700 Ωm) and distributed along 750m (from -450 to 300) along the horizontal distance and is located at a depth of about 450m from the surface and extends with depth. The obtained resistivity model coincides with the known location of the lead–zinc mineralization in the study area (Figure 7a). In addition to high conductivity signature there is a high Percent frequency effect (PFE) anomalous zones showed in (Figure 6b) which confirms the location of IP anomalous zone. Based on all these inversion models along with the geological model (Figure 7b) and the ore specimen’s electrical resistivity (table 1), we interpret that the high conductive body (C2) is due to deep disseminated metal lead-zinc. This anomaly is clearer and more strongly supported by the joint inversion technique than by AMT or DFIP s individually.

7. Conclusions

The mineralization in the study area is known to occur within areas of silicification, which is a resistive geophysical target. However, two individual geophysical surveys AMT and DFIP over the Dongjun deposit image the resistivity signature differently. To solve the problem and comprehensive understanding of the subsurface structures we undertook a joint inversion approach. The joint inversion resistivity model shows surficial conducting layer (C1), moderately resistive body (R1) highly resistive body (R2) and resolve area of moderate to high conductive body (C2) extending along the profile and with depth. The final joint inversion resistivity model clearly images the large silica alteration zone and area of disseminated Pb-Zn sulfide mineralization based on low resistivity signature and high Percent frequency effect PFE anomaly. Further analysis from this final model highlights potential areas of sulfide and other mineralization within the silica alteration zone in the form of small or large conductive anomalies. Collectively, this study illustrates that integration of geophysical methods through joint inversion is a valuable approach for improved characterization and refinement of exploration-targeting Pb-Zn mineralization zones in future studies.

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org. Figure S1: title; Table S1: title; Video S1: title.

Author Contributions

L.C., C.R., designed the project; I.A., J.A., and O.A. conducted the original literature reviews; F.U., O.R., and S.S, write the methodology. S.F., L, C, C.R, writing the original paper with a careful discussion and revision by C.R. H.E.K., critically revised the manuscript for important intellectual and technical content and contributed in the interpretation of the inversion results. All authors have read and agreed to the published version of the manuscript.

Funding

This research is funded by Basic Science Center Project of National Natural Science Foundation of China, grant number 72088101.

Data Availability Statement

The Data files for AMT, DFIP and coordinate information of all sites can be obtained on request by contacting the corresponding author.

Acknowledgments

We would like to thank the anonymous reviewers and academic editors for taking the time and effort to review the manuscript. We appreciate all valuable comments and suggestions, which improved the quality of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gonzalez-Alvarez, I.; Goncalves, M.; Carranza, E.J.M. Introduction to the special issue challenges for mineral exploration in the 21st century: Targeting mineral deposits under cover. Ore Geology Reviews 2020, 126, 103785. [Google Scholar] [CrossRef]
  2. Wood, D.; Hedenquist, J. Mineral exploration: Discovering and defining ore deposits. SEG Newsletter 2019, 1–22. [Google Scholar]
  3. Arnd, L.; Fo, Í.; Jeffrey, W.; Kesler, S.E.; Iel, G. D future global mineral resources. <bold>2017</bold>.
  4. Okada, K. Breakthrough technologies for mineral exploration. Mineral Economics 2022, 35, 429–54. [Google Scholar] [CrossRef]
  5. Jiang, W.; Duan, J.; Doublier, M.; Clark, A.; Schofield, A.; Brodie, R.C.; Goodwin, J. Application of multiscale magnetotelluric data to mineral exploration: An example from the east tennant region, northern australia. " Geophysical Journal International 2022, 229, 1628–45. [Google Scholar] [CrossRef]
  6. Gan, J.; Li, H.; He, Z.; Gan, Y.; Mu, J.; Liu, H.; Wang, L. Application and significance of geological, geochemical, and geophysical methods in the nanpo gold field in laos. Minerals 2022, 12, 96. [Google Scholar] [CrossRef]
  7. McMillan, M.S.; Oldenburg, D.W. Cooperative constrained inversion of multiple electromagnetic data sets. Geophysics 2014, 79, B173–B85. [Google Scholar] [CrossRef]
  8. Bagchi, A.; Singh, R.K.; Sinharay, R.K. 3d phase tensor inversion and joint inversion using impedance and phase tensor components aiming better subsurface mapping and interpretation in dhanjori basin, singhbhum craton, eastern india. Journal of Applied Geophysics 2022, 202, 104646. [Google Scholar] [CrossRef]
  9. Chen, R.; He, M.; Chang, F.; Zeng, P.; Zhao, X. Large-scale 3-d inversion of amt data with application to mineral exploration. Paper Presented at the 23rd Electromagnetic Induction Workshop—EMIW2016, Chiang Mai, Thailand, 2016; pp. 14–21. [Google Scholar]
  10. He, G.L.; Wang, G.J.; Zhou, C.; Yan, X.L.; Li, H.L. Application of audio-frequency magnetotelluric (amt) in groundwater exploration: A case of the nanmushui area in zhanjiang. Progress in Geophysics 2019, 34, 304–09. [Google Scholar]
  11. Deshmukh, V.; Kumar, P.V.; Rao, P.S.; Kumar, A.; Singh, A. Audiomagnetotelluric (amt) studies across aravali-tural-rajawadi geothermal zones, western maharastra, india. Journal of Applied Geophysics 2022, 198, 104579. [Google Scholar] [CrossRef]
  12. Ishizu, K.; Ogawa, Y.; Mogi, T.; Yamaya, Y.; Uchida, T. Ability of the magnetotelluric method to image a deep conductor: Exploration of a supercritical geothermal system. Geothermics 2021, 96, 102205. [Google Scholar] [CrossRef]
  13. He, L.; Chen, L.; Dorji, *!!! REPLACE !!!*; Xi, X.; Zhao, X.; Chen, R.; Yao, H. Mapping the geothermal system using amt and mt in the mapamyum (qp) field, lake manasarovar, southwestern tibet. Energies 2016, 9, 855. [Google Scholar] [CrossRef]
  14. Wu, G.; Hu, X.; Huo, G.; Zhou, X. Geophysical exploration for geothermal resources: An application of mt and csamt in jiangxia, wuhan, china. Journal of Earth Science 2012, 23, 757–67. [Google Scholar] [CrossRef]
  15. Revil, A.; Vaudelet, P.; Su, Z.; Chen, R. Induced polarization as a tool to assess mineral deposits: A review. Minerals 2022, 12, 571. [Google Scholar] [CrossRef]
  16. Li, J. Application of the down-hole ip method to a general survey in the jinya gold deposit. Geology & Exploration, 2016; 924–30. [Google Scholar]
  17. Kearey, P.; Brooks, M.; Hill, I. An introduction to geophysical exploration; John Wiley & Sons, 2002. [Google Scholar]
  18. Gallardo, L.A.; Meju, M.A. Joint two-dimensional dc resistivity and seismic travel time inversion with cross-gradients constraints. Journal of Geophysical Research: Solid Earth 2004, 109. [Google Scholar] [CrossRef]
  19. Albouy, Y.; Andrieux, P.; Rakotondrasoa, G.; Ritz, M.; Descloitres, M.; Join, J.L.; Rasolomanana, E. "Mapping coastal aquifers by joint inversion of dc and tem soundings-three case histories. Groundwater 2001, 39, 87–97. [Google Scholar] [CrossRef]
  20. Vozoff, K.; Jupp, D. Joint inversion of geophysical data. Geophysical Journal International 1975, 42, 977–91. [Google Scholar] [CrossRef]
  21. Haber, E.; Oldenburg, D. Joint inversion: A structural approach. Inverse problems 1997, 13, 63. [Google Scholar] [CrossRef]
  22. Chen, B.; Jahn, B.-M.; Wilde, S.; Xu, B. Two contrasting paleozoic magmatic belts in northern inner mongolia, china: Petrogenesis and tectonic implications. Tectonophysics 2000, 328, 157–82. [Google Scholar] [CrossRef]
  23. Xiao, W.; Windley, B.F.; Allen, M.B.; Han, C. Paleozoic multiple accretionary and collisional tectonics of the chinese tianshan orogenic collage. Gondwana Research 2013, 23, 1316–41. [Google Scholar] [CrossRef]
  24. Chen, Z.; Zhang, L.; Wan, B.; Wu, H.; Cleven, N. Geochronology and geochemistry of the wunugetushan porphyry cu–mo deposit in ne china, and their geological significance. Ore Geology Reviews 2011, 43, 92–105. [Google Scholar] [CrossRef]
  25. Xiao, W.; Windley, B.F.; Hao, J.; Zhai, M. Accretion leading to collision and the permian solonker suture, inner mongolia, china: Termination of the central asian orogenic belt. Tectonics 2003, 22. [Google Scholar] [CrossRef]
  26. Wu, F.-Y.; Sun, D.-Y.; Ge, W.-C.; Zhang, Y.-B.; Grant, M.L.; Wilde, S.A.; Jahn, B.-M. Geochronology of the phanerozoic granitoids in northeastern china. Journal of Asian Earth Sciences 2011, 41, 1–30. [Google Scholar] [CrossRef]
  27. Tang, J.; Xu, W.-L.; Wang, F.; Zhao, S.; Wang, W. Early mesozoic southward subduction history of the mongol–okhotsk oceanic plate: Evidence from geochronology and geochemistry of early mesozoic intrusive rocks in the erguna massif, ne china. Gondwana Research 2016, 31, 218–40. [Google Scholar] [CrossRef]
  28. Ouyang, H.-g.; Mao, J.-W.; Santosh, M.; Zhou, J.; Zhou, Z.-H.; Wu, Y.; Hou, L. Geodynamic setting of mesozoic magmatism in ne china and surrounding regions: Perspectives from spatio-temporal distribution patterns of ore deposits. Journal of Asian Earth Sciences 2013, 78, 222–36. [Google Scholar] [CrossRef]
  29. Porter, T.M. The geology, structure and mineralisation of the oyu tolgoi porphyry copper-gold-molybdenum deposits, mongolia: A review. Geoscience Frontiers 2016, 7, 375–407. [Google Scholar] [CrossRef]
  30. Ge, W.; Wu, F.; Zhou, C.; Zhang, J. Porphyry cu-mo deposits in the eastern xing’an-mongolian orogenic belt: Mineralization ages and their geodynamic implications. Chinese Science Bulletin 2007, 52, 3416–27. [Google Scholar] [CrossRef]
  31. Xie, W.; Wen, S.; Zhang, G.; Tang, T. Geochronology, fluid inclusions and isotopic characteristics of the dongjun pb-zn-ag deposit, inner mongolia, ne china. Acta Geologica Sinica-English Edition 2021, 95, 1611–33. [Google Scholar] [CrossRef]
  32. Liu, J.; Mao, J.; Wu, G.; Wang, F.; Luo, D.; Hu, Y. Zircon u–pb and molybdenite re–os dating of the chalukou porphyry mo deposit in the northern great xing’an range, china and its geological significance. Journal of Asian Earth Sciences 2014, 79, 696–709. [Google Scholar] [CrossRef]
  33. Feng, Z.; Jia, J.; Liu, Y.; Wen, Q.; Li, W.; Liu, B.; Xing, D.; Zhang, L. Geochronology and geochemistry of the carboniferous magmatism in the northern great xing’an range, ne china: Constraints on the timing of amalgamation of xing’an and songnen blocks. Journal of Asian Earth Sciences 2015, 113, 411–26. [Google Scholar] [CrossRef]
  34. Duan, T.; Wen, S.; Zhou, P.; Zhu, E.; Cui, X.; Sun, J.; Zhao, T.; Wei, C. Wall rock alteration characteristics and mass balance calculation of dongjun pb-zn-ag deposit in inner mongolia. Journal of Central South University (Science and Technology) 2018, 49, 1991–2002. [Google Scholar]
  35. Liu, H.; Yang, X.; Yu, C.; Ye, J.; Liu, J.; Zeng, Q.; Shi, K. A case study in finding concealed ores by using geophysical exploration methods in combination of vlf-em, eh4 and csamt. Progress in Geophysics 2004, 19, 276–85. [Google Scholar]
  36. Tikhonov, A. On determining electrical characteristics of the deep layers of the earth's crust. Dokl. Acad. Nauk. Sssr 1950, 73, 295–297. [Google Scholar]
  37. Cagniard, L. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics 1953, 18, 605–35. [Google Scholar] [CrossRef]
  38. Cantwell, T. Detection and analysis of low frequency magnetotelluric signals; Massachusetts Institute of Technology, 1960. [Google Scholar]
  39. Cantwell, T.; Madden, T. Preliminary report on crustal magnetotelluric measurements. Journal of Geophysical Research 1960, 65, 4202–05. [Google Scholar] [CrossRef]
  40. Vozoff, K. The magnetotelluric method in the exploration of sedimentary basins. Geophysics 1972, 37, 98–141. [Google Scholar] [CrossRef]
  41. Singh, R.K.; Maurya, V.P.; Singh, S. Imaging regional geology and au–sulphide mineralization over dhanjori greenstone belt: Implications from 3-d inversion of audio magnetotelluric data and petrophysical characterization. Ore Geology Reviews 2019, 106, 369–86. [Google Scholar] [CrossRef]
  42. Tezkan, B. Radiomagnetotellurics. In Groundwater geophysics: A tool for hydrogeology; Springer, 2009; pp. 295–317. [Google Scholar]
  43. Kalscheuer, T.; Juhojuntti, N.; Vaittinen, K. Two-dimensional magnetotelluric modelling of ore deposits: Improvements in model constraints by inclusion of borehole measurements. Surveys in Geophysics 2018, 39, 467–507. [Google Scholar] [CrossRef]
  44. Chave, A.D.; Jones, A.G. The magnetotelluric method: Theory and practice; Cambridge University Press, 2012. [Google Scholar]
  45. Nabighian, M.N. Electromagnetic methods in applied geophysics: Voume 1, theory; Society of Exploration Geophysicists, 1988. [Google Scholar]
  46. Lindsay, M.; Spratt, J.; Occhipinti, S.; Aitken, A.; Dentith, M.; Hollis, J.; Tyler, I. Identifying mineral prospectivity using 3d magnetotelluric, potential field and geological data in the east kimberley, australia. Geological Society, London, Special Publications 2018, 453, 247–68. [Google Scholar] [CrossRef]
  47. Li, F.; Zeng, Q.; Zhu, R.; Chu, S.; Xie, W.; Zhang, B.; Zhang, X. Application of the amt method to gold deposits: A case study in the qinling metallogenic belt of north china craton. Minerals 2021, 11, 1200. [Google Scholar] [CrossRef]
  48. Liu, W. Application of dual frequency ip method in southern peru copper deposit exploration. Advanced Materials Research 2012, 529, 311–17. [Google Scholar] [CrossRef]
  49. He, J. Dual frequency induced polarization method. Higher education press 2006. [Google Scholar]
  50. ZOND. "Zondmt2d." 2001-2023. Available online: http://zond-geo.com/english/zond-software/electromagnetic-sounding/zondmt2d/ (accessed on 4 December 2023).
  51. ZOND. "Zondres2d." 2001-2023. Available online: http://zond-geo.com/english/zond-software/ert-and-ves/zondres2d/ (accessed on 4 December 2023).
Preprints 108150 g001
Figure 1. (a) Simplified tectonic map showing the main units of central and eastern Asia and (b) sketch tectonic map of the Songnen block (modified from [30,33]).
Figure 1. (a) Simplified tectonic map showing the main units of central and eastern Asia and (b) sketch tectonic map of the Songnen block (modified from [30,33]).
Preprints 108150 g001
Preprints 108150 g002
Figure 2. Stratigraphic sequence of the Dongjun deposit.
Figure 2. Stratigraphic sequence of the Dongjun deposit.
Preprints 108150 g002
Preprints 108150 g004
Figure 4. 2D inversion model of AMT data with interpreted geoelectrical model and possible fault.
Figure 4. 2D inversion model of AMT data with interpreted geoelectrical model and possible fault.
Preprints 108150 g004
Preprints 108150 g005
Figure 5. (a) 2D inversion model of DFIP with interpreted subsurface geoelectrical structures and showing high resistivity zone (b) inversion model of Percent frequency effect (PFE) showing IP anomaly related to mineralization zone.
Figure 5. (a) 2D inversion model of DFIP with interpreted subsurface geoelectrical structures and showing high resistivity zone (b) inversion model of Percent frequency effect (PFE) showing IP anomaly related to mineralization zone.
Preprints 108150 g005
Preprints 108150 g006
Figure 6. (a) joint inversion resistivity model of AMT and DFIP showing area of high resistivity zone above conductive zone. (b) Inversion model of Percent frequency effect (PFE) showing IP anomaly related to mineralization zone.
Figure 6. (a) joint inversion resistivity model of AMT and DFIP showing area of high resistivity zone above conductive zone. (b) Inversion model of Percent frequency effect (PFE) showing IP anomaly related to mineralization zone.
Preprints 108150 g006
Preprints 108150 g007
Figure 7. (a) Resistivity model based on joint inversion resistivity of AMT and DFIP overlapped with the geological model shwoing location and trend of Pb-Zn ore body. (b) Percent frequency effect inversion model overlap with geological cross section showing trend and location of the ore body.
Figure 7. (a) Resistivity model based on joint inversion resistivity of AMT and DFIP overlapped with the geological model shwoing location and trend of Pb-Zn ore body. (b) Percent frequency effect inversion model overlap with geological cross section showing trend and location of the ore body.
Preprints 108150 g007
Table 1. Rock and ore specimen’s electrical resistivity in the study area.
Table 1. Rock and ore specimen’s electrical resistivity in the study area.
Rock ore specimen Resistivity range ρ (Ωm) PFE (%)
Basalt 1000 ~ 2000 0.7 ~ 1.5 %
Rhyolite 1500 ~ 2800 0.5 ~ 1.2 %
Breccia 2000 ~ 3000 0.3 ~ 1.0 %,
Pb-Zn Ore 500 ~ 1000 ~ 15.0%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated