3D Multicomponent Self-Potential Inversion: Theory and Application to the Exploration of Seaﬂoor Massive Sulﬁde Deposits on Mid-Ocean Ridges

: The marine self-potential (SP) method is currently playing an increasing role in the exploration and resource evaluation of seaﬂoor massive sulﬁde (SMS) deposits. SP surveys are conducted using autonomous underwater vehicles (AUV), which yield multicomponent electric ﬁeld datasets. By comparing with the single-component electrical ﬁeld data used to date, the inversion of these multicomponent data is expected to provide a more accurate description of the 3D structure of SMS deposits beneath the seaﬂoor (like gradiometry in gravity surveys). We introduce an inversion al-gorithm speciﬁcally adapted to multicomponent SP data. A synthetic model demonstrates that the inversion of multicomponent datasets allows us to be�er recover the amplitude of the current density and the morphology of the ore bodies compared to using a single component of the electrical ﬁeld. Next, we apply our approach to a multicomponent SP dataset collected during the DY58 oceanic cruise at the Yuhuang hydrothermal ﬁeld on the Southwest Indian Ridge. Subsequently, we reconstruct the three-dimensional (3D) geometry of the SMS deposits beneath the seaﬂoor. The AUV-based SP system with the collection of multicomponent SP data inversion appears to be a powerful tool in the exploration and evaluation of seaﬂoor sulﬁde resource and, in the future, could be used in concert with induced polarization data.


