Stochastic Heterogeneity Mapping around a Mediterranean salt lens

. We present the ﬁrst application of Stochastic Heterogeneity Mapping based on the band-limited von K ´ arm ´ an function to a seismic reﬂection stack of a Mediterranean water eddy (meddy), a large salt lens of Mediterranean water. This process extracts two stochastic parameters directly from the reﬂectivity ﬁeld of the seismic data: the Hurst number, which ranges from 0 to 1, and the correlation length (scale length). Lower Hurst numbers represent a richer range of high wavenumbers and correspond to a broader range of heterogeneity in reﬂection events. The Hurst number estimate for the top of the meddy (0.39) compares well with recent theoretical work, which required values between 0.25 and 0.5 to model internal wave surfaces in open ocean conditions based on simulating a Garrett-Munk spectrum (GM76) slope of − 2. The scale lengths obtained do not ﬁt as well to seismic reﬂection events as those used in other studies to model internal waves. We suggest two explanations for this discrep-ancy: (1) due to the fact that the stochastic parameters are derived from the reﬂectivity ﬁeld rather than the impedance ﬁeld the estimated scale lengths may be underestimated, as has been reported; and (2) because the meddy seismic image is a two-dimensional slice of a complex and dynamic three-dimensional ﬂow. Nonetheless, parameters, in the Garrett-Munk spectrum (horizontal wavenumber spectrum), can provide an estimate of different internal wave scales from seismic data alone. We hence introduce Stochastic Heterogeneity Mapping as a novel tool in physical oceanography.


