Geostatistical Seismic Inversion for Temperature and Salinity in the Madeira Abyssal Plain

A two-dimensional multichannel seismic reflection profile acquired in the Madeira Abyssal Plain during June 2016 was used in a modeling workflow comprising seismic oceanography processing, geostatistical inversion and Bayesian classification to predict the probability of occurrence of distinct water masses. The seismic section was processed to image in detail the fine scale structure of the water column using seismic oceanography. The processing sequence was developed to preserve, as much as possible, the relative seismic amplitudes of the data and enhance the shallow structure of the water column by effectively suppressing the direct arrival. The migrated seismic oceanography section shows an eddy at the expected Mediterranean Outflow Water depths, steeply dipping reflectors, which indicate the possible presence of frontal activity or secondary dipping eddy structures, and strong horizontal reflections between intermediate water masses suggestive of double diffuse processes. We then developed and applied an iterative geostatistical seismic oceanography inversion methodology to predict the spatial distribution of temperature and salinity. Due to the lack of contemporaneous direct measurements of temperature and salinity we used a global ocean model as spatial constraint during the inversion and nearby contemporaneous ARGO data to infer the expected statistical properties of both model parameters. After the inversion, Bayesian classification was applied to all temperature and salinity models inverted during the last iteration to predict the spatial distribution of three distinct water masses. A preliminary interpretation of these probabilistic models agrees with the expected ocean dynamics of the region.