Introduction
Seafloor massive sulfide (SMS) deposits are rich in metals such as copper, zinc, silver, and gold, which are regarded as important seafloor mineral resources for humankind [1]. These deposits can be explored and evaluated with geophysical methods. Electromagnetic methods have been broadly applied to SMS deposits. They include the magnetic method, the transient electromagnetic method, the controlled-source electromagnetic method, the direct current resistivity method, and the self-potential method [2][3][4][5][6][7][8]. The self-potential method is a passive method requiring simple equipment, including a voltmeter and a set of non-polarizing electrodes [9]. It can be easily mounted on an AUV. In recent years, the SP method has been the focus of intensive research in the context of the SMS-deposit exploration [10][11][12][13]. Meanwhile, SP can be conducted with active electromagnetic methods such as the transient electromagnetic method and/or galvanometric and induced polarization techniques [14][15][16].
SMS deposits associated with active and inactive hydrothermal systems are the setting of distinct redox reactions [17] responsible for negative self-potential anomalies just above the seafloor. Inactive hydrothermal fields characterized by the absence of water column geochemical anomalies are more mature (older) than active hydrothermal fields [18]. A silica cap is usually formed above these deposits which is, in turn, beneficial regarding the additional accumulation of sulfides leading at the end to the formation of massive sulfide deposits [19].
The use of non-polarizing (e.g., Ag/AgCl) electrodes to map the distribution of the natural electric field near the seafloor is called the marine SP method [9,20]. On shore, a reference electrode is usually located far away from the ore body, where the electric potential vanished asymptotically to zero. The potential difference between this reference electrode and the scanning electrode is generally negative over sulfide deposits [21]. For marine self-potential surveys, it is far easier to measure the potential difference between a pair of electrodes, i.e., a component of the electrical field obtained according to     E , where φ denotes the electrical potential (the potential difference divided by the corresponding distance between the two electrodes). Therefore, the collected data are a collection of electric field components' time series (the analogy can be drawn to gradiometry in airborne gravity surveys [22]).
The towed marine SP surveys were mostly carried out in the early exploration phase [23][24][25]. The SP system was mounted on a deep-tow camera or transient electromagnetic systems. These types of measurements could only obtain the horizontal or the vertical component of the electrical field. These data were initially used to identify the location of the ore body and not to determine its 3D morphology [23][24][25]. The obvious disadvantage of the towed measurements is that the relative position of a pair of electrodes cannot be guaranteed to be constant (see discussion in [26]). At the same time, the cable moving in the Earth's magnetic field may generate significant levels of noise [10], therefore reducing the signal-to-noise ratio of the recorded SP data. Recently, we took advantage of the rapid development of investigation equipment, for instance, AUV, Remotely Operate Vehicle (ROV), and Human-Occupied Vehicle (HOV). These new pieces of equipment have been used for SMS-deposit surveys [27][28][29]. Several Ag/AgCl electrodes could be mounted on AUVs in different directions, allowing for large-scale self-potential surveys [10,13], and multicomponent self-potential data could be obtained.
In order to obtain the spatial distribution of an ore body below the seafloor, the collected SP data need to be inverted. Patella (1997) [30] proposed a cross-correlation imaging method to obtain the locations of the causative sources of self-potential anomalies below the ground surface. However, the proposed method was only applicable to the case of monopolar charge aggregation. That being said, self-potential anomalies caused by redox reactions such as SMS deposits are usually dipolar in nature [31]. Revil et al. (2001) [32] extended the probability density tomography to the dipole source case, which was successfully applied to ore deposits. For instance, Kawada et al. (2018) [11] used this method at the Izena hydrothermal field of the mid-Okinawa Trough. Biswas and Sharma (2017) [33] proposed a simulated annealing algorithm for the interpretation of natural electric field data from a complex, thick slab-like model. Based on the finite element method for forward simulation, Jardani et al. (2008) [34] obtained the distribution of source current density and thus determined the flow state of the subsurface fluid by performing 3D inversion of natural potential data. They presented various deterministic and stochastic approaches useable for the SP inverse problem, sometimes based on the analogy with electroencephalographic (EEG) signals. Mendonca (2007) [35] used Green's function to simplify the geoelectric model while considering the effects of the surrounding rock conductivity, redox potential gradient, and the conductivity of the ore body on the SP signal and finally obtained the source current density distribution by charge conservation constrained inversion. Their approach was successfully applied in the detection of a gold deposit in Peru (the Yanacocha case study). Zhu et al. (2016) [36] used the least-squares method to invert multiple sources. Miller et al. (2018) [37] obtained the distributions of the hydrothermal fluid flow under the Tongariro volcano in New Zealand by 3D SP inversion. Zhu et al. (2020) [7] inverted the horizontal component electric field to obtain the structure of the ore body, and the 3D inversion of data was initially applied to evaluate a sulfide resource. That being said, using only a single component of the electric field is insufficient for recovering the direction and magnitude of the current density associated with the occurrence of SMS deposits below the seafloor. As depicted in Figure 1, using both the horizontal component (Ex) and vertical component (Ez) of the electric field is expected to improve the 3D characterization of SMS deposits. Figure 1. Sketch of the geo-ba ery model associated with a massive ore body at a hydrothermal field. The redox reactions occur at the surface of the ore body, and the distribution of the gradient of the redox potential drives the flow of charge carriers inside the ore body, so electrical field anomalies are established that are associated with this current density (modified from [7,9,21]).
In this paper, the inversion algorithm of multicomponent self-potential data is first described in detail. Our algorithm was developed thanks to the open-source framework SimPEG [38,39]. A synthetic model was performed to demonstrate that the inversion of multicomponent data yields a source model that be er mimics the true shape of SMS deposits compared with the result obtain with a single component, as performed so far. Finally, the multicomponent self-potential data collected at the Yuhuang hydrothermal field (Southwest Indian Ridge) during the China DY58 cruise were inverted. The spatial distribution of the associated SMS deposits of the Yuhuang hydrothermal field was imaged in 3D.