Introduction
Mediterranean Water eddies, or "meddies", are large, warm, isolated lenses of highly saline Mediterranean Water that are found in the North Atlantic ocean. Mediterranean Water flows through the Strait of Gibraltar as an undercurrent (Bower et al., 2002), cascades down the continental shelf, while entraining less dense North Atlantic Central Water (Bower et al., 1997) and settles at depths between 500 and 1500 m (Richardson et al., 2000). The undercurrent then rounds the corner of the Iberian peninsula at Cape St. Vincent, directed north by the Coriolis force. It is here that meddies form, spinning off the main undercurrent, translating westward and rotating anti-cyclonically (Serra et al., 2002). Meddies were first reported in the western North Atlantic ocean by McDowell and Rossby (1978). Since that time they have been found to be a common feature in the North Atlantic ocean, (Richardson et al., 2000). Many different aspects of meddies are currently being researched to understand their properties as well as their influence on large-scale mixing and climate (e.g. Bashmachnikov et al., 2009).
Meddies have traditionally been studied using established oceanographic techniques, such as CTD (Conductivity-Temperature-Depth) probes to measure salinity and temperature (Ruddick, 1992) and acoustically tracked SOFAR floats Published by Copernicus Publications on behalf of the European Geosciences Union. (Armi et al., 1989). Recently, multi-channel seismic (MCS) reflection profiling has been employed as a thermohaline imaging tool starting with the work of Holbrook et al. (2003) (e.g. Nandi et al., 2004;Nakamura et al., 2006;. Holbrook and Fer (2005) used this technique to study internal waves in the Norwegian Sea. They found that horizontal wave number spectra derived from digitizing seismic reflector horizons in the open ocean compared favorably with the Garrett-Munk spectrum (Garrett and Munk, 1975), which describes oceanic internal wave displacements. Biescas et al. (2008) performed the first detailed MCS analysis of a meddy. They found distinct seismic reflectivity differences in the upper and lower bounds of the meddy that are consistent with differences seen in historical temperature and salinity data. The upper boundary of the meddy was characterized by a few high-amplitude, laterally continuous reflections, whereas the lower boundary exhibited more numerous, shorter, lower amplitude reflectors.
To further understand meddy processes, we apply Stochastic Heterogeneity Mapping to the study of a meddy in the Gulf of Cadiz (Fig. 1). This method of statistically analyzing the reflection field of seismic data has been used extensively to study complex acoustic impedance variability in the solid earth (e.g. Goff et al., 1994;Holliger et al., 1994;Hurich and Kocurko, 2000;Carpentier and Roy-Chowdhury, 2007) even in areas where there is a predominance of diffuse reflectivity (as opposed to specular reflectivity) observed. This research represents the first application of Stochastic Heterogeneity Mapping to the ocean.
Stochastic Heterogeneity Mapping is based on the premise that the seismic reflection wave field contains information on the spatial properties of the reflecting bodies and could thereby be used to extract quantitative information about thermohaline finestructure. By extracting these parameters from a stacked seismic image of a meddy we can estimate the range of scales of its reflectivity patterns. In this way, we provide an estimate of the characteristic scales of internal waves of physical oceanographic processes directly from seismic data on zones of particular interest to oceanographers: the top and bottom (Biescas et al., 2008) and the sides of a meddy (Armi et al., 1989;Ruddick, 1992). Since Stochastic Heterogeneity Mapping operates on the reflectivity field of a processed seismic reflection profile, we focus here on that methodology and its interpretive implications, not on MCS data acquisition or processing. For an introductory treatise of MCS methodology as it relates to physical oceanography, we refer the reader to Ruddick et al. (2009). See Yilmaz (1987) for a more complete discourse on seismic data processing.

The stochastic model
The 1-D von Kármán model is described by three parameters: the correlation length (a), which is the upper limit for the scale invariance in heterogeneity (Carpentier, 2007), the Hurst number (ν), a measure of surface roughness or equivalently, the richness of the range of scales in the power law distribution (having values between 0 and 1), and statistical variance. We apply unit variance to standardize the distribution. For scale sizes smaller than the correlation length, the von Kármán model describes a power law (fractal) process, where ν represents its exponent. The parameter, ν relates to the fractal dimension (D) by where E is the Euclidean dimension. For scales longer than the correlation length, the von Kármán model represents a process that is uncorrelated, such as white noise (Hurich and Kocurko, 2000). The structure of the impedance field, and hence its autocorrelation, are in accordance with the defined power spectrum. This spectrum could take on a variety of forms. However, we choose a von Kármán stochastic distribution because it is capable of describing a band-limited power law process and has been thoroughly tested for this algorithm, albeit for deep crustal studies (Carpentier and Roy-Chowdhury, 2007). However, the original work of Theodore von Kármán was to characterize random fluctuations in the velocity field of a turbulent medium (von Kármán, 1948), indicating that the method is also suited to the characterization of ocean fluid dynamics. We express the analytic radial 2-D von Kármán power spectrum as (Carpentier, 2007), where a x and a z are the correlation lengths in the lateral and vertical directions, respectively, ν is the Hurst number, k is the weighted radial wavenumber, k 2 x a 2 x + k 2 z a 2 z . In the space Ocean Sci., 6, 423-429, 2010 www.ocean-sci.net/6/423/2010/ domain, we express a 2-D autocorrelation function (Goff and Jordan, 1988), where G ν (r) = r ν K ν (r), K ν (r) being the second modified Bessel function of fractional order, and r the weighted radial autocorrelation lag, defined as x 2 /a 2 x + z 2 /a 2 z . G ν (0) is defined according to Carpentier (2007), where the gamma function, (ν)=(ν − 1)!. From here we shall refer to correlation lengths as scale lengths because, in band-limited fractal fields, these lengths give the threshold below which scaling is determined by a power law (Carpentier and Roy-Chowdhury, 2007) and due to the fact that the decorrelation length is the dominant visual scale of the fabric morphology. Furthermore, we restrict our analysis to estimating horizontal scale lengths because determination of the vertical scale length is complicated by the estimation of source wavelet characteristics. This is because the source impulse is not a perfect spike, but is instead contaminated by side lobe energy, especially in the near offset traces (Yilmaz, 1987, pp. 17-20).
We estimate a x and ν by performing a 2-D Fast Fourier Transform, from which we derive a 2-D power spectrum. We sum over the frequency direction and apply an inverse Fourier transform to obtain an autocorrelation function of NxM samples Carpentier and Roy-Chowdhury, 2007). We choose window sizes to be representative of the areas of interest. Constraints are placed on window sizes to cover a minimum of 10 times the horizontal scale length, and at least 2 cycles vertically. Next, the 1-D analytic von Kármán autocorrelation function (Eq. 4) is fitted to the calculated autocorrelation function using a model-space grid search and an L2 norm misfit, providing the estimate of a x and ν. Temporal and spatial band-limiting of the broadband von Kármán spectrum of the seismic data affects the accuracy of estimating ν, which, based on impedance contrast fields and Primary Reflectivity Sections can be overestimated by a factor of 2. Likewise, a x can be underestimated by a factor of 3-6 (Carpentier and Roy-Chowdhury, 2007). More recently, for complex scattered visco-elastic reflection data, Carpentier et al. (2009a) found underestimation factors for a x of between 6 and 10 and overestimation factors of up to 10 for ν. These errors are due to the fact that, although the impedance field is highly correlated to the reflectivity field, they are not equivalent. Nonetheless, they represent the current state of the art of Stochastic Heterogeneity Mapping.
The Stochastic Heterogeneity Mapping algorithm accounts for the fact that reflectors may have some apparent dip by performing dip searching so that the values are derived along the angle of maximum coherence. However, apparent dips are small: 3 • on average across the seismic section, with a standard deviation of 0.3 • .