INTRODUCTION
Fine-scale ocean processes happening on ranges from a few meters to a few kilometers have a profound impact on turbulent dynamics, on the ocean energy budget, on primary production and ecosystems, on gas and tracer exchange, and ultimately on the global ocean circulation and climate (e.g., Wunsch and Ferrari, 2004;Mahadevan, 2016). Yet, ocean measurements at high resolutions are limited to fixed point probes or profiling devices. Therefore, quasi-synoptic measurements of simultaneously vertical and lateral high resolutions require detailed planning and the combination of several probing devices, which are not globally available to sample the ocean (Pascual et al., 2017).
In SO, multichannel seismic reflection data is processed to boost the amplitudes corresponding to seismic reflections occurring at interfaces between water layers of different temperature and salinity, where water temperature is the most important property and contributes on average 80% to the reflection coefficient (Ruddick et al., 2009;Sallarès et al., 2009). A challenge for SO data processing is the ability to image reflections near the sea surface, as traditional seismic processing sequences fail to successfully mitigate the effect of the direct wave traveling from source to receivers on the very weak reflections in the near-surface water layer (e.g., Ruddick et al., 2009;Pinheiro et al., 2010). Several processing workflows have been proposed to tackle this limitation. Huang et al. (2012) used an adaptive subtraction scheme. Hardy et al. (2007) and Jones et al. (2008) combined linear moveout with dip filtering. Ristow et al. (2017) used a combination of a linear Radon transformation with adaptive subtraction. Often, another source of coherent noise originates from the echo of the previous shot, but this noise is not addressed herein. The objective of the seismic oceanography processing shown herein is twofold: (1) to effectively attenuate the direct arrival effect using a combination of linear moveout, horizontal median filtering and adaptive subtraction; and (2) to preserve the relative seismic amplitudes of the seismic data. The attenuation of the direct wave arrival ensures a good image of the water structure in the first few hundred of meters, whereas preservation of original amplitudes is required to invert the seismic data for the ocean's physical properties.
Seismic oceanography data represent an indirect measurement of the physical properties of the water column such as temperature and salinity. In fact, seismic reflections present in SO data originate from the interfaces between water masses with distinct properties. From an oceanographic perspective, having the ability to predict the spatial distribution of such properties from SO data would provide insights about oceanographic processes not detected by conventional oceanographic sampling techniques. The spatial prediction of such properties from seismic oceanography data is an inverse problem. Mathematically, seismic oceanography inversion can be expressed as: where F is the forward operator through which the recorded seismic amplitudes (d obs ), with d obs R d , are obtained from an ocean's model, m R m , and e represents the error term associated with the observations and modeling uncertainties errors present in the seismic oceanography data. In the seismic oceanography case, the ocean's acoustic properties (P-wave propagation velocity and density), m, can be computed from temperature and salinity using the thermodynamic equation of seawater (IOC, SCOR, and IAPSO, 2010). Seismic inversion methods can be broadly divided in two different classes: deterministic or statistical (Bosch et al., 2010). Deterministic approaches are based on regression models of optimization algorithms providing a single best-fit solution. In deterministic seismic inversion the uncertainty assessment is limited and defined as a linearization around the best-fit inverse solution, which is normally retrieved by least squares, and in this sense, the uncertainty is strictly represented by a local multivariate Gaussian (Tarantola, 2005). In statistical seismic inversion, the solution is expressed as a probability density function in the model parameters space. Therefore, these inversion methods provide a set of alternative models as solution and allow the assessment of the uncertainty associated with the inverted models (Tarantola, 2005). Assessing the uncertainty about the predictions in seismic inversion is critical in any modeling procedure. The uncertainty represents the lack of knowledge about the system under investigation, measurement errors and physical approximation during the data processing (Tarantola, 2005). Also, accounting for uncertainty, and therefore risk, leads to better-informed decisions.
Iterative geostatistical seismic inversion methods allow predicting models at higher resolution than the observed data due to their ability to incorporate high-resolution information provided by existing direct observations (e.g., Azevedo and Soares, 2017). This is particularly of interest in oceanographic studies as it might open a window to a reality not yet known. Conventional ocean models built exclusively from the interpolation of sparse direct measurements of ocean properties, such as CTD and/or XBT, are always a smooth representation of the true ocean variability and are unable to describe complex features in detail due to the large distances between observations. The inverted models of the ocean properties retrieved from SO data are much richer from a spatial perspective as the SO data constrains the model predictions far from the location of the direct observations (e.g., Dagnino et al., 2016).
In a preliminary work, geostatistical inversion has been successfully applied to predict the spatial distribution of ocean temperature and salinity from seismic oceanography data (Azevedo et al., 2018). These authors showed how geostatistical SO inversion could retrieve a set of high-resolution temperature and salinity models that generate synthetic SO data consistent with observations. Each model represents an alternative scenario that fits equally well the observed data SO data. However, in this approach the model perturbation technique (i.e., geostatistical simulation) (Deutsch and Journel, 1992) requires the existence of contemporaneous and collocated direct measurements of temperature and salinity (e.g., CTD/XBT casts) along the SO section to be inverted. The simultaneous acquisition of CTD/XBT data is difficult due to both operational challenges and costs, and represents a major drawback in the practical use of this technique.
The main objective of this work is to propose an alternative geostatistical SO inversion method to overcome the need for contemporaneous and collocated direct observations of temperature and salinity, which might open the door to the generalization of this type of inversion method in SO studies. Low-resolution models of temperature and salinity are extracted from large-scale ocean simulations and integrated as part of the objective function within the geostatistical SO inversion method. The low-resolution models represent a background model with the expected spatial trend of temperature and salinity (Pereira et al., 2019). In practice the background models act as spatial constraints in the inversion procedure. These models are not included as part of the model parameter space to avoid limiting the exploration of the model parameter space and therefore in a limited uncertainty assessment. The marginal and joint distributions of temperature and salinity, necessary to perform the geostatistical simulations, are borrowed from quasicontemporaneous and quasi-collocated ARGO floats profiles (Argo, 2000) which were acquired approximately simultaneous during the acquisition of the MCS profile.
As part of the proposed workflow, and to interpret the inverted models obtained in the last iteration of the inversion procedure, we classified the set of inverted models generated during the last iteration of the geostatistical SO inversion into distinct water masses using Bayesian classification (Avseth et al., 2005). From the classified models we computed the probability of occurrence of each water mass. The most likely depths of the different water masses agree with the expected values for the area of interest as proposed by other authors (Comas-Rodríguez et al., 2011).
The following section presents a detailed description of the proposed geostatistical SO inversion methodology and the Bayesian classifier used to predict the spatial distribution of the expected water masses. The proposed inversion method is then applied to a MCS profile acquired over the Madeira abyssal plain (MAP) that was specifically reprocessed for this purpose. The tailored seismic processing workflow is described in the subsequent section. We then show a preliminary and global interpretation of the oceanographic insights provided by the inverted models. The workflow proposed in this work can be applied to other locations worldwide where no contemporaneous direct measurements of salinity and temperatures and global ocean models exist in the vicinity of the seismic profiles.

