Geostatistical characterization of internal structure of mass-transport deposits from seismic reﬂection images and borehole logs

down-slope disaggregation of the mass-ﬂow during transport. The dominant scale lengths can be used as a proxy for strain history, which can improve understanding of post-failure dynamics and emplacement of subacqueous mass-movements, important for constraining the geohazard potential from future slope failure.


I N T RO D U C T I O N
Subacqueous mass-movements such as slides, slumps and debris flows are capable of rapidly mobilizing and transporting large volumes of sediment in marine and lacustrine slope environments. They represent a significant geohazard to seafloor infrastructure (Piper et al. 1999;Carter et al. 2014) and to shoreline populations from slide-induced tsunami (Assier-Rzadkieaicz et al. 2000;Satake 2012). Such events are preserved in the geological record as mass-transport deposits.
One focus of geohazard research is to characterize the internal structure of mass-transport deposits to better understand the failure dynamics and emplacement of subacqueous mass-movements. Outcrop studies of exhumed 'fossil' mass-transport deposits identify a wide variety of internal structural fabrics, often showing complex and intense deformation (Pini et al. 2012). Lucente & Pini (2003) document low-angle thrusting, recumbent folding and progressive down-flow deformation within mass-wasting deposits outcropping in the Marnoso-Arenacea Formation in central Italy. They also identify kinematic indicators such as asymmetric folding and imbricated duplexes. Ogata et al. (2014) report soft-sediment deformation structures, slumping and intact blocks of substrate, ripped up and incorporated into the flow during sliding in mass-transport deposits caused by the collapse of carbonate platforms. These different internal fabrics reflect differing modes of slope failure, sediment properties and flow dynamics. Internal structure informs strain history and can thus characterize different flow regimes, enabling division of mass-transport deposits into, for example, headscarp (extensional), translational and toe (compressive) domains. It can also constrain flow kinematics such as run-out distance and flow acceleration, which play a large role in governing the geohazard potential of an event.
Acoustic reflection techniques are currently the only geophysical methods able to image submarine mass-transport deposits in situ. In recent years it has become increasingly common to study mass-transport deposits using seismic reflection imaging (Martinez et al. 2005;Moscardelli & Wood 2008;Berndt et al. 2012;Sun et al. 2017). Scientific drilling and coring is also commonly performed to estimate geotechnical and petrophysical parameters, such as undrained shear strength and excess pore pressure (Camerlenghi et al. 2007;Sawyer et al. 2009;Strasser et al. 2011;Dugan 2012). Sediment cores can give a high resolution 3-D reconstruction of strain fabric within mass-transport deposits (e.g. Strasser et al. 2011), but only for centimetre-scale structure at single point locations. Bull et al. (2009) catalogue a variety of internal structures seen in mass-transport deposits and mass-transport complexes from 3-D seismic reflection data. Jackson (2011) documents internally coherent rafted megablocks emplaced within more chaotic sediments in a mass-transport deposit from a seismic survey in the Santos Basin, offshore Brazil. Steventon et al. (2019) estimate the overall strain distribution within a mass-transport complex offshore Uruguay from 3-D seismic data by identifying individual seismic reflectors and using a backstripping approach.
Evidently conventional seismic interpretation techniques (horizon tracking) can be used to interpret internal structure of masstransport deposits. But this is only possible when (i) the deposit is well-imaged; (ii) there is sufficient internal reflectivity to generate seismic reflections and (iii) the scale of internal structure is above the seismic resolution. Very often, however, internal reflectors can appear chaotic, disordered or low-amplitude (e.g. Diviacco et al. 2006;Moscardelli & Wood 2008;Badhani et al. 2019). In fact, identifying apparently chaotic or transparent seismic facies is a common way to discriminate failed from unfailed sediments. Many studies use seismic attributes which are sensitive to discontinuous reflectors, for example the coherence attribute or the chaos attribute (Chopra & Marfurt 2016) to identify mass-transport deposits (Martinez et al. 2005;Gafeira et al. 2010).
This common lack of laterally continuous reflectors within masstransport deposits makes conventional interpretation of internal structure difficult. Instead, this study aims to characterize internal structure from seismic reflection data using a geostatistical approach. This avoids the subjectivity inherent to approaches such as horizon tracking and is applicable even when the chaotic nature of a deposit changes laterally, such as a progressive down-slope loss of horizon continuity (e.g. Badhani et al. 2019). The goal is to go beyond using non-dimensional seismic attributes such as the chaos attribute (Chopra & Marfurt 2016) and estimate geostatistical parameters that are quantitative and physically meaningful.
Numerous studies have shown evidence that heterogeneous geology can be described as a band-limited, self-similar medium (sometimes referred to as having fractal characteristics). Examples include (i) seafloor bathymetry from multibeam measurements (Goff & Jordan 1988); (ii) exhumed lower continental crust analysed from geological maps (Holliger & Levander 1992); (iii) acoustic and elastic numerical modelling of teleseismic waves recorded by earthquake seismology arrays (Frankel & Clayton 1986) and (iv) analysis of borehole logs from the upper crust (Holliger 1996;Dolan & Bean 1997;Browaeys & Fomel 2009;Cheraghi et al. 2013). Self-similarity means that the statistical properties of the medium do not change with scale. Specifically, medium properties in power-spectral domain will follow an inverse power-law (Dolan & Bean 1997). In this context, band-limited means that there exists a so-called dominant scale length, the scale above which the medium stops showing self-similar characteristics.
There is also a long history of characterizing self-similarity in complex geology using geophysical reflection images. These include (i) investigating partial saturation in freshwater acquifers from ground-penetrating radar images (Irving et al. 2009); (ii) modelling random media heterogeneities to characterize the seismic response of the crust and mantle at different scales (Carcione et al. 2005) and (iii) characterizing the geostatistics of complex turbidite systems from 3-D seismic reflection volumes (Caers et al. 2001). Some studies have explored the link between the spatial statistics of the geological medium and the power spectrum of the reflected wavefield.  present an analytical relationship between band-limited, self-similar random media and a corresponding idealized reflection image. They demonstrate that it is possible to use this relationship to estimate geostatistical parameters characterizing the P-wave velocity heterogeneity, such as the aspect ratio of lateral and vertical dominant scale lengths and the Hurst number (a self-similarity coefficient related to the roughness of the medium). This approach relies on the assumption that the reflection image approximates a so-called primary reflectivity section, an idealized seismic image in depth. Irving et al. (2009) demonstrate that this relationship can recover geostatistical parameters for zero-offset ground-penetrating radar images of shallow, partially saturated acquifers. Scholer et al. (2010) use a similar approach to estimate the correlation structure of P-wave velocity heterogeneity in the crystalline crust from seismic reflection images, including a term to compensate for the theoretical lateral resolution limit of migrated reflection images.
There is currently little published literature investigating the selfsimilar characteristics of the internal structure of mass-transport deposits. Micallef et al. (2008) document scale invariant characteristics of the Storegga Slide, a retrogressive 'megaslide' from the mid-Norwegian margin. They use multibeam bathymetry data to perform a statistical analysis of sub-bodies within the slide and infer that the slide exhibits self-organized critical behaviour. They observe an inverse power-law scaling in their frequency-area distribution and find that headwalls are self-similar from small to large scales. The authors hypothesize that the fractal statistics of submarine landslides could be related to the physics of slope failure.
This study represents the internal structure of mass-transport deposits as a specific type of band-limited, self-similar medium, an anisotropic von Kármán random medium (Von Kármán 1948). In two dimensions the random medium can be characterized by three geostatistical parameters: lateral and vertical dominant scale lengths and the Hurst exponent (roughness). The dominant scale lengths describe the degree of continuity of the medium in horizontal and vertical directions respectively. The Hurst number is a dimensionless parameter related to the degree of self-similarity, which characterizes the roughness or texture of the medium .
The aim of this study is to demonstrate a method to constrain the geostatistics of the internal structure of mass-transport deposits directly from seismic reflection images (after . A further goal is to integrate information from a vertical borehole log, where available. The method is first validated on a synthetic model representing a typical submarine mass-transport deposit scenario, with a synthetic multichannel seismic reflection image and a colocated synthetic vertical borehole. Then, the method is applied to a real data case study from the Nankai Trough, offshore Japan.

M E T H O D O L O G Y
This method inverts seismic reflection images of mass-transport deposits for geostatistical parameters which can characterize their internal structural fabric. This is achieved by assuming the velocity heterogeneity within the mass-transport deposit is a random field defined by an anisotropic von Kármán random medium. Under this assumption it is straightforward to forward model the spatial power spectrum of a corresponding idealized seismic reflection image for a given seismic source spectrum. A likelihood function is defined from the residual between the forward modelled and observed power spectra. Then, a Bayesian Markov Chain Monte Carlo (MCMC) inversion is used to estimate the posterior probability distribution (expected value and uncertainty) for each geostatistical parameter. When borehole log information is available, a constraint on the vertical dominant scale length and Hurst number can be included in the inversion as a prior.

Spatial power spectrum of a random field
Here the velocity field, v, is represented by two components, a smoothly varying background component, v 0 , and a zero-mean, where 1 (i.e. the stochastic component is small relative to the background). In general terms, the background velocity is well resolved by geophysical techniques such as traveltime tomography. At the bandwidth of conventional marine seismic data (approximately 10-100 Hz), however, the small-scale stochastic component generates the vast majority of observed reflections in a seismic image. The small-scale stochastic velocity structure is generally poorly resolved by seismic reflection experiments except perhaps by fullwaveform modelling techniques, which can require significant acquisition effort, model conditioning and computational power, with little measure of uncertainty in the final result.
We make the assumption that the internal heterogeneity (smallscale stochastic structure, v ) of a mass-transport deposit can be approximated as an anisotropic von Kármán random medium.
The normalized 2-D spatial power spectrum of an anisotropic von Kármán random medium (eq. A1) is where c is a normalizing constant, k x and k z are the horizontal and vertical spatial wavenumbers, a x and a z are the dominant lateral and vertical scale lengths and γ is the Hurst number (see Appendix A).

Migrated seismic image
This section follows the methodology presented in  which links the random medium parameters to the 2-D power spectrum of a resulting idealized seismic image, sometimes referred to as a primary reflectivity section (Irving et al. 2009;Scholer et al. 2010). The idealized seismic image is a convolutional, zero-offset, normal-incidence, constant density approximation. The formulation in depth-domain is: where s(x, z) is the idealized seismic image in depth, r(x, z) is the normal-incidence acoustic reflectivity, w(z) is the source wavelet in depth and h(x) is a horizontal filter to account for the lateral resolution of a migrated seismic section (Scholer et al. 2010). The choice of the lateral resolution operator h(x) is arbitrary but should reflect the lateral resolution of the analysed reflection image, which after proper migration is on the order of the dominant wavelength of the seismic source (Chen & Schuster 1999). This study follows Scholer et al. (2010) in using a Gaussian low-pass filter with width (distance between the two points where the filter is 1 per cent of the maximum value) equal to the dominant wavelength: Assuming that variation in density is small relative to velocity, the normal-incidence reflectivity can be approximated as the derivative of the velocity field: If reflections from the smooth background velocity, v 0 , are negligible (i.e. the only contribution to acoustic reflectivity is the smallscale stochastic component of P-wave velocity) and the source wavelet is stationary in depth within the analysis window, the idealized seismic response s(x, z) depends only on the stochastic velocity component, v : The power spectrum of the seismic image can then be related to the spatial power spectrum of the stochastic component by the Fourier transform : where P w is the power spectrum of the source wavelet, w, and P h is the power spectrum of the lateral resolution filter, h. It follows that the power spectrum of the seismic image can be directly related to the random medium parameters by eq. (2): Therefore it is possible to forward model an idealized spatial power spectrum which is comparable to a window of an observed seismic image under the following assumptions: (i) The analysed window of the observed seismic image approximates a noise-free, zero-offset, true-amplitude, convolutional image in depth-domain.
(ii) The stochastic component of P-wave velocity heterogeneity, v , within the analysed window is an anisotropic von Kármán random medium parametrized by a x , a z and γ .
(iii) The geostatistical parameters and source wavelet are stationary over the analysed window.
Only physically realizable models are considered (i.e. dominant scale lengths are non-negative and non-zero).

Borehole log
For geohazard studies borehole logs and cores are often acquired to estimate geotechnical or petrophysical information about the masstransport deposit (Strasser et al. 2011;Dugan 2012). As these logs have spatial power spectra, we can better constrain geostatistical parameters in the direction of the borehole.
Normally, boreholes are approximately vertical, so we can estimate a z and γ independently from a vertical borehole log alone (Browaeys & Fomel 2009). As borehole logs generally directly measure physical parameters we do not need to account for the effect of the seismic source wavelet on the geophysical response of the medium. The 1-D form of eq. (2) is where the exponent is modified for a field with Euclidean dimension N = 1 (see Appendix A).

Inversion for geostatistical parameters
This study uses a Bayesian MCMC approach to invert for the geostatistical parameters. The output of the method is a chain of 'accepted' models whose joint distribution is proportional to the posterior probability density of the model. The chain is sampled using the Metropolis-Hastings algorithm (detailed in Appendix B). Bayesian approaches have the advantage of using prior probability density functions, so prior geological information can be easily incorporated if it can be expressed in terms of the model parameters.
The likelihood function assumes Laplacian errors (doubleexponential distribution) for each observation (Mosegaard & Tarantola 1995): where g i (m) represents the forward modelled power spectrum at a given wavenumber, σ is a parameter proportional to the magnitude of the combined modelling and observation error, N is the number of observations (total number of points in the observed power spectrum) and d obs,i is the observed power spectral density for a given wavenumber k i . For this study, the model parameters considered are the geostatistical parameters (lateral and vertical dominant scale lengths and the Hurst number) and the error parameter. Multiple parallel chains are run to measure convergence and ensure that individual chains are well-mixed and have truly converged (not simply sampling a lowprobability area). For this study, convergence is measured using the Gelman-Rubin statistic (Brooks & Gelman 1998),R, and the weighted mean absolute error (WMAE, Pirot et al. 2017). Details of the convergence measures are given in Appendix B. Chains are generally assumed to have converged whenR < 1.2 for all parameters (Brooks & Gelman 1998). The weighted mean absolute error should oscillate around 1 when the chain is sampling the posterior distribution. To ensure that none of the pre-convergence 'burn-in' samples are included in the posterior distribution, the first half of each chain is discarded from the final posterior distributions.

Seismic reflection image
For the chosen window of the 2-D image (the chaotic mass-transport deposit zone), calculate the following: (i) P obs (k x , k z ), the 2-D spatial power spectrum of the chaotic window (using a 2-D fast Fourier transform).
(ii) P w (k z ), the power spectrum of the seismic source wavelet.
(iii) P h (k x ), the power spectrum of the lateral resolution filter (eq. 4).
Each iteration of the Metropolis-Hastings algorithm proposes a new candidate model, m = [a x , a z , γ , σ ]. For each proposal, forward model the idealized 2-D spatial power spectrum (eq. 8); compute the likelihood of the proposal given d obs = P obs (k x , k z ) (eq. 10) and accept or reject the model according to the acceptance criterion (Appendix B).

Borehole log
Borehole logs generally attempt to directly measure a physical property of the subsurface. Specifically for a sonic log, the measured velocity (or slowness) will include both the background velocity trend, v 0 , and the small-scale stochastic component, v (eq. 1). The background trend must be removed prior to estimation of the dominant scale lengths and Hurst number (Cheraghi et al. 2013). The choice of method for detrending is arbitrary and depends on the complexity of the geology. As this study uses relatively small windows of data from borehole logs, we simply remove the first-order background trend by finding a line of best fit and subtracting it. Borehole logs from different geology may instead require detrending with a more sophisticated approach such as subtracting a loworder best-fitting polynomial. The resulting detrended log should be approximately zero-mean and contain only information from the small-scale stochastic component.
As for the seismic inversion, the power spectrum of the detrended borehole log, P b (k z ), should be computed using a fast Fourier transform.
Each iteration of the Metropolis-Hastings algorithm proposes a new candidate model, m = [a z , γ , σ ]. For each proposal, forward model the idealized 1-D spatial power spectrum (eq. 9); compute the likelihood of the proposal given d obs = P b (k z ) (eq. 10) and accept or reject the model according to the acceptance criterion (Appendix B).  show that under typical experimental conditions the two dominant scale length parameters a x and a z are strongly dependent on each other, such that it may not be possible to resolve each parameter individually from conventional reflection images. However, they show analytically that is possible to reliably estimate the aspect ratio of heterogeneity α = ax az . With an external estimate of one of the dominant scale lengths, for example a z from a vertical borehole log, it should be possible to resolve a x and a z individually.

Seismic image and borehole log
Because the probabilistic inversion approach uses prior probability density functions as an input, we can alter these prior probability density functions to reflect our a priori knowledge of the subsurface. For this inversion, prior probability density functions for the dominant vertical scale length, a z , and Hurst number, γ , are chosen to be Gaussian, with mean and standard deviation calculated from the marginal posterior distributions from the borehole log inversion. The inversion proceeds as for the seismic reflection image.

Synthetic benchmark -buried submarine mass-transport deposit
This synthetic example is designed to benchmark the inversion for a typical marine geohazard survey. The data acquisition simulates a multichannel, marine, towed-streamer acquisition over a chaotic  Table 1. Background elastic parameters and geostatistical parameters for each unit in the synthetic model (Fig. 1). z is the depth below the waterbottom, v P and v S are the P-and S-wave velocities, respectively, and ρ is the density. mass-transport deposit body buried under a water layer and heterogeneous sediment overburden. The aim of this test is to estimate geostatistical parameters from the seismic reflection image with and without an a priori estimate of the vertical dominant scale length from a synthetic borehole velocity log. The model is divided into two layers, a water layer and a sediment layer, both 350 m thick (see Fig. 1). Background elastic parameters and geostatistical parameters for the small-scale stochastic component are given in Table 1. The sediment layer has linearly increasing background P-and S-wave velocity to approximate the effect of increasing compaction with depth on the seismic velocities. It includes a zone with significantly shorter lateral dominant scale length and distinct Hurst number to represent a buried, chaotic mass-transport deposit. Otherwise, the mass-transport deposit zone has the same background elastic parameters as the host sediment layer. The random medium zones are realized on a regular (staggered) 2-D mesh (Ikelle et al. 1993).

Background elastic parameters
This synthetic benchmark simulates a typical 2-D multichannel marine acquisition geometry. The modelled source wavelet is a 40 Hz Ricker. For this synthetic test we use a pseudospectral, isotropic, viscoelastic scheme (Carcione et al. 2005;Carcione 2014) to forward model the seismic reflection response. Sources and receivers are located in the first row of grid points (z = 0 m). For this experiment free surface multiples are not modelled; perfectly absorbing boundary conditions are imposed on all four boundaries of the mesh. P-and S-wave quality factors are set to Q P = Q S = are given in Table 2. In total, 50 shots are modelled which required 25 hr computation time on a quad-core Intel R Core TM i7-6700 3.40 GHz CPU.
As the background velocity model is known and does not vary laterally, the seismic processing follows a basic marine imaging flow, with a pre-stack true-amplitude Kirchhoff time migration (to 60 • maximum angle), outer angle mute (to eliminate refracted arrivals), stack and time-to-depth conversion using the background P-wave velocity model. The image is cut to the full-fold area, with maximum depth equal to the maximum depth in the synthetic model (Fig. 2).

Borehole log inversion
The synthetic P-wave velocity borehole log is shown in Fig. 2. The window analysed is the mass-transport deposit zone between 500 and 650 m depth. For the inversion, uniform priors are used: 0 < a z ≤ 50 m, 0 ≤ γ ≤ 1 and 0 < σ ≤ 2.
The MCMC is run with 12 parallel chains until 1 × 10 4 samples are accepted to the chain ( Table 6). The final Gelman-Rubin statisticR < 1.02 for all parameters. The WMAE begins to oscillate around 1 after approximately 200 accepted samples.
Marginal posterior probability distributions for a z , γ and σ are shown in Fig. 3. Summary statistics of the distributions are shown in Table 3. Both geostatistical parameters are centred closed to their true values.

Seismic image inversion
Two inversions were run on the seismic reflection image, with and without estimates of the vertical scale length Hurst number from the borehole as priors. The synthetic seismic image is shown in Fig. 2. The window analysed is the mass-transport deposit zone highlighted in Fig. 1.
The second inversion (seismic image with borehole) is parametrized as the first, but includes a constraint for a z and γ from the borehole inversion results. The prior probability density functions for a z and γ are Gaussian, with mean and standard deviation from the results of the borehole-only inversion ( Table 3). The priors for a x and σ are uniform, as above: 0 < a x ≤ 500 m, 0 < σ ≤ 2. The priors for a z and γ are truncated Gaussian distributions: for a z , mean 15.9 m and standard deviation 3.5 m (truncated at 0 < a x ≤ 50 m); for γ , mean 0.37 m and standard deviation 0.09 m (truncated at 0 ≤ γ ≤ 1).
Both MCMCs are run with 12 parallel chains until 1 × 10 4 samples are accepted to the chain (Table 6). For the first inversion, the final Gelman-Rubin statisticR < 1.01 for all parameters. For the second inversion, the final Gelman-Rubin statisticR < 1.02 for all parameters. For both inversions the WMAE begins to oscillate around 1 after approximately 100 accepted samples.
Marginal posterior probability distributions for a x , a z , γ and σ , alongside a distribution representing the aspect ratio of heterogeneity, α = ax az , are shown in Fig. 3. Summary statistics of the distributions are shown in Table 3.
With respect to the first inversion (seismic image only) the second inversion (seismic image with constraint from borehole) shows marginal posterior distributions that are closer to the true values.  Table 1). Location of the synthetic borehole is shown in solid red. The mass-transport deposit zone (dashed red outline) shows a more disordered, chaotic seismic character compared to the more stratified unfailed sediments. (b) P-wave velocity log sampled at 0.25 m.