Results
Under the premise that the statistical properties of the scattered wave field are highly correlated to the properties of the acoustic impedance field (Carpentier and Roy-Chowdhury, 2007), we analyzed a meddy in the seismic line GO-LR-05 acquired in the Gulf of Cadiz during the GO (Geophysical Oceanography) cruise of April and May, 2007 ( Fig. 1) for the distribution of the stochastic parameters Hurst number (Fig. 2) and correlation length (Fig. 3) around its margins. Four zones are described: The top of the meddy (A), the bottom of the meddy (B) and its sides (C and D). The stochastic parameters were extracted from the seismic section and overlaid using two different color schemes. The mapping procedure reveals the stochastic heterogeneity of the thermohalinerelated fabric observed in the reflectivity field. Median and average Hurst numbers and scale lengths and their statistical distribution are reported in Fig. 4.
Zone horizontal wavenumber spaces (Garrett and Munk, 1975). To simulate the GM76 slope they used Hurst numbers in a range of 0.25-0.50 with horizontal scale lengths of 5-10 km. Internal waves develop in the ocean when density stratification is disturbed, thereby driving turbulent processes and diapycnal mixing (Müller and Briscoe, 2000;Garrett, 2003). The recorded seismic signal is reflected from ocean density stratification whether it is strong specular reflectivity or weak diffuse reflectivity. Acoustic impedance boundaries, which are what give rise to the reflectivity as revealed by the multichannel seismic (MCS) method, are expressed as the product of in-situ density and sound speed. Density and sound speed, in turn are functions of the relative proportions of temperature and salinity of the reflecting interface . Ruddick (1992) reported that the upper boundaries of meddies are dominated by diffusive convection processes, whereas the lower boundaries are susceptible to salt fingering. Salt finger scales are several orders of magnitude smaller than that directly resolvable by MCS methods (typically, in the tens of centimeter range, Linden, 1973). However, we may expect to detect them indirectly through our estimation of the Hurst number. The reason for this is that the Hurst number is a measure of surface roughness, low Hurst numbers being representative of a richer range of high wavenumbers. As seen in Fig. 4, Hurst numbers calculated for the bottom of the meddy are found to be dominated by lower numbers, whereas for the top and sides of the meddy, a lack of these lower Hurst numbers is found. This seems to indicate that, although we can not hope to actually resolve the salt fingers themselves with seismic frequencies, the richer range of high wavenumbers may point to the presence of such finer structures at the bottom as opposed to at the meddy's top and sides.
The scale lengths estimated from the seismic section are significantly lower than those used to model internal waves by Vsemirnova et al. (2009), even for Zone A (ca. 1 km compared to their 5-10 km). This difference is possibly due to the factor of 3-6 underestimation that was reported by Carpentier and Roy-Chowdhury (2007). Alternatively, the lower estimated scale lengths of the meddy may be explained by the fact that it is a dynamic, three-dimensional structure with currents that may be moving obliquely to the two-dimensional acquisition path. This is supported by an observation by Klaeschen et al. (2009), who estimated reflector motion by using an in-situ sound speed model near the NE edge of Zones A and B. They found longer horizontal wavelengths on the NE side of Zone A, as opposed to the SW side, and conclude that this is an indication of a different movement between the respective sides. Direct LADCP (Lowered Acoustic Doppler Current Profiler) measurements during the seismic acquisition confirmed that different parts of the meddy have different distributions of velocities, some along the path of the acquisition, some oblique to it . That is, the meddy was not simply in solid-body rotation, but was stretching slightly in a SW direction at the time of acquisition. This motion could disturb the reflector undulations that we observe, as measured by the extracted horizontal scale lengths. Since we are making a two-dimensional observation in a three-dimensional domain, the stochastic values we obtained are measurements of the component of the meddy motion in the plane of the seismic profile, rather than in the direction of meddy motion. Thus, we can expect our measurements of scale length to be shorter than theoretical predictions based on two-dimensional geometries.
Due to the band limiting of the seismic data, the scales observed are most likely confined to the lower extreme of mesoscale (2-200 km) size features. It follows that the stochastic parameters seen here are a reflection of the effects of internal waves, as seen in the Garrett-Munk spectrum, rather than smaller turbulent, ie. Batchelor scales (Batchelor, 1959). The notably lower Hurst numbers seen at the bottom of the meddy (Zone B) may be partially a result of the fact that the frequency content of the source used in this survey (10-70 Hz) was too narrow-band to recover the smaller thermohaline staircases known to occur at the base of meddies (Ruddick, 1992). Thermohaline staircases are well-known structures in the ocean that range from tens to hundreds of meters thick and are found in regions where both temperature and salinity Ocean Sci., 6, 423-429, 2010 www.ocean-sci.net/6/423/2010/ 17 to bottom represent analyzed parts of the meddy, as outlined in analysis boxes (Figure 2).
Values of ν were grouped into families incremented by 0.05 between 0 and 1. Values of a x were grouped into families incremented by 500 m between 0 and 6500 m. Median and average ν and a x were then calculated. increase upward in a manner that promotes salt fingering processes (Schmitt, 1994). Given the low frequency, narrowband source, we are only able to theoretically recover structures on the order of 5-75 m (Widess, 1973). Moreover, in the ocean there is both heat and mass diffusion across interfaces that produce a gradual change over a finite boundary width that may extend for several tens of meters ). Hence, due to the thickness and gradient of a reflecting interface, a safer upper estimate for maximum resolvable thickness would be about double this estimate, between 10 and 150 m, the reflectivity being frequency dependent. As a result, it is plausible that the scales of some structures at the base of the studied meddy are smaller than those resolvable by this seismic source. Horizontal resolution is determined by the Fresnel zone width (Yilmaz, 1987, p. 470), itself a function of frequency, sound speed and depth.
For the upper and lower limits of the meddy, considering a dominant frequency of 50 Hz, the horizontal resolution thus ranges between 12 and 30 m. The seismic data presented herein are unmigrated. Given the near horizontal reflectivity of the seismic data (dips average at about 3 • with a standard deviation of 0.3 • ) migration does not appreciably change the position of reflectors. With this in mind and due to the fact that the spatial bandlimit that is imposed by the migration operator smoothes the data, thus removing the information in the first few lags of the autocorrelation function, we chose to analyze the unmigrated stacks. Intriguingly, in recent studies by Carpentier et al. (2009a, b) Hurst numbers were found to actually exceed 1 in both synthetic and real data in order to obtain von Kármán www.ocean-sci.net/6/423/2010/ Ocean Sci., 6, 423-429, 2010 best fits. These exceptionally high Hurst numbers were said to be a result of the lateral smoothness of the data, caused by first-Fresnel-zone averaging during migration. Although migration, in principal, should completely collapse the 2-D seismogram's first-Fresnel zone, it is complicated by factors such as the lack of a sufficient recording aperture (zone of seismic illumination), complex scattering effects, wavelet attenuation and wavefront healing, and the uncertainty inherent in choosing the correct migration velocity (sound speed). Therefore there will always be some residual first-Fresnelzone averaging. Furthermore, 2-D migration in a 2-D seismic profile will only collapse the first-Fresnel-zone and does not address out-of-plane reflections. This result supports the fact that the Hurst numbers we obtained from unmigrated data are not in the upper end (nor do they exceed) the range of theoretical Hurst numbers, and therefore may be closer to the true values than those extracted from migrated data, which are biased by many of the factors aforementioned. While we cannot definitively state that the low-end Hurst numbers are not a simply result of some remaining uncorrelated noise (thereby representing high-wavenumber artifacts, not actual sub-wavelength thermohaline finestructure) we can neither dismiss the possibility that the information contained in the first few autocorrelation lags actually has a positive effect on the von Kármán fitting process. Uncertainty in both Hurst number and scale length could be further improved by inverting in-situ sound speed functions for acoustic impedance, thereby improving the certainty of correlation to the reflectivity field, upon which the stochastic parameters are extracted. Improving this certainty would also be useful to Stochastic Heterogeneity Mapping studies of the solid earth where sufficiently sampled sound speed functions are in practice unachievable.