The Causative Source of Negative SP Anomalies above the Seafloor
We first discuss the occurrence of source current density associated with marine SMS deposits. The occurrence of a source current density associated with an ore deposit requires the presence of electron donors and electron acceptors, as well as an electronic conductor facilitating electron transport over long distances [9]. In the case of ore body corrosion, electrons are directly provided by the ore body, moving from the bo om of the ore body to the shallower part. Figure 1 provides a sketch of the redox reaction associated with an ore body. SMS deposits are rich in pyrite (FeS2), chalcopyrite (CuFeS2), and sphalerite (ZnS). When chalcopyrite and pyrite come into contact, chalcopyrite undergoes easier oxidation than pyrite. Consequently, chalcopyrite decomposes first, with the S (−2) in chalcopyrite being oxidized to S (0), thereby releasing Fe 2+ and Cu 2+ (reaction Equation (1) below). Until chalcopyrite is completely oxidized, pyrite begins to be oxidized: it releases SO4 2− and Fe 2+ (reaction Equation (2) below). The ions released by the reaction further react with the oxygen at the diving surface by flow, diffusion, and electromigration. The oxidation reaction on the surface of the deep ore body can be summarized by the following reaction equations: The half-reaction in the shallow part of the seafloor is as follows: Sato and Mooney (1960) [21] explained that the electrons are provided by other oxidation reactions, such as microbial catalysis [40,41]. In fact, the corrosion of the ore body and the oxidation reaction of other substances may exist at the same time [9], and both of them can be available to act as the anode. As a consequence, electric field anomalies are usually observed above the SMS deposits ( Figure 1).

Governing Equations for the Self-Potential Anomalies
The total current density, J , contains the conduction current density and the source current density [42], where  (in S/m) denotes the conductivity below the seafloor, E (V/m) denotes the quasi-static electric field, and s J (A/m 2 ) denotes the source current density below the seafloor. The conservation of charge in the low-frequency limit of the Maxwell equations yields , so Equation (6) can be wri en as where  (V) denotes the potential field, and v q (A/m 3 ) denotes the charge density per unit volume. Moreover, Equation (7) is the govern equation of self-potential forward modeling. We discretize Equation (7) by using the finite volume method and obtain the discrete equation (see [43] for details), where K is the linear operator between the current density and the electrical field (kernel), v is the volume of the cell, u is the calculated potential filed at the cell center, is inner product matrix (see [44] for more details), D is the divergence matrix, and s q is the source term.
As we know, only a limited amount of self-potential data can be collected above the seafloor during the SP survey. Therefore, the number of observed data, d  , is significantly smaller than the computational data,  , defined at the cell centers. An pretation matrix, i P , that interpolates the field data from the computational point to the measured points is introduced. As we mentioned earlier, we used a pair of electrodes to collect potential different data above the seafloor; thus, the potential different  (8) can be wri en as follows [45]: where [ ] surface current density in the different face. Since the redox potential changes with the depth below the seafloor, the SMS deposits are regarded as the source of an electrical current density (primary source). From Equation (9), it is clear that the self-potential data, d  , are related to the conductivity contrast (secondary source, the 3 3 M M  matrix K is related to the conductivity). Therefore, it is desirable to conduct active source electromagnetic surveys to obtain the conductivity distribution below the seafloor and then carry out self-potential inversion. Indeed, combined TEM and SP data have the ability to separate primary and secondary (ghosts) source current densities in SP tomography to be er image SMS deposits. This was demonstrated at the TAG Hydrothermal Mound by [16]. When the conductivity distribution below the seafloor is not known, it can be considered uniform. In this case, the calculation of the Jacobian matrix, J , is only related to the source current density in the model, and the process of generating the electric field by the source current density can be regarded as a linear system, and the sensitivity matrix only needs to be calculated once in the inversion, so we explicitly computed the sensitivity matrix by using the adjoint problem in this paper.

Jacobian Matrix
Combining the definition of the Jacobian matrix, J , with the governing equation, Equation (9), of the SP method, we obtain the Jacobian matrix as follows: Since the matrix K is a large dense matrix, the inverse of matrix 1  K is difficult to compute. Thus, the Jacobian matrix is not usually calculated according to Equation (10). Directly, we can solve the adjoint problem to explicitly calculate the transpose of the sensitivity matrix, The interpretation matrix, T P , can be wri en as where the element "±1" is in the vector N I ; the location of element "1" and "−1" denotes the location of two different electrodes. Therefore, solving Equation (11) is equivalent to solving the adjoin forward equation N times, Thus, all the solutions to Equation (11) consist of