Real data case study -Nankai Trough, offshore Japan
The Nankai Trough (offshore southwest Japan) is an oceanic trench formed by the subduction of the Philippine plate under the Eurasian Plate. Associated accretion, seismicity and slope-steeping have resulted in significant mass-wasting during the last 3 Ma (Strasser et al. 2011;Lackey et al. 2018). A large mass-transport deposit is identified in a 3-D seismic volume (Fig. 4). Here we consider a 2-D profile extracted from the 3-D volume, chosen to show the maximum extent and thickness of the mass-transport deposit. The body has a chaotic internal seismic character, with little visible coherent structure.
The survey acquisition parameters are documented in Table 4 (Uraki et al. 2009). The maximum observed thickness (at the point where the mass-transport deposit intersects the edge of the seismic volume) is approximately 180 m (Strasser et al. 2011).
Also available are logging-while-drilling borehole logs from nearby International Ocean Discovery Programme (IODP) borehole C0018B (Expedition 338, Henry et al. 2012), which penetrates the mass-transport deposit (Fig. 4). No sonic log was acquired, so the gamma ray log is used to estimate the vertical dominant scale length and Hurst number. Whilst the gamma ray log is not a measure of the P-wave velocity, it is sensitive to changes in lithology (specifically shale fraction), which should correlate well with the P-wave velocity. It is expected that both gamma ray and sonic velocity logs should have similar geostatistics within a local interval of a 1-D borehole log.

