TY - JOUR AU - Minshull, Timothy, A AB - SUMMARY Deep-towed geophysical surveys require precise knowledge of navigational parameters such as instrument position and orientation because navigational uncertainties reflect in the data and therefore in the inferred geophysical properties of the subseafloor. We address this issue for the case of electrical conductivity inferred from controlled source electromagnetic data. We show that the data error is laterally variable due to irregular motion during deep towing, but also due to lateral variations in conductivity, including those resulting from topography. To address this variability and quantify the data error prior to inversion, we propose a 2-D perturbation study. Our workflow enables stable and geologically reliable results for multicomponent and multifrequency inversions. An error estimation workflow is presented, which comprises the assessment of navigational uncertainties, perturbation of navigational parameters, and forward modelling of electric field amplitudes for a homogeneous and then a heterogeneous subseafloor conductivity model. Some navigational uncertainties are estimated from variations of direct measurements. Other navigational parameters required for inversion are derived from the measured quantities and their error is calculated by means of error propagation. Some navigational parameters show direct correlation with the measured electric fields. For example, the antenna dip correlates with the vertical electric field and the depth correlates with the horizontal electric field. For the perturbation study each standard deviation is added to the navigational parameters. Forward models are run for each perturbation. Amplitude deviations are summed in quadrature with the stacking error for a total, laterally varying, data error. The error estimation is repeated for a heterogeneous subseafloor model due to the large conductivity range (several orders of magnitude), which affects the forward model. The approach enables us to utilize data from several components (multiple electric fields, frequencies and receivers) in the inversion to constrain the final model and reduce ambiguity. The final model is geologically reasonable, in this case enabling the identification of conductive metal sulphide deposits on the seafloor. Electrical properties, Instrumental noise, Controlled source electromagnetics (CSEM), Marine electromagnetics 1 INTRODUCTION Marine controlled source electromagnetic (CSEM) experiments using a deep-towed electromagnetic field transmitter (e.g. Sinha et al. 1990; Constable 2013) and towed electric field receivers (e.g. Constable et al. 2016) are an effective geophysical tool to map the electrical conductivity up to depths of several hundred metres below the seafloor. The system is commonly applied to detect and delineate shallow small-scale anomalies in the seafloor that are either resistive, such as caused by gas-hydrates (Constable et al. 2016) or conductive, such as caused by seafloor massive sulphides (Gehrmann et al.2019) compared to the surrounding host rock. In areas of relatively flat bathymetry, that is most continental margins, navigational parameters such as towing altitudes and dipole alignments are measurement quantities that remain fairly stable throughout surveying (Myer et al. 2012). As a result, the systematic error is commonly relatively small and assumed to be constant along the profile. Consequently, towed CSEM data have been successfully analysed with a constant relative error that account for the systematic uncertainties of the measurement system, for example, to detect shallow gas hydrate occurrences on the continental shelf that could be correlated with active gas venting (Goswami et al. 2017). Constable et al. (2016) converted navigational uncertainties into associated errors in the measured electrical fields, which, when the conductivity structure does not vary much laterally, can be estimated with a perturbation study of the navigational parameters across a simple 1-D model. When surveying areas with rough seafloor topography and laterally changing conductivity, the 1-D assumption is not adequate. The electric field amplitude and its error depend on the conductivity model. Furthermore, towing the equipment at a (relatively) constant altitude above the seafloor is challenging and involves larger associated measurement errors which are laterally variable due to regular adjustments of the acquisition system. The Trans-Atlantic Geotraverse (TAG) hydrothermal field in 26°N at the Mid-Atlantic Ridge is an example of an area with large seafloor depth variability. The seafloor rises steeply from the neovolcanic zone at about 4.5 km depth to the Eastern Wall at about 2.5 km depth over just ∼10 km. Mound-like structures may reach heights of ∼50 m over diameters of ∼200 m with slopes up to ∼45°. Additionally, ridge-parallel faults and fissures cut and deform the mound-like structures with up to ∼10 m vertical displacement (Petersen & Shipboard Scientific Party 2016). The TAG area is affected by long-term hydrothermal venting leading to alteration of oceanic crust and creation of seafloor massive sulphide (SMS) deposits (Humphris et al. 2015). The electrical properties differ greatly between gabbros, basalt, altered basalt and SMS, and the electrical conductivity is several orders of magnitude larger for SMS than for basalt (Spagnoli et al. 2016). Hence, conductivities will vary laterally not only due to the bathymetry, but also because of the small-scale conductivity anomalies beneath the seabed. In addition, amplitude fluctuations due to navigational uncertainty of the measurement system may result in patchy inversion models that limit geological interpretation. Here, we present a two-stage work flow adapted and improved from an earlier approach by Haroon et al. (2018) to reduce the effects of varying navigational uncertainties along a measured transect. The navigational uncertainty of each measurement parameter is either directly derived from positioning instrumentation or through error propagation theory of the measured parameters through the measurement system (Fig. 1). Subsequently, the derived navigational parameters are used in a 2-D perturbation study that utilizes the bathymetry mapped onto a 2-D grid and a homogeneous seafloor of 10 Ωm to derive error limits attributed to navigational fluctuations along the profile. Variable data errors for the electric fields are obtained by adding the amplitude differences from each parameter in quadrature. A 2-D inversion is conducted for the preliminary navigational uncertainty estimates, and in a second step, the estimates are improved by re-applying the 2-D perturbation study using the final laterally varying conductivity model from the first inversion run. The final conductivity model used for the geological interpretation is derived by re-running the 2-D inversion with the updated error estimates of the second perturbation study. The proposed work flow allows a 2-D inversion of electric field components giving less weight to areas with larger navigational uncertainties (e.g. due to rapid movement) and more weight to components which are less sensitive to the navigational uncertainties. Figure 1. Open in new tabDownload slide Top panel: sketch of the deep-towed electromagnetic experiment across a sulphide mound showing the electric current streamlines (white arrows) at the antenna and the two three-axis electric field receivers. Bottom panel: detailed sketch of sensor position along instrument array with specific legend. Two ultra-short baseline (USBL) transponder beacons were placed at the DASI body and at the end of the array. An altimeter measured DASI’s altitude over the seafloor. A conductivity, temperature and pressure/depth (CTD) sensor was placed at DASI. Pressure was also recorded at the end of the antenna and at the two Vulcan receivers. The two Vulcan receivers additionally recorded instrument heading and tilt. Figure 1. Open in new tabDownload slide Top panel: sketch of the deep-towed electromagnetic experiment across a sulphide mound showing the electric current streamlines (white arrows) at the antenna and the two three-axis electric field receivers. Bottom panel: detailed sketch of sensor position along instrument array with specific legend. Two ultra-short baseline (USBL) transponder beacons were placed at the DASI body and at the end of the array. An altimeter measured DASI’s altitude over the seafloor. A conductivity, temperature and pressure/depth (CTD) sensor was placed at DASI. Pressure was also recorded at the end of the antenna and at the two Vulcan receivers. The two Vulcan receivers additionally recorded instrument heading and tilt. 2 METHODOLOGY In July 2016, deep-towed CSEM profiles were acquired in the TAG hydrothermal field across known SMS targets (Gehrmann et al. 2019). The deep-towed active source instrument (DASI, Sinha et al. 1990) transmitted an electromagnetic (EM) square wave signal on average of about 84 Ampere with a fundamental frequency of 1 Hz through a 50-m-long dipole antenna. Two three-axis electric field receivers (Vulcan, Constable et al. 2016) were towed at 350 and 503 m offset behind the antenna centre (Fig. 1), where the measured inline electric field component is referred to as Ey, the horizontal component orthogonal to the antenna Ex and the vertical Ez. The Ey data for the nearest Vulcan receiver (Vulcan 1) were affected by strong amplitude distortions caused by a damaged electrode and, therefore, excluded from data analysis. Navigational parameters were measured continuously and included DASI’s position with an ultra-short baseline (USBL) transponder, pressure at each instrument and at the end of the antenna, orientation of the Vulcan receivers and DASI’s altitude above the seafloor. Additionally ship parameters such as winch speed, deep-sea cable length, GPS position and speed as well as environmental parameters such as water conductivity and temperature were recorded. Positions of the sensors along the instrument array are shown on the bottom panel of Fig. 1. Based on the approach of Myer et al. (2011), processing of the three-axis electric field data involved the following steps: Applying a Fourier transform on 1-s long non-overlapping time windows of the transmitter and receiver time-series. Calculating the Earth response through deconvolution of the receiver signal by the source signal (division in frequency domain) and normalising the data by the antenna and receiver dipole length. Stacking the Earth response over 30-s-long time intervals (averaging over about 15 m of the seafloor) to derive a mean Earth response and its standard deviation, which we refer to as stacking error. Correcting for the time lag of the source and receiver clocks due to the drift of the Vulcan quartz clocks compared to the GPS triggered source. To infer a conductivity model of the subseafloor we applied the 2-D inversion algorithm MARE2DEM (Key 2016), a linearized algorithm favouring gradual conductivity changes over abrupt changes between model cells. Preliminary 2-D inversion results (shown in Fig. 2) of Vulcan 1 and 2 Ez data, and Vulcan 2 Ey data using a constant error floor of 5 per cent in addition to the stacking error, contain erroneous conductivity structure that is geologically unreasonable. On other profiles the assumption of a constant error floor even prohibits the inversion for multiple components and frequencies from converging and is therefore insufficient. Data residuals (Fig. 2d), defined as the difference between the predicted and measured data for the final conductivity model, which were standardized by dividing by their corresponding error, are laterally heterogeneous. The reasons for this variability are that a constant error floor overestimates the error in some places and underestimates it in others, and that navigational uncertainties affect the data and vary along the profile. Here, we present a work flow to include laterally variable and frequency dependent errors based on navigational uncertainties in electric field amplitude. Unfortunately, an additional time lag with unknown source has prevented us from using phase data. The amplitude, however, is affected more strongly by navigational errors than phase (Myer et al. 2012), and our work flow for amplitude will be suitable for real and imaginary parts of the electric field individually to also encompass phase. Figure 2. Open in new tabDownload slide (a) Final conductivity model for a profile across Southern and Double mound using a constant relative error floor of 5 per cent and a target misfit of 1.0. Colours are faded in areas of lower data sensitivity. Antenna and receiver positions are shown as circles and triangles, respectively. (b) Horizontal, and (c) vertical electric field predicted and observed data with error bars for three frequencies (1–5 Hz) are plotted at antenna-receiver midpoints. (d) Standardized residuals. Figure 2. Open in new tabDownload slide (a) Final conductivity model for a profile across Southern and Double mound using a constant relative error floor of 5 per cent and a target misfit of 1.0. Colours are faded in areas of lower data sensitivity. Antenna and receiver positions are shown as circles and triangles, respectively. (b) Horizontal, and (c) vertical electric field predicted and observed data with error bars for three frequencies (1–5 Hz) are plotted at antenna-receiver midpoints. (d) Standardized residuals. 2.1 Background noise Electric field noise data were collected when the source was switched off. The system was towed at a depth of approximately 2500 m below sea level at a speed of approximately 1.2 kn (∼0.62 m s–1). The amplitude spectra (Fig. 3) were computed for every 60-s long time window in an 1-hour long interval. Following Djanni et al. (2016) the electric field noise En(f) [V m–1|$\sqrt{\rm {Hz}}$|–1] in the frequency domain can be modelled with: $$\begin{eqnarray} E_n(f)= E_T (f) + E_i (f) + V_r(f)/l_r , \end{eqnarray}$$ (1) where ET(f) is the contribution from the motionally induced noise caused by water currents and sudden cable tugging, Ei(f) is the contribution from the naturally occurring magnetotelluric (MT) signals and Vr/lr is the combined voltage noise of electrodes and the instrument electronics divided by the dipole length. These three terms can be classified into two groups: external (ET(f) + Ei(f)) and internal noise (Vr(f)/lr). The MT signal for frequencies higher than 1.2 Hz is mostly attenuated at 2500 m water depth, which is >10 times the skin depth of these frequencies, and is therefore negligible. Motionally induced noise is visible at specific frequency ranges (black arrows in Fig. 3). We observed correlated noise bursts on all the channels at frequencies between 1.5 and 3 Hz (dashed circle in Fig. 3) and noise bursts on some channels centred at about 0.2, 0.4, 3.5 and 7 Hz, respectively. These bursts raise the noise floor by approximately a factor of 4, and are due to the vibration and strumming of the electrodes at the end of each electric dipole (also observed by Constable et al. 2016). Figure 3. Open in new tabDownload slide Amplitude spectrum density for inline Ey, orthogonal to the antenna Ex and vertical Ez component of the electric field for two Vulcans towed at a speed of 1.2 kn (0.62 m s–1) calculated from 1 hr of data when the source was switched off. Noise peaks at certain frequencies indicated by black arrows are identified as strumming and vibration noise. The dashed rectangle marks a typical range of operation for marine CSEM. We have used 1, 3 and 5 Hz only. Figure 3. Open in new tabDownload slide Amplitude spectrum density for inline Ey, orthogonal to the antenna Ex and vertical Ez component of the electric field for two Vulcans towed at a speed of 1.2 kn (0.62 m s–1) calculated from 1 hr of data when the source was switched off. Noise peaks at certain frequencies indicated by black arrows are identified as strumming and vibration noise. The dashed rectangle marks a typical range of operation for marine CSEM. We have used 1, 3 and 5 Hz only. In general, the noise term is mainly made up of the third term in eq. (1). The Ey component in Vulcan 2 (orange curve in Fig. 3) has the lowest noise level compared to the other components. The amplitude levels for Ex and Ez are twice (high frequency band) to ten times (low frequency band) higher than for Ey. This difference represents an improvement of the signal-to-noise ratio proportional to the dipole length (Ey dipole is twice as long). At frequencies above ∼3 Hz, the spectrum is flat for all channels. We can conclude that the flat spectrum above ∼3 Hz is due to white noise while the spectrum below ∼3 Hz is due to pink noise (amplitude falling off by f−1) from the instrument electronics and electrodes. 2.2 Navigational uncertainty analysis The error analysis involves two methods illustrated below. For parameters for which direct measurements exist such as pressure and instrument tilt, their mean and standard deviations are calculated over the same 30-s-long time window that was used for data stacking. For parameters which are functions depending on other parameters such as the position and orientation of the instruments along the profile used for the inversion, their error is calculated through error propagation theory. 2.2.1 Averaging over stacking time window Parameters such as the pitch and the roll of the Vulcans (from the tilt sensor, shown in Fig. A1) are averaged over the stacking window of 30 s giving a mean value and its standard deviation. The depths of DASI, the end of the antenna and the two Vulcans are converted from pressure assuming hydrostatic conditions at mid-latitude oceans (Saunders 1981). The pressure measurement uncertainty is neglected as it is much smaller than the depth variation in the stacking time window due to the constant readjustment of the instrument height over rough terrain and a bias of the pressure sensor in over 3 km water depth. A second estimate of the depth of DASI comes from the seafloor depth directly underneath DASI minus the distance to the seafloor from the altimeter reading, which for the TAG data set differs from the converted pressure reading by about 5 m. We adjust the depth and add a base standard deviation of 2 m (Fig. 4). The instrument’s position is estimated using a USBL transponder installed on DASI. The standard deviation along profile and perpendicular to the profile (in a rotated coordinate system) is estimated over the stacking window. Figure 4. Open in new tabDownload slide Amplitude, antenna dip and altitude (estimated from instrument and seafloor depth) at antenna-receiver midpoint. Single values of the Tx dip and altitude are plotted at the midpoint of the Tx-V1 location. The choice of position is designed to illustrate the correlation between for example the dip and the Ez amplitude. The inset on the right shows the cross-correlation between Vulcan 2 Ez amplitude for 5 Hz and dip, normalized by the square root of the product of the autocorrelations of both time series at zero lag. Plotting amplitude and navigational parameters at the same location is equivalent to plotting them over time. Figure 4. Open in new tabDownload slide Amplitude, antenna dip and altitude (estimated from instrument and seafloor depth) at antenna-receiver midpoint. Single values of the Tx dip and altitude are plotted at the midpoint of the Tx-V1 location. The choice of position is designed to illustrate the correlation between for example the dip and the Ez amplitude. The inset on the right shows the cross-correlation between Vulcan 2 Ez amplitude for 5 Hz and dip, normalized by the square root of the product of the autocorrelations of both time series at zero lag. Plotting amplitude and navigational parameters at the same location is equivalent to plotting them over time. 2.2.2 Error propagation Other navigational parameters are not directly measured and their uncertainty is a result of error propagation from independent parameters (e.g. Taylor 1997). The standard deviation ∂s of any parameter s = f(a, b, ...) can be estimated by adding individual errors times the derivative of the parameter in quadrature (neglecting correlations between parameters). $$\begin{eqnarray} \delta s= \sqrt{\left(\frac{\partial f}{\partial a} \delta a\right)^2+\left(\frac{\partial f}{\partial b} \delta b\right)^2+...} . \end{eqnarray}$$ (2) In our case, the dip of the antenna γ is computed by differencing the depth as measured at the two ends (DASI |$z_\rm {\small {DASI}}$| and the end of the antenna |$z_\rm {\small {EOA}}$|⁠) divided by the distance l between them. $$\begin{eqnarray} \gamma =\arcsin \left(\frac{z_\rm {\small {DASI}}-z_\rm {\small {EOA}}}{l}\right) = \arcsin (A) . \end{eqnarray}$$ (3) The error, or standard deviation, in the dip has to be propagated (eq. 2) from the error in the two depths and distance, which are presumed to be independent, and is calculated using $$\begin{eqnarray} \delta \gamma = \frac{|A|}{\sqrt{1-A^2}}\sqrt{\left(\frac{\delta z^2_\rm {\small {DASI}}+ {\delta }z^2_\rm {\small {EOA}}}{(z_\rm {\small {DASI}}-z_\rm {\small {EOA}})^2}+\frac{\delta l^2}{l^2}\right)}. \end{eqnarray}$$ (4) Similarly to the dip we calculate the angle between the end of the antenna and the first Vulcan as well as between the first and the second Vulcan to then estimate the horizontal distances between the instruments and their standard deviations (see Appendix A1 and A2). The orientation in the horizontal plane requires the heading of the Vulcans. The heading measurement is biased by the battery pack (Constable et al. 2016). A polynomial correction function was constructed which is approximately sinusoidal and based on matching the Vulcan heading with DASI’s heading (see Appendix A3). The heading uncertainty is the sum of the standard deviation from the stacking time window and the polynomial curve fit with the assumed heading (given Vulcan and DASI heading are similar). The absolute position of each instrument is then calculated using the estimated headings. The profile for the 2-D conductivity model is set along a straight line which is fit between the Vulcan and antenna positions. New Vulcan and antenna positions are now calculated with respect to a coordinate system along the profile direction and perpendicular to it. Navigational uncertainties for the Vulcan pitch and roll (often less than one degree), which are averaged over the time window, are relatively small, while they are larger for the Vulcan heading (few degrees). The horizontal position of the instruments also has a larger uncertainty of 20–40 m due to error propagation of the uncertainties in the USBL position, heading and antenna dip. Some parameters such as altitude/depth and antenna dip already show a direct correlation with the amplitude data (Fig. 4). The horizontal electric field amplitude decreases with increasing altitude. The vertical electric field varies with the antenna dip. In fact, the dependency of the amplitudes on the navigational parameters precludes qualitative interpretation of subseafloor conductivities directly from the data. 2.3 Integration of navigational uncertainty into the data analysis Data analysis with a 2-D inversion algorithm (e.g. MARE2DEM, Key 2016) requires a 2-D starting model, the data and associated error, and the transmitter and receiver geometry. The latter is calculated through the procedure illustrated above. The starting model is built using the seafloor depth for the 2-D profile interpolated from high-resolution (2 m) bathymetry data obtained from GEOMAR’s Abyss system (Petersen & Shipboard Scientific Party 2016). The conductivity for air (10−12 S m–1) and for several layers of water (5 S m–1 at the surface to 3.28 S m–1 at DASI’s towing altitude) are fixed in the model. The initial seafloor conductivity is chosen to be 0.1 S m–1 as an average value for oceanic crust (Evans & Everett 1994). The root-mean-square misfit between the observed (d) and predicted (with forward operator F on model m) data divided by the individual data error ϵ is calculated with $$\begin{eqnarray} \chi = \sqrt{\frac{1}{N}\sum _{i}^{N}{\frac{(d_{i} -F_{i}({\bf m}) )^2}{\epsilon _{i} ^2}}}. \end{eqnarray}$$ (5) The inversion algorithm compares χ to a chosen target misfit and the inversion stops when the target misfit is reached. A target misfit of 1.0 means that the inversion algorithm fits the observed data within the data error and a smaller target misfit allows the algorithm to fit the data better than the data error. A rigorous choice of the data error for the inversion is, therefore, crucial. To estimate a data error that incorporates variability dependent on position and the conductivity model ϵ([x, y], m), we have developed the work flow illustrated in Fig. 5. First, the standard deviations of the navigational parameters are added to the parameters individually. For each perturbed parameter a forward model is run to calculate synthetic data. Amplitude differences in percentage are then calculated by comparing the data to the non-perturbed case (Fig. 6). For the horizontal electric field Ey the largest amplitude differences are for variations in instrument depth, horizontal position and antenna dip. Amplitude differences in the vertical electric field Ez are largest for variations in antenna dip and horizontal position. The standard deviation from the stacked time series is orders of magnitude larger for Ez than Ey which is partly related to the smaller dipole length and to Ez being more sensitive to changes in dip and pitch. Figure 5. Open in new tabDownload slide Work flow to integrate navigational uncertainties in the error model for CSEM data inversion: First, the navigational uncertainties are calculated from direct measurements. A 2-D conductivity model is constructed including a seawater conductivity depth profile from the CTD and a homogeneous subseafloor (Model 1). In the next step (Forward), synthetic data are modelled. The amplitude error is then estimated summing up amplitude differences from modelling and standard deviation from stacking in quadrature. The error model with individual standard deviations for each datum are then fed into the inversion algorithm to obtain an optimal heterogeneous conductivity model (Model 2). The forward modelling, error estimation and inverse modelling is then repeated for the heterogeneous model to obtain the final conductivity model (Model 3). Figure 5. Open in new tabDownload slide Work flow to integrate navigational uncertainties in the error model for CSEM data inversion: First, the navigational uncertainties are calculated from direct measurements. A 2-D conductivity model is constructed including a seawater conductivity depth profile from the CTD and a homogeneous subseafloor (Model 1). In the next step (Forward), synthetic data are modelled. The amplitude error is then estimated summing up amplitude differences from modelling and standard deviation from stacking in quadrature. The error model with individual standard deviations for each datum are then fed into the inversion algorithm to obtain an optimal heterogeneous conductivity model (Model 2). The forward modelling, error estimation and inverse modelling is then repeated for the heterogeneous model to obtain the final conductivity model (Model 3). Figure 6. Open in new tabDownload slide Amplitude differences for synthetic data of Vulcan 2 for Ey (left-hand panel) and Ez (right-hand panel) at 1, 5 and 9 Hz (top to bottom) for changing each navigational parameter with their standard deviation for the heterogeneous seafloor model (Model 2 in workflow, Fig. 5). Navigational parameters are (see also Fig. A1): Receiver azimuth θRx, pitch β, roll α, antenna azimuth θTx, dip γ, array horizontal position y and array depth z. The total relative error is the sum of the navigational errors and the stacking error σst in quadrature. Figure 6. Open in new tabDownload slide Amplitude differences for synthetic data of Vulcan 2 for Ey (left-hand panel) and Ez (right-hand panel) at 1, 5 and 9 Hz (top to bottom) for changing each navigational parameter with their standard deviation for the heterogeneous seafloor model (Model 2 in workflow, Fig. 5). Navigational parameters are (see also Fig. A1): Receiver azimuth θRx, pitch β, roll α, antenna azimuth θTx, dip γ, array horizontal position y and array depth z. The total relative error is the sum of the navigational errors and the stacking error σst in quadrature. Second, the individual amplitude differences are added up together with the stacking error in quadrature to estimate a total relative difference that varies for each receiver position and frequency. Total relative errors increase with frequency to values larger than 10 per cent. Including data points with large errors will allow the target misfit to be reached even when some data ranges with smaller data errors (low frequency data) are not fit yet, which results in an increase in ambiguity and in this case, due to the chosen regularization (Occam’s razor), in a smooth/homogeneous model. To allow data ranges with small data errors to have a higher weight in the inversion frequencies larger than 5 Hz were excluded. The total amplitude difference is then converted to absolute values using the synthetic data and a first inversion is started with the initial model, geometry and the newly estimated data error. The inversion converges to an optimal model (Model 2 in Fig. 5) with large conductivity variability given the conductive SMS deposits and the resistive oceanic crust. Due to the complexity of the model and conductivity varying over a couple of magnitudes the error estimation is repeated including forward modelling and adding up individual amplitude differences. The new absolute data error is included in a second inversion run. Residual errors [di − Fi(m)] for Ey tend to be smaller than the navigational errors (Fig. 7) indicating that the errors may be overestimated. The mismatch between residual and navigational errors is smaller for larger frequencies for Ey, and absent for Ez. The navigational errors match better for some profiles (example shown in Fig. 3 in Gehrmann et al. 2019) but are overestimated for others. Reasons for this could be that navigational effects on amplitude are not necessarily additive, but may cancel each other out. The relationship between uncertainties in navigational parameters and amplitude is more complex and also depends on the subsurface model. Adding up the standard deviations in quadrature therefore yields an upper error bound. Another criterion for convergence is the residual statistics. Ideally the residuals should be Gaussian distributed [intrinsic to eq. (5)], but navigational errors and 3-D effects (Haroon et al. 2018) introduce an unavoidable bias. Therefore, we have step-wise reduced the target misfit until the inversion won’t be able to reach the target misfit any more. Figure 7. Open in new tabDownload slide Standard error distribution for Vulcan 2 for Ey (left-hand panel) and Ez (right-hand panel) at 1 Hz (top panel) and 5 Hz (bottom panel) from the standard deviation derived from stacking over a 30-s-long time window, the residuals between observed and predicted data for the final model, and from navigational error which depends on the position and the model, and which is estimated using the workflow described in Section 2.3. Figure 7. Open in new tabDownload slide Standard error distribution for Vulcan 2 for Ey (left-hand panel) and Ez (right-hand panel) at 1 Hz (top panel) and 5 Hz (bottom panel) from the standard deviation derived from stacking over a 30-s-long time window, the residuals between observed and predicted data for the final model, and from navigational error which depends on the position and the model, and which is estimated using the workflow described in Section 2.3. 3 RESULTS AND DISCUSSION A profile across two SMS deposits (Southern and Double mound, Gehrmann et al. 2019) has been inverted with a target misfit of 0.8 using our proposed work flow (Fig. 8). For comparison, the profile was also inverted using a constant minimum relative error floor of 5 per cent (Fig. 2). Data points with larger errors (>3800 m along profile) had to be excluded from the inversion to be able to converge by reaching a target misfit of 1.0. Our workflow therefore reveals areas with large uncertainties which are then given less weight in the inversion (Figs 8b and c) due to the increased data error. Biases due to systematic navigational uncertainties or 3-D topography (Haroon et al. 2018) exist for both inversion results (compare data fit in Figs 2 and 8b and c). Both inversion models show conductive features greater than 1 S m–1 at the two mounds which can be explained by seafloor massive sulphide occurrences (more conductive than seawater-saturated sediments, e.g. Spagnoli et al. 2016; Müller et al. 2018; Gehrmannet al.2019). The comparison also highlights the advantages of the proposed work flow, as additional conductivity features that are introduced by overfitting amplitude fluctuations (Fig. 2) are less prominent in the final inversion model (Fig. 8). Figure 8. Open in new tabDownload slide (a) Final conductivity model for profile across Southern and Double mound using a laterally variable data error estimated from navigational uncertainties and a target misfit of 0.8. Colours are faded in areas of lower data sensitivity. Antenna and receiver positions are shown as circles and triangles respectively. (b) Horizontal and (c) vertical electric field predicted and observed data with error bars for three frequencies (1–5 Hz) are plotted at antenna-receiver midpoints. (d) Standardized residuals. Figure 8. Open in new tabDownload slide (a) Final conductivity model for profile across Southern and Double mound using a laterally variable data error estimated from navigational uncertainties and a target misfit of 0.8. Colours are faded in areas of lower data sensitivity. Antenna and receiver positions are shown as circles and triangles respectively. (b) Horizontal and (c) vertical electric field predicted and observed data with error bars for three frequencies (1–5 Hz) are plotted at antenna-receiver midpoints. (d) Standardized residuals. Overall, the electric fields at these short offsets are quite sensitive to navigational changes and little can be derived about the subsurface conductivity structure by means of data inspection. An inversion is necessary to make robust estimates of the subseafloor conductivity structure. For example, the horizontal electric field decreases in amplitude with increasing height above the seafloor. Hence, an amplitude decrease of the inline electric field which correlates with an altitude increase is probably not related to a geological structure in the seafloor. Similarly, the vertical electric field component is correlated with the antenna dip, which changes due to instrument depth changes (Fig. 4). A few metres uncertainty in height above the seafloor is the largest error source in the Ey amplitude while Ez shows a dip dependency that is largest for 1 Hz and reduces for higher harmonics, which have a larger stacking error (Fig. 6). In fact, stacking errors become a larger component of the residual errors with increasing frequency (see 1 Hz against 5 Hz error histogram in Fig. 7). Overall, the navigational uncertainties for Ey are overestimated, especially for low frequencies. In comparison, the navigational uncertainty estimate for Ez seems suitable as the navigational errors match the residual errors (Fig. 7). We have not found a global navigational uncertainty estimate which works for every component, frequency or profile in a way that the target misfit becomes one (the predicted and observed data are fit within the data error). There are too many dependencies and therefore the weights in the inversion change for each profile. Instead we find that a slight overestimation of the navigational error and a reduction of the target misfit below 1.0 allows each component to be fit without introducing erroneous structure in the final inversion model. The main indicator for the convergence of the inversion is the residual fit (Fig. 8d). The ideal results are random and Gaussian distributed standardized residuals, which cannot be achieved here for all components due to the uncertainties resulting from navigation and 3-D effects. Including our new data error model results in residuals with little systematic variation along the profile, and the good data fit validates the navigation assumptions. The final conductivity models are geologically reasonable and consistent with interpretations of other geophysical data (Gehrmann et al. 2019). 4 CONCLUSIONS Deep-tow surveying in great depth over rough terrain causes unpredictable instrument motion due to changing currents and rapid instrument movements. In general, it is more difficult to control the towed instrument from the vessel (for example, through the cable winch speed) than if towed at the surface or in shallow water. Moreover, all navigational parameters, such as depth derived from pressure readings or absolute position from acoustic transponder readings, have additional uncertainties and are less reliable at great depth. Due to the seafloor topography and associated adjustments made to the measurement system to prevent seafloor contact, navigational uncertainties vary along the survey profile and cannot be accounted for using simple 1-D perturbation studies. Erroneous navigational parameters translate directly into systematic uncertainty. Assuming a constant error along the profile, may prevent the inversion from converging and can introduce erroneous conductivity structure in the final inversion model. In our example application, some parts of the profile have much larger errors than others. The effects on amplitude also depend on the topography and the conductivity structure of the subseafloor. Therefore, it is important to estimate the influences of navigational uncertainties in a rigorous 2-D perturbation analysis. We have presented a work flow to estimate the uncertainty of direct measured navigational parameters by averaging over the stacking time window, and how these uncertainties propagate into required navigational parameters for the inversion algorithm. Dependencies between certain navigational parameters and electric field data appear during the process. For example, the horizontal electric field data are sensitive to the instrument’s altitude above the seafloor. In contrast, the vertical field data are generally more sensitive to changes in transmitting antenna dip. Additionally, this sensitivity is frequency dependent, so that larger frequencies are less sensitive to the dip and instead the background noise/stacking error dominates. The dependencies of the electric field components on different navigational parameters and stacking error implies a weighting in the inversion that cannot be avoided. Our estimates of navigational uncertainties and summation of resulting data uncertainties with stacking errors in quadrature enables convergence of a multicomponent, multifrequency inversion, reaching an acceptable data misfit of 1.0 and for some profiles even below 1.0 (indicating that data errors are slightly overestimated). The general aim of our approach is to reach random and Gaussian distributed standardized residuals, but some areas express unavoidable bias. In part, this bias can be caused by 3-D effects that we do not consider in the scope of this paper, but which have been shown to cause erroneous model structure and data biases (Haroon et al. 2018). The inversion algorithm assumes that errors are random and Gaussian, but is stabilized with the smoothness regularization. To validate our results, we have rigorously compared them with other geophysical data and geological prior knowledge (Gehrmann et al.2019). Our final inversion results are geologically meaningful since occurrences of seafloor massive sulphides appear at the expected sites. The suggested approach limits additional erroneous structure while enabling us to include as many data components as possible and therefore information that is sensitive to different model regions and depths. ACKNOWLEDGEMENTS We would like to thank the European Commission for funding the Framework 7 project Blue Mining (Grant number 604500). We thank all cruise participants and crew from M127 and JC138 for their support in data acquisition especially Ian Tan and Laurence North for their engineering expertise, Sebastian Hölz for his real-time navigation script, and Eric Attias and chief scientist Bram Murton for their proactive involvement. We thank Karen Weitemeyer and Steven Constable for their advice, and Steven for internally reviewing the manuscript. We also thank editor Ute Weckmann, reviewer Rob Evans and one anonymous reviewer for their constructive feedback. We thank Kerry Key for the inversion code MARE2DEM and supporting matlab scripts as well as David Myer for processing routines. Tim Minshull was supported by a Wolfson Research Merit award. Amir Haroon was in part funded by the German Research Foundation, Grant number 389727048. Data are available from Gehrmann (2019). REFERENCES Berndt C. et al. . , 2017 . RV MARIA S. MERIAN Fahrtbericht / Cruise Report MSM63 - PERMO, Southampton Southampton (U.K.) 29.04.-25.05.2017. GEOMAR Report N.Ser. 037 . OpenURL Placeholder Text WorldCat Constable S. , 2013 . Review paper: Instrumentation for marine magnetotelluric and controlled source electromagnetic sounding , Geophys. Prospect. , 61 , 505 – 532 . 10.1111/j.1365-2478.2012.01117.x Google Scholar Crossref Search ADS WorldCat Crossref Constable S. , Kannberg P.K. , Weitemeyer K. , 2016 . Vulcan: a deep-towed CSEM receiver , Geochem. Geophys. Geosyst. , 17 ( 3 ), 1042 – 1064 . 10.1002/2015GC006174 Google Scholar Crossref Search ADS WorldCat Crossref Djanni A.T. , Ziolkowski A. , Wright D. , 2016 . Electromagnetic induction noise in a towed electromagnetic streamer , Geophysics , 81(3) , E187 – E199 . 10.1190/geo2014-0597.1 Google Scholar Crossref Search ADS WorldCat Crossref Evans R.L. , Everett M.E. , 1994 . Discrimination of hydrothermal mound structures using transient electromagnetic methods , Geophys. Res. Lett. , 21 ( 6 ), 501 – 504 . Google Scholar Crossref Search ADS WorldCat Gehrmann R.A.S. , North L.J. , Graber S. , Szitkar F. , Petersen S. , Minshull T.A. , Murton B.J. , 2019 . Marine mineral exploration with controlled-source electromagnetics at the TAG hydrothermal field, 26N Mid-Atlantic Ridge , Geophysical Research Letters , 46 , 11 , 5808 – 5816 . 10.1594/PANGAEA.899073 Google Scholar Crossref Search ADS WorldCat Crossref Goswami B.K. , Weitemeyer K.A. , Bünz S. , Minshull T.A. , Westbrook G.K. , Ker S. , Sinha M.C. , 2017 . Variations in pockmark composition at the Vestnesa Ridge: insights from marine controlled source electromagnetic and seismic data , Geochem. Geophys. Geosyst. , 18 ( 3 ), 1111 – 1125 . 10.1002/2016GC006700 Google Scholar Crossref Search ADS WorldCat Crossref Haroon A. , Hölz S. , Gehrmann R.A. , Attias E. , Jegen M. , Minshull T.A. , Murton B.J. , 2018 . Marine dipole-dipole controlled source electromagnetic and coincident-loop transient electromagnetic experiments to detect seafloor massive sulphides: effects of three-dimensional bathymetry , Geophys. J. Int. , 215 ( 3 ), 2156 – 2171 . 10.1093/gji/ggy398 Google Scholar Crossref Search ADS WorldCat Crossref Humphris S.E. , Tivey M.K. , Tivey M.A. , 2015 . The Trans-Atlantic Geotraverse hydrothermal field: a hydrothermal system on an active detachment fault , Deep Sea Res. Part II: Topical Stud. Oceanogr. , 121 , 8 – 16 . 10.1016/j.dsr2.2015.02.015 Google Scholar Crossref Search ADS WorldCat Crossref Key K. , 2016 . MARE2DEM:A 2-D inversion code for controlled-source electromagnetic and magnetotelluric data , Geophys. J. Int. , 207 ( 1 ), 571 – 588 . 10.1093/gji/ggw290 Google Scholar Crossref Search ADS WorldCat Crossref Key K. , Lockwood A. , 2010 . Determining the orientation of marine CSEM receivers using orthogonal Procrustes rotation analysis , Geophysics , 75 ( 3 ), F63 – F70 . 10.1190/1.3378765 Google Scholar Crossref Search ADS WorldCat Crossref Müller H. , Schwalenberg K. , Reeck K. , Barckhausen U. , Schwarz-Schampera U. , Hilgenfeldt C. , von Dobeneck T. , 2018 . Mapping seafloor massive sulfides with the Golden Eye frequency-domain EM profiler , First Break , 36 ( 11 ), 61 – 67 . OpenURL Placeholder Text WorldCat Myer D. , Constable S. , Key K. , 2011 . Broad-band waveforms and robust processing for marine CSEM surveys , Geophys. J. Int. , 184 , 689 – 698 . 10.1111/j.1365-246X.2010.04887.x Google Scholar Crossref Search ADS WorldCat Crossref Myer D. , Constable S. , Key K. , Glinsky M.E. , Liu G. , 2012 . Marine CSEM of the Scarborough gas field, Part 1: Experimental design and data uncertainty , Geophysics , 77 ( 4 ), E281 – E299 . 10.1190/geo2011-0380.1 Google Scholar Crossref Search ADS WorldCat Crossref Petersen S. , Shipboard Scientific Party , 2016 . RV METEOR Cruise Report M127, Metal fluxes and Resource Potential at the slow-spreading TAG Mid-Ocean Ridge Segment (26oN, MAR) - Blue Mining Sea, Bridgetown (Barbados) - Ponta Delgada (Portugal), 25.05.-28.06.2016 , Tech. rep., GEOMAR, Helmholtz Centre for Ocean Research Kiel . OpenURL Placeholder Text WorldCat Saunders P.M. , 1981 . Practical conversion of pressure to depth , J. Phys. Oceanogr. , 11 ( 4 ), 573 – 574 . Google Scholar Crossref Search ADS WorldCat Sinha M. , Patel P. , Unsworth M. , Owen T. , MacCormack M. , 1990 . An active source electromagnetic sounding system for marine use , Mar. Geophys. Res. , 12 ( 1-2 ), 59 – 68 . Google Scholar Crossref Search ADS WorldCat Spagnoli G. , Hannington M.D. , Bairlein K. , Hördt A. , Jegen M. , Petersen S. , Laurila T. , 2016 . Electrical properties of seafloor massive sulfides , Geo-Mar. Lett. , 36 ( 3 ), 235 – 245 . 10.1007/s00367-016-0439-5 Google Scholar Crossref Search ADS WorldCat Crossref Taylor J.R. , 1997 . An Introduction to Error Analysis—The Study of Uncertainities in Physical Measurements , 2nd edn, University Science Books . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC APPENDIX: NAVIGATIONAL PARAMETERS AND THEIR UNCERTAINTY To estimate the location of each instrument, the angles from the horizontal to the direct path between the instruments, the horizontal distance and the orientation of the instruments were required. Each of these parameters is a combination of other parameters whose errors propagate following eq. (2). Here, we go through each parameter in detail. A1 Angles We calculate the angle of a right triangle beneath the Vulcan receivers through the trigonometric ratio between the opposite, the depth difference between the pressure sensors of two adjacent instruments, and the hypotenuse, the direct distance between them. Analogous to the antenna dip (eq. 3) the angles become $$\begin{eqnarray} \gamma _{\rm v _1}=\arcsin \left(\frac{z_\rm {\small {EOA}}-z_{\rm v _1}}{l_{\rm v _1}-l_\rm {\small {EOA}}}\right) = \arcsin (A) , \end{eqnarray}$$ (A1) and $$\begin{eqnarray} \gamma _{\rm V_2}=\arcsin \left(\frac{z_{\rm v _1}-z_{\rm v _2}}{l_{\rm v _2}-l_{\rm v _1}}\right), \end{eqnarray}$$ (A2) where |$l_\rm {\small {EOA}}$| is equivalent to l in eq. (3), and |$l_{\rm v _2}$| and |$l_{\rm v _1}$| are the distances between the Vulcans and DASI. The angle uncertainty becomes $$\begin{eqnarray} \delta \gamma _{\rm V_1} = \frac{|A|}{\sqrt{1-A^2}}\sqrt{\left(\frac{\delta z^2_\rm {\small {EOA}}+ {\delta }z^2_{\rm V_1}}{(z_\rm {\small {EOA}}-z_{\rm V_1})^2}+\frac{\delta l_{\rm V_1}^2+\delta l_\rm {\small {EOA}}^2}{(l_{\rm V_1}-l_\rm {\small {EOA}})^2}\right)}. \end{eqnarray}$$ (A3) Here, we refer to pitch as β in Fig. A1. One may suggest that the horizontal angles are similar to the measured pitch at the Vulcans. The trend is quite similar but the values may differ by a few degrees without showing a systematic offset over several profiles. Therefore, the difference is likely explained by the pull of the instruments and the parachute in different directions which affect the measured pitch more strongly than the horizontal angles. Figure A1. Open in new tabDownload slide Sketch of Vulcan orientation: Heading θ, pitch β and roll α angles (adapted from Constable et al. 2016; Key & Lockwood 2010, courtesy of the SIO Marine EM Laboratory). Note: the pitch and roll (positive downward and anticlockwise) are opposite direction than recorded. Figure A1. Open in new tabDownload slide Sketch of Vulcan orientation: Heading θ, pitch β and roll α angles (adapted from Constable et al. 2016; Key & Lockwood 2010, courtesy of the SIO Marine EM Laboratory). Note: the pitch and roll (positive downward and anticlockwise) are opposite direction than recorded. A2 Horizontal distances and antenna depth The horizontal distances between instruments are defined as the distances when the instruments are projected to the same depth. These are calculated using the cosine of the respective angle (which was calculated in Appendix A1). First, we sum up the horizontal distances between instruments, but because the directional vectors of the instruments are not the same, that is the instruments are not in a straight line, we later subtract the horizontal distances again to only use the horizontal distance between adjacent instruments to calculate the absolute position (Appendix A4). For example, the horizontal distance h of the second Vulcan is $$\begin{eqnarray} h_{\rm V_2} = h_{\rm V_1} + (l_{\rm V_2}-l_{\rm V_1})\cos (\gamma _{\rm V_2}). \end{eqnarray}$$ (A4) The uncertainty becomes $$\begin{eqnarray} &&\delta h_{\rm V_2} \\ && = \sqrt{\delta h_{\rm V_1}^2 + (\delta l_{\rm V_2}^2+\delta l_{\rm v _1}^2)\cos ^2(\gamma _{\rm v _2}) + (l_{\rm v _2}-l_{\rm v _1})^2\sin (\gamma _{\rm v _2})^2\partial {\delta \gamma _{\rm v _2}}}. \\ \end{eqnarray}$$ (A5) The antenna dip can be used to calculate the antenna midpoint depth |$z_\rm {T{\small x}}$| with $$\begin{eqnarray} z_\rm {T{small X}}=z_\rm {\small {DASI}}- l_\rm {T{small X}}\sin (\gamma ) \end{eqnarray}$$ (A6) and its corresponding standard deviation $$\begin{eqnarray} \delta z_\rm {T{small X}}= \sqrt{\delta z_\rm {\small {DASI}}^2 + \delta l_\rm {T{small X}}^2 \sin ^2(\gamma ) + l_\rm {T{small X}}^2 \cos ^2(\gamma )\delta \gamma ^2 }. \end{eqnarray}$$ (A7) A3 Vulcan heading As described in Constable et al. (2016) the reading of the internal compass has a heading dependent bias due to the magnetic field of the battery pack. We have attempted to calibrate the heading of the Vulcans with the heading of DASI obtaining calibration functions shown in Fig. A2 which are 7th order polynomials that resemble sinusoidal functions. In a survey in 2017 external compasses were mounted in the Vulcan nose (Berndt et al. 2017). The shape of the correction functions are very similar, supporting the polynomial function. The amplitudes, though, are quite different, suggesting that the calibration function is survey dependent. The standard deviation for the heading is obtained from the difference between DASI’s heading minus the Vulcan internal compass reading and the polynomial. The heading of the antenna mid point is taken from a linear relationship between the heading of DASI and the first Vulcan and the distance between them. Figure A2. Open in new tabDownload slide Heading calibration polynomial. Left-hand panel: heading correction from comparison of DASI’s heading with the internal compass reading for all six profiles with four different orientations. Right-hand panel: heading correction for data set acquired one year later from comparison of external and internal compass. Figure A2. Open in new tabDownload slide Heading calibration polynomial. Left-hand panel: heading correction from comparison of DASI’s heading with the internal compass reading for all six profiles with four different orientations. Right-hand panel: heading correction for data set acquired one year later from comparison of external and internal compass. A4 Absolute positions Knowing the heading and the horizontal distance the absolute position in UTM coordinates (Universal Transverse Mercator projection with zone 23 R) can be calculated using trigonometry again. For example, for the antenna mid point: If the heading θ is between 90° and 180°, the angle between the horizontal distance and the direct distance becomes φ = θ − 90. The x and y coordinates in the UTM system therefore become $$\begin{eqnarray} \begin{aligned} {\begin{bmatrix}x\\ y \end{bmatrix}} &= &{\begin{bmatrix}x_\rm {\small {DASI}}-h_\rm {T{small X}}\cos (\varphi )\\ y_\rm {\small {DASI}}+h_\rm {T{small X}}\sin (\varphi ) \end{bmatrix}} \end{aligned} \end{eqnarray}$$ (A8) and the standard deviation is derived similarly to eq. (A7). The position of the Vulcans are calculated with their headings and relative to the instrument in front of them using the difference in the horizontal positions. A5 Profile selection A representative profile line (y = mx + b) is now chosen using a linear regression between the Vulcans’ and DASI’s UTM coordinates. This best-fitting line is used to extract the bathymetry and build a 2-D model for the inversion algorithm. The UTM positions are now converted into x and y position in respect to the line by calculating the nearest point (x, y) for each instrument on the line. $$\begin{eqnarray} \begin{aligned} {\begin{bmatrix}x\\ y \end{bmatrix}} &= \frac{1}{1+m^2}&{\begin{bmatrix}x_\rm {T{small X}}+ my_\rm {T{small X}}-mb\\ mx_\rm {T{small X}}+ m^2y_\rm {T{small X}}+b \end{bmatrix}} .\end{aligned} \end{eqnarray}$$ (A9) If neglecting the standard deviation of the linear regression fit, the standard deviation for the point on the line is a combination of the standard deviation of the instrument positions and the slope of the line m. $$\begin{eqnarray} \begin{aligned} {\begin{bmatrix}\delta x\\ \delta y \end{bmatrix}} &= \frac{1}{1+m^2}&{\begin{bmatrix}1\\ |m| \end{bmatrix}} &\sqrt{(}\delta x_\rm {T{small X}}^2 + (m\delta y_\rm {T{small X}})^2) . \end{aligned} \end{eqnarray}$$ (A10) The estimate in eq. (A10) may lead to an underestimate. For our study we added the standard deviation of the linear regression fit to eq. (A10), but did not use the error propagation theory from eq. (A9) as it results in an unreasonable high estimate. The x and y position along the profile are now calculated using the point on the line and the defined start coordinate of the profile using Pythagoras’ theorem. © The Author(s) 2019. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Seafloor massive sulphide exploration using deep-towed controlled source electromagnetics: navigational uncertainties JF - Geophysical Journal International DO - 10.1093/gji/ggz513 DA - 2020-02-01 UR - https://www.deepdyve.com/lp/oxford-university-press/seafloor-massive-sulphide-exploration-using-deep-towed-controlled-vPcNicvYPj SP - 1215 VL - 220 IS - 2 DP - DeepDyve ER -