Multicomponent SP Data Inversion
The goal of the SP inversion is to recover the current density distribution below the seafloor. To reach this goal, an objective (cost) function, ( )  m , is defined as follows: where the regularization parameter,  , is used to balance the data misfit term (first term) and model constraint term (second term), and the optimal value of the regularization parameter is determined using the L-curve criterion (see [34,47] for details). represent the inversion model and the reference model, respectively. The self-potential inversion, similar in essence to other potential field techniques, such as the gravity and magnetic methods, looks to retrieve a model vector of causative current sources; the sources located close to the observation points have a lager contribution to the inversion. In order to mitigate this spurious effect, the deep or long-distance anomalies, the Jacobian matrix is normalized in such a way that the weight of the deep sources is increased. (For the construction of the depth weighting matrix, please refer to [48].) The inversion problem is a linear problem; the solution o p t m corresponds to the minimum of the cost function (13) [9].

The Characteristics of Different Electric Field Components
To demonstrate the advantages of the inversion of multicomponent electric field data, we designed a simple 3D current density model (Figure 2). The mesh of the synthetic forward model shown in Figure 2 is discretized by blocks 5 × 5 × 5 m in the source volume, and the outer padding cells have an increasing width out to the boundaries of the source. The total number of meshes is 98 × 98 × 45 = 432,180. Typical SMS deposits have a shape that is similar to that shown in Figure 1 [49]. Therefore, to mimic such a shape, we use an overturned trapezoidal prism. The trapezoidal prism source volume is positioned at the center of the domain, and the vertical source current density value is set to −10 mA/m 2 to represent a realistic value associated with an SMS deposit. (The negative sign indicates that the direction of the current is downward, as discussed in Section 2, above.) The SMS deposit has a thickness of 40 m, extending from a depth of z = 30 m to z = −10 m. At the top, it has a maximum radius of 80 m, which gradually tapers down to a radius of 20 m at the bo om. Above the ore body, there is a cover layer with a thickness of 20 m. The source current density of the seawater and the surrounding rock below the seafloor was set to 0 mA/m 2 . The measurement points were set at a fixed height (30 m) above the seafloor, and the self-potential survey line is characterized by an S-shape. The distribution of the receivers is shown in Figure 2. The distance between the adjacent sampling points is 5 m. In the simulation, the horizontal and vertical component data are collected using two pairs of electrodes, and the distance of the electrode pairs is 5 m. Electrode #1 (front) and Electrode #2 (back) were used for horizontal component observation, and Electrode #3 (bo om) and Electrode #4 (up) were used for vertical observation, with the observation data arranged as a one-dimensional matrix. The observed data in Equation (2)  The direction of the source current is downward in the ore body, and the resulting anomalous electric field is shown in Figure 1. The direction of the resulting electric field is vertical, pointing downward above the ore body, and the direction of the electric field is inclined downward around the ore body. Thus, the horizontal component of the total electric field will be reversed along the survey line, which is 0 mA/m 2 above the ore body (Figure 3a). The direction of the vertical component of the total electric field is always downward, so the value is the largest above the ore body, and away from the ore body,