Borehole log inversion
The gamma ray log from IODP borehole C0018B is shown in Fig. 4. The analysis window is the mass-transport deposit zone between 3235 and 3295 m. For the inversion, uniform priors are used: 0 < a z ≤ 50 m, 0 ≤ γ ≤ 1 and 0 < σ ≤ 2.
The MCMC is run with 12 parallel chains until 5 × 10 4 samples are accepted to the chain ( Table 6). The final Gelman-Rubin statistiĉ R < 1.01 for all parameters. The WMAE begins to oscillate around 1 after approximately 50 accepted samples.
Marginal posterior probability distributions for a z , γ and σ are shown in Fig. 5. Summary statistics of the distributions are shown in Table 5.

Seismic image inversion
Two analysis windows are used on the seismic image, in the downslope and mid-slope parts of the mass-transport deposit (Fig. 4). Both windows have the same dimensions (1000 m by 60 m). The down-slope window is located towards the toe of the mass-transport deposit. The mid-slope window is located relatively further upslope, in the more proximal part of the mass-transport deposit. Two inversions are run for each window, with and without estimates of the vertical scale length a z and Hurst number γ from the borehole log.
The second inversions (seismic image with borehole) are parametrized as the first, but include a constraint from the borehole log inversion results ( Table 5). The prior for a x is uniform, as above: 0 < a x ≤ 500 m. The priors for a z and γ are Gaussian, fit to the marginal posterior probability distributions from the boreholeonly inversion: for a z , mean 5.3 m and standard deviation 1.3 m; for γ , mean 0.41 m and standard deviation 0.13 m.
Marginal posterior probability distributions for a x , a z , γ and σ , alongside a distribution representing the aspect ratio of heterogeneity α = ax az , are shown in Fig. 6 for both zones. Summary statistics of the distributions are shown in Table 5.
With respect to the first inversion (seismic image only), the second inversion (seismic image with borehole) shows better-constrained (lower standard deviation) marginal distributions for a x , a z and γ . The marginal distributions for the down-slope zone show a notably smaller mean a x and α compared to the mid-slope zone, while maintaining similar distributions for a z .