GEOSTATISTICAL SEISMIC OCEANOGRAPHY INVERSION
The proposed iterative geostatistical SO inversion method inverts SO data directly for temperature and salinity when no contemporaneous direct observations of the ocean are available. It can be considered an extension of the method introduced by Azevedo et al. (2018). It allows the simultaneous integration of high-resolution direct measurements of temperature and salinity, such as CTD and XBT data, and background models from climatology, data products or numerical ocean simulations. A relevant aspect of the proposed method is that the background models do not constrain the model generation (e.g., used for example as a local mean during the geostatistical simulation), but are included as part of a two-term objective function where they act as a spatial regularizer (Pereira et al., 2019).
The proposed iterative geostatistical SO inversion method can be divided into four main steps: (i) generation of temperature and salinity background models; (ii) generation of highresolution temperature and salinity models using stochastic sequential simulation and co-simulations; (iii) multi-objective mismatch evaluation; (iv) stochastic update and generation of a new set of models.
In the application example shown herein, two-dimensional large-scale temperature and salinity sections, describing the expected background spatial distribution, were retrieved from global numerical simulations of ocean dynamics provided by the Copernicus Marine Environment Monitoring Service  (CMEMS, 2016). These sections were extracted for the same geographical location, acquisition date and time of the available SO section. From a geophysical inversion point of view, these simulated profiles can be thought as low-frequency models in conventional seismic inversion methodologies (e.g., Sams and Carter, 2017;Pereira et al., 2019). Contemporaneous observed high-resolution temperature and salinity vertical profiles acquired close to the seismic oceanography profile were used to infer the marginal and joint distributions of both ocean properties. This information was used during the model generation and perturbation and not as spatial constraining data. These data were used as target marginal and joint distributions to be reproduced by the geostatistical simulation algorithms. In this inversion method, water temperature models are generated with direct sequential simulation (DSS; Soares, 2001) and salinity models were cosimulated with co-DSS with joint probability distributions (Horta and Soares, 2010). The sequential generation of models ensures that the observed relationship between temperature and salinity is reproduced in any given pair of models generated during the iterative procedure. This is essential to guarantee the plausibility of the predicted location and extent of water masses and for classification of distinct water masses as shown below.
The pairs of temperature and salinity models are used to compute water density and P-wave propagation velocity using the international thermodynamic equation of seawater (IOC, SCOR, and IAPSO, 2010). Then, normal incidence reflection coefficients are computed and convolved with a representative wavelet. In cases where contemporaneous and collocated vertical profiles of salinity and temperature are available, the wavelet can be  estimated by comparing synthetic seismic traces with the real data (i.e., as in conventional seismic-to-well tie). For this work, however, the wavelet was extracted from the processed SO data by averaging the primary seafloor reflection of several traces selected from a region of relatively flat bathymetry, following the approach used by Warner (1990). The resulting synthetic traces are then compared against the corresponding real traces following: where x s and y s are the real and synthetic seismic traces, respectively, with N seismic samples. The similarity, S, like the Pearson's correlation coefficient is bounded between −1 and 1, but ensures a simultaneous match of the synthetic seismic on both waveform and amplitude values of the recorded SO data. The plausibility of the inverted models depends on the reproduction of both properties of the observed data. The deviations (dev) of each single realization of temperature and salinity from the background ocean models are computed following: where m sim is a realization of temperature or salinity and m background is the corresponding background model. Finally, a two-term objective function (OF) is computed combining Eqs.
(2) and (3): where w 1 and w 2 are user defined weights that sum to 1 and control the influence of each term depending on the quality of the existing SO data and the reliability of the background model; dev T and dev S are the deviations of the simulated parameter from the background model for temperature and salinity, respectively. If the quality of the SO data is low, then w 1 should decrease. Similarly, if the reliability of the background temperature and salinity is poor, w 2 should be reduced. In the application example shown herein w 1 and w 2 were set by trial-and-error. However, these can be optimized under an optimization framework (e.g., Gennert and Yuille, 1988;Mead, 2008;Marler and Arora, 2010). Notice that OF is also bounded between −1 and 1 so it can be used as a proxy of a collocated correlation coefficient in the geostatistical co-simulation of a new set of temperature and salinity models in the subsequent iteration of the inversion (e.g., Soares, 2001). For a given iteration (j), the pairs of temperature and salinity traces (i.e., vertical columns of grid samples) that ensure the maximum OF values are stored in auxiliary temperature and salinity models along with the corresponding OF values. These models are used as secondary variables in the co-simulation of a new set of models in the subsequent iteration. In practice, regions of the seismic profile with low OF will exhibit a large variability of simulated values within the ensemble of simulated models, while region of high OF will produce similar models in the next iteration.
The proposed iterative geostatistical SO inversion method is summarized in the following sequence of steps (Figure 1): (i) Create temperature and salinity background models for the entire inversion grid. The background models might be vertical sections extracted from numerical ocean simulations collocated with the existing SO data for the same acquisition time; (ii) Stochastic sequential simulation (direct sequential simulation; Soares, 2001) of Ns temperature models for the entire inversion grid. Direct temperature measurements located nearby the location of the SO profile are used to infer the conditioning distribution; (iii) Stochastic sequential co-simulation (direct sequential cosimulation with joint probability distributions; Horta and Soares, 2010) of Ns salinity models for the entire inversion grid. Direct salinity measurements located nearby the region of interest are used to infer the conditioning distribution. Each temperature model simulated in ii is used as an auxiliary variable to ensure the relationship between both ocean properties are reproduced in each Ns pairs of simulated models; (iv) For each pair of simulated ocean models, calculation of Ns reflection coefficient models using the International  (v) Calculate the objective function (Eq. 4) on a trace-by-trace basis; (vi) Select and store the temperature and salinity traces that generated the highest OF value in best local salinity and temperature models along with the corresponding OF; (vii) Use these local best OF, temperature and salinity models as secondary variables in the co-simulation of a new set of temperature and salinity models; (viii) Return to ii and iterate until the global correlation coefficient between real and synthetic SO data is above a certain threshold or a pre-defined number of iterations is reached.
Temperature and salinity models simulated and co-simulated during the last iteration generate highly correlated synthetic SO data with the observed data. These models were classified into distinct water masses. Bayesian classification (e.g., Avseth et al., 2005;Grana et al., 2017) was trained based on the existing direct measurements of temperature and salinity profiles of spatially located nearby Argo floats taken during the acquisition of the seismic oceanography data. N w different types of water were identified including the Central Atlantic Water, the Mediterranean Outflow and the Subarctic Intermediate Waters, as described in Comas-Rodríguez et al. (2011). The statistical properties (i.e., mean, covariance and proportions) were inferred from the training data and used to compute the prior and likelihood function for the Bayesian classification according to Bayes' rule: P k|d = P d|k P(k) P(d) = P d|k P(k) N w k = 1 P d|k P(k) , k = 1, . . . , N w (5) where, d is the vector of the ocean properties used for the classification, the simulated pairs of models, and k is the number of water masses. In Eq. (5), P d|k is the likelihood function, P(k) is the prior model and P(d) is a normalization constant. The set of Ns models classified in k water masses was then used to compute the probability of occurrence of each water mass.
Finally, the pointwise average models of temperature and salinity inverted during the last iteration of the inversion and the water probability sections were used to perform a simple and preliminary interpretation of the oceanographic features observed in the data. Nevertheless, this is not the main focus of this work.