Single-Component versus Multicomponent Electric Field Results
The single-component electric field (horizontal component data) and multicomponent electric field (horizontal and vertical component data) were used to perform SP inversion. The purpose was to investigate whether the inversion results of multicomponent data can be used to obtain a more accurate spatial distribution of SMS deposits?
In the SP inversion, the regularization parameter,  , was chosen using the classical cooling strategy [50]. The value  is equal to the estimated ratio of the largest eigenvalue of T J J and T m m W W multiply a scalar factor, to weight the relative contributions of the data misfit term and the model constraint term. In this case, the initial value of the scalar factor is set to 10 to ensure that the value of β is large enough at the beginning so that the model term dominates the data misfit term in the objection function (14). Then, we reduce the value of β until a small misfit is found [51]. The termination condition is RMS equal to 1 (here, RMS is defined as Here, we present the inversion results for both the single-component data ( Figure 4a) and multicomponent data (Figure 4b). A comparison between the two reveals that the inversion results for the multicomponent data exhibit a more focused distribution of current density, closely resembling the morphology of the ore body. Histograms of the inverted current density are shown in Figure 4c (single-component data) and Figure 4d (multicomponent data). It is evident from these histograms that the anomalies obtained via the inversion of the multicomponent data are closer to the true model (−10 mA/m 2 ). Furthermore, the predicted data obtained from both inversions show a good fit with the observed electric field values (Figure 4e,f). Consequently, the inversion results based on the synthetic 3D model show that the inversion of multicomponent data can be er recover the structure of SMS deposits.

Comparing the Computation Times and Memory for the Different Inversion
Both SP inversions were implemented on the PC computer equipped with an Intel Core i7-7820HQ 2.90 GHz CPU. Table 1 illustrates a comparison of the time and memory usage for the different inversions. Both inversions used the same mesh, but the multicomponent data consisted of twice as many data points as the single-component data. Although the memory usage was similar for both cases, the inversion for the multicomponent data took longer to complete.

Data Acquisition
We now recount how we applied our methodology to a real case study. For this purpose, we used the multicomponent electric field collected at Yuhuang hydrothermal field on Southwest Indian ridge during the Research DY58 cruise in 2020. Six Ag/AgCl nonpolarized electrodes were mounted on the Qianlong II AUV (Figure 5a). The speed of the AUV was 0.5-1.0 m/s, and the AUV was maintained at an average altitude of ~50 m above the seafloor. The recorded SP survey line is reported in Figure 5b. Three pairs of non-polarized electrodes (channels 1-6) were used to measure the potential differences: . The position of the electrodes is shown in Figure 6. From the potential differences and the distance between the electrodes, we determined the electric field, , along the measurement survey. Two horizontal components (E34, E56) and one vertical component electric field (E12) were simultaneously collected, and their combination can be used to assess the electrical field in regard to both amplitude and direction. The distance of Electrodes #1 and #2 is 1 m, the distance of Electrode #3 and #4 is 3.5 m (first horizontal component), and E56 (the second horizontal component) is characterized by a spacing of 1.2 m. Compared to Electrodes #3 and #4, which collected the data E34, the electrode dipole consisting of Electrodes #5 and #6 is located far away from the body of the AUV itself; therefore, the signal-to-noise ratio of the electric field E56 is higher than it is for the other components. At the same time, the E34 and E56 are both part of the horizontal electric field. Therefore, the data processing and analysis in horizontal and vertical directions were focused on the horizontal electric field E56 ( Figure 7c) and vertical E12 (Figure 7d). Figure 6. Schematic diagram of AUV-based self-potential system. The numbers on the AUV show the locations of six Ag/AgCl electrodes. Electrodes #1 and #2 were used to measure the vertical electric field E12, Electrodes #3 and #4 were used to measure the horizontal electric field E34, and Electrodes #5 and #6 were used to measure the horizontal electric field E56. The CTD sensor mounted on the AUV was used to measure in the situ conductivity and temperature of seawater along the survey.

Data Processing
The data from the a itude sensors on the AUV showed that the a itude of the AUV was relatively stable above the seafloor during the SP survey, with the pitch angle roughly within ±10° and the angle of cross roll controlled within ±3° [13]. The source of noise in horizontal component data E56 and vertical component data E12 is different. The noise in the data E56 comes from the drift and variable temperature. For the data E12, the distance of Electrodes #1 and #2 is short (approximately 1 m), and they are closer to the propellers, so the noise of the vertical component is mainly from the AUV propellers' rotation. Su et al. [13] showed the denoising flow of horizontal component data E56; here, we first used the same procedure to process the vertical component data E12. Then, the corrected vertical component data eliminated the noise introduced by the motor rotation. The corrected and filtered vertical component (E12) and corrected horizontal component (E56) electric field data displayed in Figure 7a,b were obtained, and the consistent anomalous appeared in the horizontal and vertical components' data at ~GMT 6:00 (where the maximum value of the horizontal component reaches 0.5 mV/m and the maximum value of the vertical component is about −1.8 mV/m. We matched the data with the ultra-short baseline (USBL) data on the AUV, and the 2D horizontal and vertical electric field data used to invert are shown in Figure 7c,d. As for the inversion of multicomponent data in the synthetic model test, the different components were rearranged (i.e., the horizontal component is in the first place, and the vertical component follows behind).

Inversion Results of Multicomponent Data
The 3D inversion of the horizontal and vertical electric field data was performed to obtain the geometry of SMS deposits at the Yuhuang hydrothermal field. A value of 0.1 S/m was chosen as the conductivity of the surrounding rocks below the seafloor based on measurements of laboratory rocks' electrical properties [7]. The in situ conductivity of the seawater was measured by CTD sensors (average conductivity at 3.2 S/m during the survey). The seawater and base rock below the seafloor were dissected with a 3D tensor mesh. In the core area of the SP anomaly, the fine mesh has a cell size of 10 × 10 × 5 m, and the outer padding cells have an increasing width that extends out to the boundaries of the model. There are 532 observations data in the inversion. The observation noise is set to be 5% and the maximum number of iterations is 15 (Figure 8a). The measured electric field and the predicted electric field are shown in Figure 8b. We note that the measured electric field is well reproduced by the inverted source current density distribution. The inverted 3D current density is shown in Figure 9a. Three compact negative-current sources that are associated with ore bodies below the seafloor were recovered. The largest current anomaly extends 80 m below the seafloor, which shows to be about 100 m long in the north-south direction and 100 m long in the east-west direction, and it corresponds well with the previous geological survey (Figure 9b). The fine geometry of the ore was mapped, and we provide key elements for resource evaluation.

Conclusions
Usually, the inversion of self-potential data near the seafloor for the localization and source current density inversion is performed with the measurement of a single component of the electrical field. A synthetic model test was first used to quantify how the inversion results of multicomponent data can be er recover the distribution and amplitude of the source current density caused by ore bodies. Then, our 3D inversion approach was applied to a multicomponent self-potential dataset collected at the Yuhuang hydrothermal field during the Research DY58 oceanic cruise. We obtained the spatial distribution of ore bodies by inverting for the source current density distribution. The advantage of the AUV-based self-potential method is that it overcomes the influence of the seafloor currents in order to stabilize the position of the electrodes along a smooth path. In addition, the simultaneous acquisition of multicomponent self-potential data represents an advantage since this additional information can effectively be used to be er image ore bodies.
The AUV-based SP method and 3D inversion of multicomponent data are expected to play an important role in the future exploration and resource evaluation of SMS deposits. However, current collection methods of SP data have drawbacks. The electrodes are too close to the propeller of the AUV, and the collected signals are subject to electromagnetic interference from the AUV, so we need to make sure that the electrodes are far away from the AUV to further improve the signal-to-noise ratio of the data. Such an approach will also be valuable in collected induced-polarization data that could be used in concert with SP data in order to improve SMS-deposit imaging.
Author Contributions: Z.Z., A.R. and Z.S (Zhaoyang Su). wrote the paper; C.T. was the project leader; Z.S. (Zhigang Shan) was the supervision; Z.Z., J.S., X.D. and J.Z. designed the AUV-SP system; Z.Z., Z.S (Zhaoyang Su) and Z.N. participated in the data collection and processing. All authors have read and agreed to the published version of the manuscript.