D I S C U S S I O N
This study applies a geostatistical inversion to characterize the internal structure of mass-transport deposits from seismic reflection images, with and without a constraint from a borehole log. We first demonstrate the method on a synthetic model representing a typical buried submarine mass-transport deposit scenario and then on a real data case study from the Nankai Trough, offshore Japan. The method gives probabilistic estimates of lateral and vertical dominant scale lengths and the Hurst number of the internal heterogeneity. To the authors' knowledge, this is the first time that this technique has been validated with a synthetic test on multichannel, stacked seismic reflection data. This is also the first published example to demonstrate how to condition the inversion using priors derived from a vertical borehole log to better constrain the lateral and vertical dominant scale lengths. We suggest that this technique could be a useful tool to better constrain internal structure of mass-transport deposits as it can be applied even to chaotic seismic reflection images of mass-transport deposits, which are common but difficult to interpret using conventional horizon-tracking methods.

Synthetic inversion results
For the inversion performed on the synthetic seismic image with uniform priors, the estimated aspect ratio of heterogeneity, α = ax az , is close to the true model value (Fig. 3). This result is expected from previous studies, which suggest that the 2-D power spectrum (equivalently the 2-D autocorrelation function) is strongly sensitive  to the aspect ratio of heterogeneity rather than to the individual dominant scale lengths or the Hurst number (Irving et al. 2009;Scholer et al. 2010). This synthetic test, however, shows relatively good resolution of separate lateral and vertical scale lengths from the seismic image alone. The Hurst number is still poorly constrained. Repeating the inversion with priors for vertical scale length and Hurst number estimated from a synthetic borehole log improves the result, but only slightly. This is in constrast to the conclusions of , who predict that the 2-D power spectrum should be sensitive only to the aspect ratio of heterogeneity. Our result is likely because the bandwidth of the seismic source overlaps both the 'white noise' and self-similar parts of the random medium in power-spectral domain for this test. Another contributing factor is that this is a synthetic experiment. Seismic images created from field data contain noise from (i) environmental noise (ii) instrument noise (iii) multiple arrivals and (iv) processing artefacts. Future studies should investigate the reliability of this method to discriminate lateral and vertical dominant scale lengths under a range of noise conditions and source bandwidths, with respect to the spatial power spectrum of the medium.