Dataset Description
The proposed iterative geostatistical SO inversion method was applied to a seismic profile (WM-MAD01-003) acquired over the MAP with conventional MCS reflection methods (Figure 2). A summary of the main acquisition parameters is shown in Table 1. This seismic profile was acquired by the Portuguese Task Force for the Extension of the Continental Shelf (EMEPC in its Portuguese acronym) between June 6 and June 8, 2006 and is part of a larger seismic dataset located within the MAP. The typical thermohaline structure of the water column in this area is characterized by surface waters of subtropical type (warmer and saltier) over central waters of subpolar origins (Central Atlantic Waters) with lower temperature and salinity. Below the Central Atlantic Waters between about 500-1,500 m the water becomes saltier and warmer due to the presence of Mediterranean Outflow Water (MOW). Deeper in the water column, at the lower intermediate levels, temperature and salinity decrease with the presence of subpolar type intermediate waters (e.g., Segade et al., 2015). This vertical structure of water masses is unique as far as a considerable thermohaline structure is enclosed in the upper 2,000 m of the water column and a clear structure in the SO data was expected. Besides, this area is characterized by recurrent eddy activity associated with the energetic and unstable Azores Current jet on the upper ocean (Barbosa Aguiar et al., 2011), and with the main path of propagation of the Mediterranean Water Eddies (Richardson et al., 2000;Barbosa Aguiar et al., 2013), which carry very distinct salty and warm water anomalies within its core at intermediate levels.
The MCS profile was processed to image the fine-structure in the water column. Particular attention was given to mitigate the effect of the direct arrival and enhance shallow reflections while preserving true amplitudes by applying processing parameters that minimize amplitude and phase distortion. Different processing sequences would produce different results and impact the models predicted with the SO inversion. The detailed description of the processing sequence is presented in the following sub-section.
As collocated and contemporaneous direct measurements of temperature and salinity were not available, we used three vertical profiles of temperature and salinity (from two ARGO floats- ARGO,4,660 and 44,909), profiling close to the acquisition period (in May 21 and May 31, 2006) and located in the surroundings of the SO section during the data acquisition (Figure 3). The ARGO profiles were used to infer the marginal (Figure 4) and joint distributions (Figure 5) of both ocean properties and were used as conditioning distributions for the stochastic sequential simulation and cosimulation of temperature and salinity models. We assume that the statistical properties of temperature and salinity measured by these floats hold for the location of the SO profile. According to the temperature and salinity measurements, three distinct water masses could be inferred: central Atlantic water; MOW; and Subarctic Intermediate Water, as described in Comas-Rodríguez et al. (2011). Low-frequency temperature and salinity background models were built using a global ocean dynamics   model (CMEMS, 2016) for the dates of seismic acquisition (Figure 6). These two-dimensional sections are collocated with the MCS data but are smooth and low resolution.
While the background temperature model shows a vertical trend of high temperature at shallower water depths and low temperature at deeper depths, the background salinity model follows the description above and considers the effect of the MOW (i.e., a saltier water layer) around 1,125 m of water depth.