Conclusions
We partition our observation of Hurst number distributions for the meddy into three distinct zones: top, bottom and sides. Our calculations of Hurst number for the top of the meddy agree with recent theoretical work, which used values between 0.25 and 0.5 to model internal wave surfaces in open ocean conditions based on simulating a Garrett-Munk (GM76) slope spectrum of −2. The corresponding scale lengths (correlation lengths) mapped over the same reflectivity field, however, do not fit as well to specific seismic reflection events. We suggest two possible explanations for this discrepancy: (1) because the stochastic parameters are derived from the reflectivity field rather than the impedance field the estimated scale lengths may be underestimated, as has been reported; and (2) because the meddy seismic image is a two-dimensional slice of a complex and dynamic three-dimensional object, the estimated scale lengths are likely skewed to the direction of flow, where the theory is most applicable. Nevertheless, this work illus-trates the potential of Stochastic Heterogeneity Mapping as a tool within the growing discipline of Seismic Oceanography because it allows an estimate of lateral scale ranges of reflection events, and therefore actual physical oceanographic processes, such as internal waves. The multi-channel seismic (MCS) method records seismic signal that is reflected from ocean density stratification (i.e. boundaries of acoustic impedance). These boundaries give rise to the reflectivity field. Acoustic impedance boundaries are functions of the relative proportions of temperature and salinity of a given reflecting interface, whether strictly specular reflectivity or more diffuse reflectivity. So we may conclude that the estimations of the Hurst numbers and scale lengths are representative of actual thermohaline finestructure. We therefore introduce Stochastic Heterogeneity Mapping as a new tool of physical oceanography. To improve upon errors in estimation of Hurst number and scale length in-situ sound speed functions should be inverted for acoustic impedance. In this manner stricter constraints on the degree of correlation to the reflectivity field, upon which the stochastic parameters are extracted could be obtained. These constraints, in turn, could be applied to solid earth studies, where such a sound speed function is practically unattainable.