Nankai Trough case study inversion results
For the Nankai Trough experiment we consider two identically sized data windows, a down-slope zone and a mid-slope zone (Fig. 4). The down-slope zone is located towards the toe of the mass-transport deposit. The mid-slope zone is more proximal, located towards the middle of the deposit. The seismic character in both windows is chaotic, lacking laterally coherent seismic reflectors. First, we invert for the geostatistical parameters in both windows with uniform priors (Fig. 6). In the down-slope zone, the aspect ratio of heterogeneity, α, is significantly smaller than in the mid-slope zone. Including priors for a z and γ based on the nearby IODP borehole C0018B (Fig. 4), we still see a reduction in α from mid-slope to down-slope, but we see the distributions for lateral and vertical dominant scale lengths, a x and a z , are much better constrained.
Mass-transport deposits often show extensional structures near the headwall, little deformation in the central translational zone and compressional structures in the toe region, where the flow may be confined (Fig. 7). The observed reduction in lateral dominant scale length from mid-slope to down-slope is consistent with this interpretation of the mass-transport deposit. More compression will result in increased stratal disruption, giving a shorter lateral dominant scale length compared to relatively undeformed sediments. Alternatively, the reduction is lateral dominant scale length could be due to progressive down-slope deformation of the mass-flow (Lucente & Pini 2003). Both models could explain the reduction in lateral dominant scale length and aspect ratio of heterogeneity.
The velocity heterogeneity within the mass-transport deposit should be closely related to lithological heterogeneity. For masstransport scenarios, this heterogeneity could be predominantly due to included megaclasts, intact blocks or intense folding from stratal disruption. The observed reduction in lateral scale length is consistent with most conceptual models of the variation in internal structure from proximal to distal within the depositional part of mass-transport deposits (e.g. Bull et al. 2009, see Fig. 7).