Seismic Processing of Line WM-MAD01-003
The processing sequence (Figure 7) included data resampling and recording length reduction to comprise data exclusively above the seafloor. Bad traces, both due to poor signal-to-noise ratio or bad readings, were edited and those unrecoverable were removed from the dataset. The direct arrival was tackled by applying a horizontal median filter with adaptive amplitude subtraction, similar to ocean bottom seismometers data processing (Duncan and Beresford, 1995). This kind of amplitude subtraction aims at minimizing the effects on the resulting amplitudes. This process was performed sequentially in a four-step approach (Figure 8): first, field records were flattened using a linear moveout correction with a constant velocity of 1,500 m/s ( Figure 8A); the flattened records were doubled to avoid edge effects when applying a median filter ( Figure 8B); a horizontal median filter was applied to those records to preserve horizontal coherent reflection ( Figure 8C). Finally, the resulting record was subtracted to the doubled flatten record (Figure 8E) to eliminate the effect of the direct arrivals and keep reflections in the records (Figure 8F). The filtered gather, after applying a normal moveout (NMO) correction with a constant velocity of 1,500 m/s, is shown in Figure 8G where the reflections within the water column around 500 ms are enhanced.
The second part of the processing sequence comprises bandpass filtering between 10 and 80 Hz, and surgical mutes to remove FIGURE 15 | Probability of occurrence of: (A) Central Atlantic Water; (B) Mediterranean Outflow Water; (C) Subarctic Intermediate Water. Probability models calculated from the classification of the thirty-two temperature and salinity models generated during the last iteration of the inversion using the training data shown in Figure 5.
bursts of energy at the smallest offsets (Pinheiro et al., 2010). True amplitude recovery was applied to compensate spherical divergence. A detailed velocity analysis was performed along the MCS profile. The resulting velocity field was used in the normal moveout correction of the records. After CMP sorting, the CMP gathers were stacked considering all offsets. Finally, a constant velocity (1,500 m/s) phase-shift migration (Stolt, 1978) and time-to-depth conversion were carried out using the same constant velocity model.
Due to computational constraints, the MCS profile was processed in swaths of 400 field records, with an overlap of 50 records. This computational limitation results in the vertical stripping observed in the final time migrated section, which also affect the inverted temperature and salinity models (Figure 9). All sections are plotted in the depth domain to ease interpretation by assuming an average P-wave velocity of 1,500 m/s. However, the inversion of the SO section was performed in the time domain (i.e., prior to depth conversion).

Preliminary Interpretation of Line WM-MAD01-003
A preliminary interpretation of the SO profile allows its vertical division in two layers: the top one, down to approximately 1,875 m, comprises bright and coherent reflections with different seismic signatures and dips; the bottom one, below this depth, which is relatively reflection-free maybe due to the relatively homogeneous North Atlantic Deep Water. For this reason, only the top layer of the SO profile was considered for the geostatistical inversion (i.e., from 0 to 2.5 s in TWT).
The top layer shows the presence of an eddy at depths where Meddies are typically expected 600-1,500 m ( Figure 10A). A smaller lenticular feature is also observed at shallow depths (around 750 m) that could be associated with the transition to Mediterranean waters or a secondary feature associated with that eddy (Figure 10B). Oblique reflections between 750 and 2,250 m might be related to oceanic fronts within the MOW or an inclined eddy of smaller dimensions as imaged in Tang et al. (2020). A possible double diffusion phenomenon can be detected by the continuous, parallel and bright reflections at approximately 1,500 ms for almost all the profile (Figure 10C). Double diffusion generates staircase thermohaline structure. The seismic signature of this oceanic structure has already been observed and investigated in the area near the Lesser Antilles in the Caribbean Sea (e.g., Fer et al., 2010) and in the Gulf of Cadiz associated to Mediterranean Outflow Water eddies (e.g., Biescas et al., 2008).

INVERSION RESULTS
The geostatistical SO inversion ran with six iterations, where at each iteration thirty-two pairs of temperature and salinity models were generated using direct sequential simulation (Soares, 2001) and direct sequential co-simulation with joint probability distributions (Horta and Soares, 2010). The stochastic sequential simulation and co-simulation were conditioned to the histograms and bi-histograms from Figures 4, 5, respectively. Since the ARGO floats were not located along the seismic profile, no spatial conditioning was considered. Vertical variograms for each property were modeled from both ARGO floats while the horizontal variogram was modeled directly on the seismic amplitudes. This is a conventional approach in seismic reservoir characterization and often results in overestimation of the horizontal ranges of the variogram model (Azevedo and Soares, 2017). The objective function (Eq. 4) used to drive the inversion is based on the mismatch between synthetic and real seismic traces used for the inversion and the deviation of each realization of temperature and salinity from the background models shown in Figure 6.
The synthetic SO data computed from the pointwise average temperature and salinity models generated during the last iteration is shown in Figure 11. These data reproduce the location of the main oceanographic features as interpreted from the observed data as well as their amplitude content. As expected, the synthetic SO data calculated from these models is less noisy than the observed seismic (e.g., Avseth et al., 2005) and, consequently, increases the spatial continuity of the seismic reflections at the bottom part of the section (∼1,500 m). Iterative geostatistical seismic inversion methods are known for the ability to remain unmatched in noisy areas of the observed data (Azevedo and Soares, 2017). This effect is illustrated by the lack of vertical artifacts in the synthetic data. To illustrate the local convergence of the synthetic data, Figure 11C shows the trace-by-trace S between true and inverted SO data.
The inverted temperature and salinity models (Figure 12) capture the oceanic features of interest at finer scale when compared with the observed SO data. This effect is related to the use of geostatistical simulation as model perturbation technique, the geostatistical simulation fills-in the frequency band related to high-frequencies (Azevedo and Soares, 2017). Modeling oceanographic features at these scales is not possible with either deterministic seismic inversion methods or conventional interpolation techniques of direct observations of the ocean properties as represented by CTDs or XBTs. This is one of the main benefits of using geostatistical seismic inversion methods. These models show the filamentation structures around the eddie's core and in particular the warm intrusions around the homogenous nucleus. The use of background temperature and salinity models (Figure 6) allows reproducing the expected largescale vertical distribution of both properties as interpreted from the global ocean models.
Additionally, the benefit of using geostatistical inversion methods is related to the ability to assess the uncertainty associated with the model predictions. Figure 13 shows the pointwise variance of temperature and salinity computed from the ensemble of models generated for each property during the last iteration of the inversion procedure. It is interesting to discuss the spatial distribution pattern of these models. As temperature is the main contributor for the existence of reflection coefficients (Ruddick et al., 2009;Sallarès et al., 2009) the spatial uncertainty (i.e., the variance) is smaller in regions where the observed SO data exhibits coherent seismic reflections (i.e., above approximately 375 m and below 1,500 m). On the other hand, the region of lower variance for salinity, between the 750 and 1,125 m, matches the depths associated with the saltier layer as observed in the background model (Figure 6). The reason for this phenomenon still needs to be further investigated but might be related to: (i) the differences in signal-to-noise ratio in different parts of the seismic section; (ii) the local influence of the background models.
When contemporaneous direct measurements of temperature or salinity along the SO profile are available, one could assess the local performance of the inversion by retaining one observation out of the conditioning data and comparing the inverted traces with the observed data (i.e., a blind-well test in subsurface modeling). In this application example we compare the depth trend of the inverted two-dimensional sections of temperature and salinity of all the realizations generated during the last iteration of the inversion with the vertical one-dimensional profiles acquired by the ARGO floats (Figure 14). Note that the ARGO floats are not used to spatially constrain the inversion. In this simple exercise we aim at evaluating if the vertical trend of both properties is reproduced. As expected the direct measurements are not exactly reproduced but we consider the reproduction of the main trends to be a positive result.To illustrate the potential of geostatistical seismic inversion methods we used the ensemble of pairs of temperature and salinity models generated during the last iteration to generate two-dimensional sections of probability of occurrence of different water masses. First each pair of models was classified into three distinct water masses using Bayesian classification.
The a priori probabilities were inferred from the ARGO floats profiles, which were used as a training dataset (Figure 5). After classification of each pair, the ensemble of 32 models was used to compute the probability of occurrence of each water type (Figure 15). The resulting probability sections agree with the overall knowledge of this oceanic basin and with previous results obtained using exclusively large-scale oceanographic observations (Comas-Rodríguez et al., 2011). It is relevant to highlight that the order relationship of the different water masses (i.e., the Central Atlantic Water above the MOW above of the Subarctic Intermediate Water) is reproduced in the probability models, but it is not imposed by any other information rather than the SO data and the background models. As expected the regions of higher uncertainty are located at the boundaries of each water mass.