Internal structure from geostatistical parameters
How should these geostatistical parameters be interpreted in the context of mass-transport deposits? These parameters are abstract and set in terms of a statistical model, not in terms of geological structure. We suggest that the dominant scale lengths can be proxies for relative deformation from both mass-transport processes and tectonic stresses. Increasing deformation (e.g. folding from compression, reduction in size of intact blocks due to progressive disaggregation) should reduce the lateral dominant scale length and also the aspect ratio of heterogeneity.
Here we only consider heterogeneity of the P-wave velocity field, as we believe this should capture much of the geological heterogeneity that controls the seismic response. In fact, this method could be used to constrain any kind of geological heterogeneity, so long as it can be related to the acoustic impedance (the idealized seismic image approximation only models normal-incidence reflections). For the mass-transport deposit case, for example, one could consider the geological medium as a mixture of two component lithologies with distinct acoustic impedances (e.g. matrix and clasts). Thus estimating the geostatistical parameters can inform the geostatistics of the geology directly.   Bull et al. 2009). Note increasing deformation due to confinement towards the toe of the slide. (b) Illustration of two mechanisms for reducing the lateral dominant scale length by mass-transport -disaggregation of large coherent intact blocks and stratal disruption of soft sediments. In general increased deformation will result in a decrease in lateral dominant scale length (and aspect ratio of heterogeneity). (c) Outcrop example of variation in lateral dominant scale length due to a reduction in size of included megaclasts (Vernasso Quarry, NE Italy). For unfailed sediments, one would expect very long lateral dominant scale lengths due to the presence of laterally continuous beds. After failure, sediments may become deformed due to shearing and disaggregation, reducing the lateral dominant scale length. Therefore the lateral dominant scale length is a useful structural parameter that can be a proxy for lateral shortening from deformation. The vertical dominant scale length is more closely related to the average thickness of beds, and therefore may be less affected by mass-transport.