CONCLUSION
This paper presents the first seismic oceanography images of the Madeira abyssal plain region. The work focuses on two main aspects: (i) the introduction of a simple but efficient way to mitigate the effect of the direct arrival in the data; and (ii) the development of a geostatistical SO inversion that can be used when contemporaneous and collocated direct measurements are not available. The interpretation of the inverted models from an oceanographic perspective is limited, as it would benefit from the processing and inversion of the other two adjacent SO sections.
The processing sequence applied to these data was able to effectively attenuate the effect of the direct wave, revealing reflections in the first few hundred meters below the sea surface (Figure 10). The time-migrated section clearly shows fine-scale structure down to 2,000 m, below this depth the SO data shows no reflection. The upper part of the section exhibits a series of interesting oceanographic features that might be interpreted as eddies associated with the MOW and as double diffusive phenomena. However, these features need to be further explored to provide insights about the complex dynamics of the study area (Figure 10). The detailed interpretation of these oceanographic features will be performed after inverting the neighbor SO profiles existing in the region.
The processed time-migrated section was inverted using the proposed geostatistical SO inversion. We show how this inversion technique can be applied when no direct and contemporaneous observations of the ocean are available. We leverage common models of large-scale ocean dynamics and existing vertical profiles of the ocean properties measured by ARGO floats. This method was proved to be a useful tool to characterize submesoscale oceanic features and has demonstrated a potential to invert for temperature and salinity. From the inverted models we also propose Bayesian classification of water masses. The probability of occurrence of the different water masses along the profile clearly agrees with the known vertical distribution (Figure 15).
As a final remark, it is worthwhile to highlight that the ensemble of inverted temperature and salinity fields generated at the last iteration still exhibits large variability (Figure 14). While this uncertainty is related to the lack of knowledge about the system under investigation (i.e., the ocean properties), it may be reduced and better predictions can be obtained by including additional constraints from other oceanographic variables (e.g., density, water column stability) during the inversion procedure. This approach would increase the reliability of the inverted models.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
LA conceptualized and implemented 50% of the geostatistical SO inversion and wrote of the manuscript. LM conceptualized and supervised the SO processing data and wrote part of the manuscript. FT implemented 50% of the geostatistical SO inversion, ran the inversion on the application example, and wrote the manuscript. RT applied the SO processing flow. ÁP contributed for the preliminary interpretation of the results from an oceanographic point of view and wrote the manuscript. All authors contributed to the article and approved the submitted version.