Limits in generalization
Using a synthetic example we show that an idealized seismic image approximation (Section 2.2.1) is valid for one multichannel marine seismic experiment, with a specific overburden and seismic character. This allows a computationally inexpensive inversion method (Table 6) to estimate random medium parameters from a window of a reflection image. The validity of the approximation will depend on the local geology and on the seismic imaging performed. Multiple scattering, attenuation and seismic noise will all reduce the validity of the idealized seismic image approximation.
The method presented in this study uses the spatial power spectrum to evaluate random media models and to estimate the misfit between a corresponding theoretical and observed seismic reflection image. For a given spatial power spectrum there exist infinite physical realizations of the corresponding random medium. It is important to note that this method only constrains the statistics of the heterogeneity, not the medium properties directly. It is possible that there are better statistical representations, especially for small window sizes which may suffer from edge-effects from the fast Fourier transform when calculating the power spectrum. Some previous studies have used the autocorrelation function instead (e.g. Irving et al. 2009;Scholer et al. 2010).
This study only considers 2-D seismic profiles. Mass-transport is an inherently 3-D geological process, so strong lateral heterogeneity observed in the plane of the profile implies that strong heterogeneity perpendicular to the profile is also likely. This 3-D heterogeneity could generate strong out-of-plane reflections. For a chaotic seismic reflection image, it may be impossible to identify or remove these out-of-plane reflections during imaging or interpretation. It is presently unclear how the results of the inversion may be affected if these spurious reflections contaminate the analysis window. This topic could be addressed in a future 3-D numerical modelling study by performing geostatistical inversion on 2-D profiles which include out-of-plane reflections.
Is the anisotropic von Kármán random medium a suitable statistical representation of the internal structure of mass-transport deposits? There exist many studies suggesting that geology in general has fractal-like properties (band-limited self-similarity; e.g. Goff & Jordan 1988;Turcotte 1997;Browaeys & Fomel 2009;Nelson et al. 2015). There exist, however, few studies investigating the fractal properties of internal structure of mass-transport deposits (Micallef et al. 2008). Analysis of mass-transport deposits in outcrop is necessary to determine whether an anisotropic von Kármán random medium could be a broadly applicable geostatistical model.
The formulation used in this study (eq. 2) assumes no dominant dip direction. This could be reasonable for mass-transport deposits deposited in the deep ocean, for example, but not if there has been post-depositional deformation from tectonics. In future work it should be straightforward to include dominant dip direction as an extra parameter in the inversion (see Yuan et al. 2014, for an example).

C O N C L U S I O N S
In this study we demonstrate a method to invert for geostatistical parameters (lateral and vertical dominant scale lengths and Hurst number) which can describe the internal structure of mass-transport deposits from chaotic multichannel seismic reflection images and borehole logs. This approach assumes that the internal structure can be approximated as an anisotropic von Kármán random medium . The results are probability distributions which provide the expected value and uncertainty of each geostatistical parameter.
The method is first validated on a synthetic scenario containing a buried chaotic body, representing a submarine mass-transport deposit, imaged with a typical multichannel marine seismic acquisition and penetrated by a synthetic borehole. When the seismic image is inverted with a constraint from the borehole, lateral and vertical dominant scale lengths and Hurst number can be recovered.
The method is then applied to a real data case study from Nankai Trough (offshore Japan). The data considered are a seismic reflection profile and the gamma ray log from a borehole which penetrates a thick mass-transport deposit. We see a reduction in lateral dominant scale length from mid-slope to down-slope part of the masstransport deposit. This is consistent with progressively increasing deformation due to disaggregation or compression towards the toe of the slide.
Geostatistical inversion is a useful new tool to constrain the internal structure of mass-transport deposits from seismic reflection data. The geostatistical parameters can be used to validate conceptual models of internal structure and as a proxy for varying strain or degree of deformation in different domains of the slide, even when the seismic image appears chaotic or reflections lack the continuity required for horizon-tracking approaches. The lateral dominant scale length in particular could be a good proxy for strain history, as it is strongly related to the degree of sediment deformation and stratal disruption.

A C K N O W L E D G E M E N T S
The quality of this manuscript was greatly improved thanks to constructive comments from the reviewers Guilliame Pirot and Jules Browaeys and the Editor, Frederik Simons. We thank Umberta Tinavella, José Carcione and Aldo Vesnaver for feedback on an early version of this manuscript. Thanks to Philippe Cance and Davide Gei for providing the pseudospectral numerical modelling code and for discussion on the topic and to Greg Moore (SOEST) for providing the Kumano 3-D seismic data set. Jonathan Ford was supported by a Marie Curie Doctoral Fellowship through the SLATE Innovative Training Network within the European Union Framework Programme for Research and Innovation Horizon 2020 under Grant Agreement No. 721403.

B2 Convergence criteria
The Markov Chain is guaranteed to sample the posterior distribution once the chain has converged (after the so-called 'burn-in' period). One problem with MCMC methods is determining when the 'burnin' period has finished, that is, after which point to start considering samples as part of the posterior distribution. Estimating convergence is important as oversampling the chain increases the computation time, whilst undersampling the chain may bias the chain towards the starting values and not properly sample low probability regions of the posterior.
For this study, the Gelman-Rubin statistic (sometimes called the scale-reduction factor),R, is calculated (Brooks & Gelman 1998). This involves running several chains in parallel and comparing the in-chain variance to the between-chain variance for each parameter in the model. For m chains of length 2n accepted samples, W is the mean variance each chain, B is the variance of the mean of each chain, V h = B + W (n−1) n ,R = V h W . It is commonly considered that chains have converged for a parameter whenR < 1.2 (Brooks & Gelman 1998).
The weighted mean absolute error (WMAE) is given as where N is the number of observations (number of wavenumber pairs). The WMAE should oscillate around 1 when the chain is properly sampling the posterior distribution (Pirot et al. 2017).