Radial anisotropy in Valhall: ambient noise-based studies of Scholte and Love waves

We


I N T RO D U C T I O N
The presence of the seismic anisotropy in the Earth's crust and mantle is well established. Mantle anisotropy is well studied and is believed to principally reflect the lattice preferred orientation of olivine (e.g. Christensen & Lundquist 1982;Montagner & Anderson 1989b). The crustal anisotropy can be caused by a variety of mechanisms including mineral orientation, fine layering within sedimentary or magmatic rocks, or the preferred orientation of faults or cracks (e.g. Duret et al. 2010;Moschetti et al. 2010;Jaxybulatov et al. 2014;Montagner 2014;Maupin & Park 2015;Mordret et al. 2015).
S-wave velocities inferred from the Rayleigh waves often differ from those obtained with the Love waves. Such incompatibilities are generally considered as a robust diagnostic for the presence of anisotropy in the crust and upper mantle, which is commonly called as radial or polarization anisotropy (e.g. Anderson 1961;Schlue & Knopoff 1977;Nakanishi & Anderson 1983;Montagner 1985;Nataf et al. 1986;Ekström & Dziewonski 1998) or Rayleigh-Love wave discrepancy. The radial anisotropy is a property of a medium in which the speed of the wave depends on its polarization and corre-sponds to a vertical transversely isotropic (VTI) medium. Seismic anisotropy is becoming more and more important in exploration seismology (e.g. Muyzert & Kommedal 2002;Helbig & Thomsen 2005;Hatchell et al. 2009;Tsvankin et al. 2010;Mordret et al. 2013c).
With the development of ambient noise surface-wave tomography (ANSWT, e.g. Shapiro & Campillo 2004;Shapiro et al. 2005;Campillo & Sato 2011;Ritzwoller et al. 2011), surface wave dispersion curves can be measured at periods short enough to construct S-wave velocity models at crustal depths and to infer the seismic anisotropy in the crust. The range of the tomography scale based upon surface wave varies from several thousands of kilometres (Nishida et al. 2009) to a few hundreds of metres. However, there are only few studies based upon the ANSWT that investigate the depth ranges from upper tens of metres to several hundreds of metres. To investigate the depth down to 1 km, it is required to have surface waves with the wavelengths smaller than 1 km, implying that for frequencies between ∼0.5 and ∼2 Hz the phase velocities must be below 1000 m s −1 . These frequencies are still lower than the frequencies in most of the existing type of active seismic sources.

1525
In exploration seismology, studies to image the subsurface using surface waves have been developing recently (e.g. Stewart 2006;Dellinger & Yu 2009; de Ridder & Dellinger 2011). In particular, the Valhall passive seismic data sets were used extensively for the ANSWT (e.g. de Ridder & Biondi 2013;Mordret et al. 2013a,b,c;Mordret et al. 2014a,b). A high-resolution isotropic 3-D S-wave velocity model down to 600 m has been derived by Mordret et al. (2014a), using ANSWT of the fundamental mode of Scholte wave. The Love wave group velocity maps were constructed at life of field seismic (LoFS) by de Ridder (2014). Although, the Scholte and Love waves analysis is done separately and azimuthal anisotropy is estimated using Scholte wave at Valhall, we propose joint analysis of Scholte and Love waves to obtain radial anisotropy as it differs from the azimuthal anisotropy and depends upon the polarization of waves.
In this paper, we present a study of the radial anisotropy at the depths above 600 m. After presenting the data and giving a brief introduction of methods, the 2-D phase-velocity maps of Scholte and Love waves at different frequencies are constructed using the eikonal equation (Aki & Richards 1980;Lin et al. 2009). The local dispersion curves are computed from the Scholte and Love waves phase-velocity maps and are jointly inverted by using the inversion method based on the neighbourhood algorithm (NA; Sambridge 1999) to obtain local 1-D velocity profiles. All 1-D profiles are then combined to obtain a 3-D radially anisotropic S-wave velocity model of the Valhall overburden.

DATA A N D N O I S E C RO S S -C O R R E L AT I O N S
The Valhall LoFS network is a permanent ocean-bottom cable array made of 2320 four-component sensors (a three-component geophone and a hydrophone) installed on the seafloor in the North Sea ( Fig. 1), above Valhall reservoir. We use 400 min (∼6.5 hr) of continuous ambient noise records on vertical (Z), south-north (N) and west-east (E) component of geophones for this study. The data set was recorded at 250 samples per second with the low-cut filter that removes almost all energy at periods longer than 2.5 s. The interstation distance between the sensors along cable is 50 m and intercable distance is 300 m.
For three component sensors more than 10 millions of interstation cross-correlations (CCs) are computed between all possible pair of sensors for vertical-vertical (ZZ) component along with north-north (NN), east-east (EE), north-east (NE) and east-north (EN) components combinations. The noise is processed station by station, before computing CC we apply mean and trend removal to the signal, and then we do amplitude spectral whitening of the signal between 0.3 and 100 Hz. Although the signal at Valhall is clean but one-bit normalization is also applied to the data set. Finally, we follow the procedure described in Mordret et al. (2013a) for CCs computation and also partially follow the workflow discussed in Bensen et al. (2007). Scholte waves can be obtained directly from ZZ CCs, while for Love waves, we need to rotate the CCs between EE, EN, NN and NE component in radial and transverse directions. The transverse-transverse (TT), radial-radial (RR), transverseradial (TR) and radial-transverse (RT) CCs between each pair of stations are obtained using tensor rotation method discussed in Lin et al. (2008). The TT and RR components exhibit Love and Scholte waves, respectively. Here, only ZZ component CCs are used for Scholte waves because of better signal-to-noise ratio (SNR).   Fig. 2 shows the virtual gathers of ZZ and TT components where, ∼2 690 040 CCs for each component are stacked that fall within a 10 m interstation distance bin. These average CC gathers are filtered between 0.3 and 3 Hz. Stacking the CCs together increases the SNR dramatically and allow us to observe different modes for Scholte and Love waves and can be used in depth inversion to estimate average velocity model. (a) F-K spectrum of ZZ component virtual gathers, the fundamental and first higher modes of Scholte wave are clearly visible, (b) the picks of these dispersion curves in F-K spectrum. (c) Average phase-velocity dispersion curves, the magenta and blue dots are the picks from F-K spectrum of fundamental and first higher modes, respectively and curves with error bars are the third-order polynomial fit for the corresponding picks.

Measurement of average dispersion curves for Scholte and Love waves
To obtain the average dispersion curve, we use the average virtual gather for Scholte wave (see Fig. 2a). A similar procedure has been applied recently by Lin et al. (2013), Nishida (2013) and Mordret et al. (2014a). Fig. 3 shows the procedure to compute the average dispersion curves of Scholte waves. We perform frequency wavenumber (F-K) analysis (Gabriels et al. 1987) of virtual gather to measure the average dispersion curves of the fundamental mode and the first overtone (Fig. 3a). The dispersion curves are measured by picking the maximum energy for each mode (Fig. 3b) and phasevelocity dispersion curves can be estimated for a frequency f and a wavenumber k using the following equation: where m is the mode number and c is phase velocity. We then fit these measurements with third-order polynomial to obtain smoothed versions of the dispersion curves (solid lines in Fig. 3c). Fig. 4 shows the procedure to obtain average dispersion curves of Love waves; we apply the same procedure to TT component (see Fig. 2b) as developed for ZZ component. The average phase-velocity dispersion curves with their error bars are shown in Fig. 4(c). Error bar shows the confidence interval of data or deviation along a curve. The fundamental mode of the Love waves (and in a lesser extent the fundamental mode of the Scholte waves) seems to present a double branch at high frequencies and high wavenumbers. This is due to the heterogeneity of the medium (the Love wave velocity maps show clearly a dichotomy with high velocity in the north and low velocity in the south.) We chose to pick the fastest velocity dispersion curves for both Scholte and Love waves.

Inversion of the average dispersion curves: isotropic velocity model
We invert Scholte and Love waves average dispersion curves with the method described by Mordret et al. (2014a) based on an NA (Sambridge 1999;Xie et al. 2013). The NA is a Monte Carlo global direct search technique to sample the whole model space. Model space is bounded by a priori information that includes the range in which the individual model parameters are allowed to vary. In the inversion, the model is 1-D layered S-wave velocity/depth profile, where the parameters for inversion could be both thickness and velocity of each layer. In this paper, only the S-wave velocity is inverted because of the low sensitivity of surface wave to P-wave velocity and density. Therefore, the P-wave velocity and density are not considered as free parameters but taken by scaling from the empirical relation: V p = 1.16V s + 1.36 (Castagna et al. 1985) and ρ = 1.74(V p ) 2 with V p in kilometres per second (Brocher 2005). The S-wave velocity of Valhall overburden is modelled using a power-law expression of the form (Wathelet 2004): where V s is the S-wave velocity, d is the depth, V 0 is the velocity at the seafloor, α is the power-law parameter and d 0 is the water depth. Valhall bathymetry is nearly flat therefore, we take constant d 0 = 70 m, requiring only two parameters V 0 and α to fit eq. (2). This equation is valid when d ≤ 600 and we choose to invert the velocity (V s = V n when, d > 600) of half-space at the 600 m depth, as a third parameter of inversion. The NA is applied in different steps to invert average dispersion curves. In the first step, the method generates n 1 random models inside the model space. Then, a mesh of Voronoi cells (the nearest neighbour portion of space associated with each model) is created. Using Herrmann & Ammon (2004) algorithm, the theoretical dispersion curves are computed for fundamental and first higher modes of Sholte and Love waves and misfit between theoretical and observed dispersion curves is calculated in the corresponding cell. The misfit of Scholte and Love waves dispersion curves are defined as the ratio between the area of synthetic dispersion curves and the area of measured dispersion curves with its uncertainties (see Mordret et al. 2014a for definition). In the second step, we choose the best n c cells with minimum misfit and generate n s new models within each of these cells. The set of new Voronoi cells are generated using all new and previous models. The new misfit is calculated and new n c cells are chosen to be resampled.
The last step is repeated and algorithm stopped after n i iterations. When the model space is considered to be well sampled, the ensemble of tested models and their respective misfits can be used statistically to characterize the solution of the inverse problem. The selected parameters bounds are [150, 550] m s −1 for V 0 [0.08, 0.28] for α and [550, 1350] m s −1 for V n . We run NA with n 1 = 10 000, n s = 1000, n 5 = 5 and n i = 8, that is, a total of 50 000 models are tested. For that, we use computer with 24 cores and 45 GB internal memory and Matlab as a platform to run the NA code written in C language that took around 15 min to run the inversion for a point. Finally, we select 1000 the best-fit models and take average of these 1000 models for defining the final velocity profile. See Appendix A, where we show that the choice of taking the average of these 1000 cells does not change the final model. The averaged models and corresponding dispersion curves obtained with isotropic-parametrization inversion of the Scholte and the Love waves are shown in Fig. 5.

Average radial anisotropy
The 1-D velocity profiles obtained with the inversion of Scholte and Love waves are different (Fig. 5c). We performed several tests and found that the measured dispersion curves for these two types of waves cannot be predicted simultaneously with a single S-wave velocity profile. Fig. 6 shows the isotropic joint inversion of Scholte and Love waves. It can be clearly seen that synthetic dispersion The inverted and observed dispersion curves of Love wave. (c) V SH and V SV are the averaged best-fit models, there is a difference between the two models above 220 m V SH < V SV and below 220 m down to 600 m V SH > V SV . curves do not fit observed dispersion curves well in isotropic joint inversion. Suggesting the presence of anisotropy in the medium.
In a most general case of anisotropic media, seismic wave propagation is sensitive to the full elastic tensor (21 parameters), the density and the attenuation. Often, medium with hexagonal symmetry are considered, where the anisotropy is described by five elastic parameters A, C, N, L and F (Love 1927). When the axis of symmetry is vertical (so-called VTI medium), these parameters can be related to vertically and horizontally polarized compressional wave speeds (V PV , V PH , respectively), the vertically and horizontally polarized shear wave speeds (V SV and V SH , respectively) and with a parameter F = η (A -2 L), where η = 1 for an anisotropic medium. This type of medium is called as radially anisotropic (Anderson 1961;Montagner & Anderson 1989a;Babuska & Cara 1991). It is assumed that Rayleigh or Scholte waves are mainly sensitive to V SV and Love waves to V SH (e.g. Anderson & Dziewonski 1982;Montagner & Nataf 1986). The difference between these waves is called radial anisotropy that also corresponds to transverse isotropy and can be represented in percentage by a parameter ξ (e.g. Huang et al. 2010;Jaxybulatov et al. 2014;Tomar et al. 2015): V S is the isotropic component of the S-wave velocity of such medium and can be represented by the Voigt average (Dziewonski & Anderson 1981;Montagner 2014): After considering the isotropic inversion results, we therefore, introduce the radial anisotropy in the top 600 m of our velocity model to explain simultaneously the observations of Scholte and Love waves at Valhall.

Inversion of the average dispersion curves: anisotropic velocity model
We jointly invert the Scholte and Love waves dispersion curves, using NA with a different parametrization than discussed above. Based on results shown in Section 3.2 and Fig. 5(c), we consider two layers of anisotropy: a layer with negative S-wave anisotropy (V SV > V SH ) above 220 m and a layer with a positive S-wave anisotropy (V SV < V SH ) below 220 m. To reduce the number of inverted parameters, we make a simplification by scaling V PV to V SV and V PV to V SH using Castagna's relation discussed above, and take η = 1. This simplification is not fully physical, because in real-mineral assemblages, the S-wave anisotropy is different from P-wave anisotropy and η differ from unity (Babuska & Cara 1991;Erdman et al. 2013). Also, it has been shown that this simplification might result a slight underestimation or overestimation of the anisotropy (Anderson & Dziewonski 1982;Montagner 1985;Xie et al. 2013). However, the difference between the anisotropy obtained using two coefficients and with five coefficients is small and lies within the uncertainties (Xie et al. 2013). Here, the anisotropic inversion is performed with four parameters; V 0 and α for the power law and ξ 1 and ξ 2 two anisotropic parameters for two different depth. The velocity of the half-space below 600 m depth is fixed on the basis of isotropic inversion to avoid over parametrization. The inversion is performed in several steps, first, a model is randomly chosen by NA and is attributed to V SV . We compute a synthetic Scholte wave dispersion curve using computer program in seismology package (Herrmann & Ammon 2004) and the misfit m r between the synthetic and observed dispersion curves is stored. Then, the same model is perturbed using two random anisotropy parameters to obtain V SH . A synthetic Love wave dispersion curve is computed and misfit m l is stored. The final combined misfit M is computed as: We choose this Love wave misfit with 0.4 weight because the Love wave dispersion curves have smaller number of points fitted with two more parameters than the Scholte waves. We test 50 000 1-D profiles and select N = 1000, the best-fit models that are averaged for defining the final anisotropy and 1-D velocity model. The average anisotropy at each location (each gridpoint) is computed as: The V s isotropic model at each location is computed in a similar way: The results of the anisotropic inversion for the average 1-D model are shown in Fig. 7. We can see that the synthetic dispersion curves fit the observations well, indicating that the parametrization is sufficient to model the average velocity profile accurately down to 600 m depth. These average velocity structure could, therefore, be used as a reference model for 3-D depth inversion. Here, we find the anisotropy in the shallow part ξ 1 = −2.79 per cent with 0.9 per cent of uncertainty and in the deeper part ξ 2 = 9.7 per cent with 1.89 per cent of uncertainty for the average dispersion curves. The uncertainties are computed as standard deviation of the anisotropy from 1000 best-fit models. The negative anisotropy in the shallower part is really low for average model. Three parameters (considering only one layer of anisotropy in the deep part) could also fit the data well. In that case, to perform inversion with four parameters might not significantly improve the results. In order to investigate the impact of inversion with three and four parameters, we use F-test in Section 4.3.

A M B I E N T N O I S E S U R FA C E WAV E T O M O G R A P H Y O F S C H O LT E A N D L OV E WAV E S
After constructing the average 1-D velocity model, we proceed with 3-D tomography. At this stage, we retain only the fundamental-mode measurements and do not use the first overtone data because it is too noisy, when considering individual station pairs. The surface wave tomography is done using a common three-step approach (e.g. Ritzwoller et al. 2011). In the first step, the frequency-dependent phase traveltime is measured for all noise CCs computed of every individual virtual source. In the second step, the phase traveltime is inverted to construct 2-D phase-velocity maps for Scholte and Love waves at different frequencies using eikonal tomography technique (Lin et al. 2009). In the final step, we invert the local dispersion curves computed from 2-D phase-velocity maps for depth structure. The regionalized dispersion curves for every cell of the grid are inverted for a local 1-D S-wave velocity profile. The ensemble of these 1-D velocity profiles is then combined to obtain a final 3-D subsurface structure.

Scholte wave phase-velocity maps
The method discussed in Mordret et al. (2013b) is used here to construct Scholte wave phase-velocity maps in the following steps: first, the frequency-dependent phase traveltimes are computed from each interstation noise CC and interpolated on to a regular grid of 50 × 50 m across the whole Valhall array, using spline in-tension interpolation (Wessel & Bercovici 1998). The in-tension coefficient takes care of smoothness of the interpolated surface between the data points. Secondly, we compute the slowness maps for all 2320 receivers by taking the gradient of interpolated phase traveltime in this period range. Thirdly, the 2-D slowness maps are constructed by averaging all 2320 maps and its uncertainty is estimated by taking their standard deviation. Finally, we invert the slowness map to get a 2-D phase-velocity distribution for the whole Valhall array at a particular frequency. Fig. 8 shows results obtained at periods of 0.7 and 1 s. The uncertainties in phase velocities for the whole grid are very low (<20 m s −1 , see Figs 8b and d). Also, these uncertainty maps are used in the depth inversion to drive 3-D structure of Valhall. We see high-velocity anomaly indicated with the dotted curve in 0.7 s map (see Fig. 8a) that corresponds to the shallow palaeochannels as interpreted by BP author (see, Sirgue et al. 2010). In the southeast part of the area, a big palaeochannel is clearly visible in the 1 s map (Fig. 8c, indicated with the rectangle) also interpreted by BP author (see, Sirgue et al. 2010). The phase-velocity maps of Scholte wave are computed from 0.7 to 1.6 s at an interval of 0.1 s.

Love wave phase-velocity maps
We use the TT CCs to measure the Love wave phase traveltimes. From ZZ and TT average correlation gathers (Fig. 2), it can be seen that Love wave signal is weak in comparison with the Scholte wave. Therefore, for the Love wave fundamental-mode tomography, a slightly different processing from Scholte waves is designed. To increase the SNR for Love wave, for a particular virtual source the CCs from every four closely located sensors, which lie within the interstation distance of 100 m and 180 • of azimuth angle, are stacked before applying the eikonal tomography. This technique is well known in exploration seismology to increase semblance for velocity analysis (Yilmaz 2001). The horizontal resolution analysis at Valhall performed by Mordret et al. (2013b) had shown that the features could be resolved at the scale of 300-400 m. This indicates that by stacking CC from close pairs of receivers, we do not compromise the horizontal resolution in our data set. Fig. 9 shows the individual virtual source gather of 595th station (position of this station is indicated with a star in Fig. 1) before (Fig. 9a) and after (Fig. 9b) trace stacking. In the figures, black lines represent the correlations that pass our SNR criteria and for these correlations, we compute spectral phase. It can be seen that after stacking the SNR as well as the phase information are improved (Fig. 9d).
After stacking closely located traces, we apply the same procedure to compute phase-velocity distribution for Love waves as was applied for Scholte waves. Fig. 10 shows Love wave phasevelocity maps at periods of 1 and 1.4 s. Similar to the Scholte waves, the uncertainties of Love wave phase-velocity maps are very low (<20 m s −1 , see Figs 10b and d) for whole grid. In Love phasevelocity map at 1.4 s, the big palaeochannel is clearly visible (indicated with the rectangle in Fig. 10c). The phase-velocity maps of Love wave are computed from 0.8 to 1.5 s at an interval of 0.1 s. These phase-velocity maps along with the Scholte waves are used to estimate the local dispersion curves.

Depth inversion of the Scholte and Love waves local dispersion curves
We create local dispersion curves from the 2-D phase-velocity maps at different periods for Scholte and Love waves fundamental mode   and invert them jointly using the method described in Section 3 with the anisotropic parametrization. The local dispersion curve at a particular location is sensitive to 1-D S-wave velocity structure beneath that point. At every location, a total of 25 000 models are tested with n 1 = 10 000, n s = 500, n c = 5 and n i = 5. The procedure of inversion is defined in Section 3.2 in details. As mentioned in Section 3.2 for a point inversion, it takes 15 min but for local inversion it takes 7 min for one point and we run inversion for 14 550 local points. Among 25 000 models, we select 1000 best-fit models and take their average models to get the final 1-D vertical profile and depth-dependent anisotropy using eqs (6) and (7). In Appendix A, we compute the histograms (Fig. A1) of these 1000 misfit at four different locations to show that the choice of taking the average of these 1000 cells does not change the final model.
The average anisotropy results from Section 3 show that in the shallow part of Valhall the anisotropy is relatively low. Therefore, we test how statistically significant was adding this extra parameter in the inversion. For that, we perform inversion for whole grid by considering two different parametrizations: one is with four parameters, as described above in Section 3 and another is with three parameters, where we consider only one parameter of anisotropy for the layer below 220 to 600 m. We use the standard F-test (Stein & Gordon 1984;Trampert & Spetzler 2006;Pollitz & Snoke 2010) based on χ 2 distribution defined as: where, N is the total number of data points (the discrete phasevelocity measurement at different frequencies in our example), d is data, d m is the data prediction, m is the number of independent parameters and σ is uncertainty in the velocity measurement. The ratio of χ 2 from two different models is F distributed, that allows to compute a probability to test the two χ 2 are significantly different or this is just a random fluctuation in the data. If a data set is fitted with two different independent parameters p (inversion with three parameters) and q (inversion with four parameters) and q > p, then the second model should fit the data better (e.g. the inversion with four parameters should fit the data better than the inversion with three parameters in our case) or we could say that the χ 2 (q) should be less than χ 2 (p). To test the increase in number of parameters makes a significant improvement in fitting the observed data with synthetic or not, we use the statistical parameter F: To get the significance of the two parametrization from the test, we get the probability, P f (F, v 1 , v 2 ), where v 1 = q − p and v 2 = N − q are the degree of freedom. The F-value is computed for different models by using eq. (9). Thus, if we have P f = 0.01 that means we have only 1 per cent risk that the improvement in the fit is due to chance, while 99 per cent chance is the extra parameter is warranted. First, we performed inversion with three parameters (p) and then with four parameters (q). Appendix B shows the misfits between the inverted and observed dispersion curve are estimated using eq. (5) at each location for both cases (Fig. B1) and it can be seen that using four parameters misfit is improved in comparison of using three parameters. The significance of the improvement in misfit is measured using F-test in Appendix C. The probability distribution indicate that we have 99 per cent chance the inversion with the four parameters bring more information about the subsurface at Valhall than the inversion with the three parameters (Fig. C1).

Modified parametrizations for the large palaeochannel
For most of the region, the misfit (Fig. B1b) is small (<0.7), implying that the four parameters are sufficient to explain the data. However, in the southeastern part of field, the misfit remains very high. This is due to the presence of a large palaeochannel with very high velocities at depths between 180 and 240 m (Sirgue et al. 2010;Mordret et al. 2013a,b;de Ridder 2014). Therefore, we need to modify our parametrization to invert such a high-velocity anomaly at that depth. For this, we use the power law eq. (2) with the addition of the Gaussian-bell layer for varying depth, width and height to model a high-velocity anomaly: where V l is the velocity perturbation, d l is the depth and σ l is the thickness of the layer. After few tests with inverting the average dispersion curves of the palaeochannel region, we fix the depth and the thickness of anomaly to 185 and 87 m, respectively, and keep only one free parameter V l . Therefore, our final anisotropic parameterization for the palaeochannel area includes five parameters: V 0 , α, ξ 1 , ξ 2 and V l . Fig. 11 shows the comparison of the depth inversion results at one point located within the palaeochannel region (the point is indicated with black bullet in Fig. 11g), using parametrizations including four and five parameters. We clearly see that adding the fifth parameter improves the fit between the dispersion curves significantly. The final 3-D velocity model is constructed by merging the inversion results of the Gaussian layer parametrization for the palaeochannel region with the results based on four parameters elsewhere. In Appendix D, we compute the misfit map (Fig. D1) of the hybrid model where, it can be seen that the misfit anomaly disappears after the hybrid inversion taken into account.

I N V E R S I O N R E S U LT S A N D D I S C U S S I O N
We compute the 3-D Voigt average V s isotropic model for Valhall using eq. (7) and present two depth slices of velocity anomalies at 120 m (Fig. 12a) and 180 m (Fig. 12b) depths below sea level. Fig. 13 shows the vertical profiles along the section that are indicated on the depth slice at 100 m (Fig. 13e).
In the depth slices, we can clearly see some geological feature such as the shallow palaeochannels (indicated with the black dotted lines in the velocity slice at 120 m depth), the low-velocity anomaly and the large palaeochannel (indicated with the dotted circle and a rectangle in the velocity map at 180 m depth, respectively). The palaeochannels are the quaternary structure at Valhall. In the profile DD the depth of the large palaeochannel is retrieved very well and match with the results obtained in Sirgue et al. (2010) and Mordret et al. (2014a). The shallow palaeochannel (Fig. 13a profile AA ) extend up to the seabed, due to poor resolution in the shallower part and limited higher frequency measurement in this study.
With the limited vertical resolution, we are able to extract the precise information about the location, depth and the shape of the geological structure at Valhall within the 600 m depth that are previously highlighted in de Ridder & Dellinger (2011), Mordret et al. (2013a, Mordret et al. (2014a) and de Ridder (2014).
On the profiles BB and CC (Figs 13b and c, respectively), a high-velocity anomaly can be seen in the shallow part of the model overlaying a low-velocity anomaly in the deeper part. This kind of velocity anomaly corresponds to the subsidence at Valhall (Pattillo et al. 1998;Olofsson et al. 2003;Barkved et al. 2005;Zwartjes et al. 2008;Hatchell et al. 2009). It is very well related to the oil production, there is a compaction in the reservoir that leads to the subsidence the whole overburden. The subsidence varies along the rock column and at the reservoir level, it is stronger (∼10 m) than at the surface (∼6 m; Kristiansen & Plischke 2010). This subsidence differential stretches the rocks in the overburden from ∼180 m down to the reservoir and the resulting volumetric strain decreases the seismic velocities (Barkved et al. 2005). While, in the shallower part the subsidence of the seafloor creates a contraction regime that leads to the increase of the seismic velocities (Barkved et al. 2005).
We get significant negative and positive radial anisotropy at Valhall (Fig. 14). The negative radial anisotropy above 220 m depth (Fig. 14a) could be due to the vertically oriented cracks. This type of radial anisotropy is observed in the mantle below mid-ocean ridge (Ekström & Dziewonski 1998;Zhou & Nolet 2006;Nettles & Dziewonski 2008) and rarely observed in the crust (Huang et al. 2010;Xie et al. 2013;Mordret et al. 2015;Tomar et al. 2015;Xie et al. 2015). Production induced seafloor subsidence can produce radial discontinuities that result in anisotropy in top hundreds of metres. This mechanism has been suggested as a possible cause of the azimuthal anisotropy observed in shallow part of the Valhall overburden (Muyzert & Kommedal 2002;Wills et al. 2008;Hatchell et al. 2009;Mordret et al. 2013c).
The prominent positive radial anisotropy (Fig. 14b) found between 220 and 600 m can be due to the finely layered medium (stratification) (Backus 1962;Huang et al. 2010;Erdman et al. 2013;Wang & Montagner 2013;Jaxybulatov et al. 2014;Thomsen & Anderson 2015). During the inversion, the azimuthal anisotropy, which is significant in Valhall (Mordret et al. 2013a) is not corrected. The velocities used to characterize the radial anisotropy are isotropic Scholte and Love wave phase velocities. Azimuthal anisotropy in Valhall is thought to be caused by subvertical cracks in the shallow overburden with a dominant orientation. Negative radial anisotropy, on the other hand, can be caused by subvertical cracks with random orientations. It is possible that the small magnitude of negative radial anisotropy in the shallow part of the model indicates that there are few randomly oriented cracks and that azimuthal anisotropy dominates at this level. However, the deeper positive radial anisotropy attributed to the sediments horizontal layering can hardly be biased by azimuthal anisotropy, therefore we are confident about our deep positive radial anisotropy results, even without azimuthal anisotropy correction.

C O N C L U S I O N S
Only 6.5 hr of continuous recording of noise data at Valhall LoFS network allowed us to image Valhall structure down to 600 m depth and to observe high-resolution 3-D S-wave velocity structure with the depth-dependent radial anisotropy information. A set of 2-D phase-velocity maps of Scholte wave at periods 0.7-1.6 s and of Love wave phase-velocity maps at periods 0.8-1.5 s were constructed to image the near-surface structure of Valhall using the eikonal tomography method. Near-trace stacking helps us to improve SNR for Love wave signal. A Monte Carlo inversion is applied to invert average dispersion curves of Scholte and Love waves to get 1-D isotropic velocity (V S ). Our depth inversion results show good agreement with the geological structure already found using controlled source seismic observations (Sirgue et al. 2010). The ensemble of the data could not fit with an inversion based on an isotropic parametrization. Therefore, we included the depthdependent radial anisotropy into the inversion. The local dispersion curves obtained from phase-velocity maps, are inverted to obtain the 3-D distribution of V SV and V SH in the uppermost crust (down to the depth of ∼600 m). Summarizing these distributions with their mean and standard deviation at each location, we observed a significant negative radial anisotropy (V SV > V SH ) with the average amplitude of −10 per cent ± 1.5 per cent in the shallow layers and the positive radial anisotropy (V SV < V SH ) with the average amplitude of 15 per cent ± 2.5 per cent in the deeper layers at Valhall.

A C K N O W L E D G E M E N T S
The authors gratefully acknowledge the support of ADEME, Schlumberger and Total in this project. We give our sincere thanks BP Nogre AS and Valhall partners Hess Norge AS for granting access the seismic data [and their permission to publish this work]. The work of NS was supported by the Russian Science Foundation (grant no. 14-47-00002).

A P P E N D I X A : H I S T O G R A M O F 1 0 0 0 B E S T M I S F I T
We compute histogram of the 1000 best misfits (Fig. A1) for four cells that are chosen randomly from different position. It can be seen that there are few outliers in these 1000 models but these outliers are also very close to the best misfit and therefore they do not change the final results. Also we show the what misfit span is shown at each point. By misfit span, we mean (max (1000 misfits) -min (1000 misfits))/(max (absolute -misfit) -min (absolute -misfit)). Here, absolute minimum misfit is 0 and absolute maximum misfit is 0.8. So, for example, for the first point (Fig. A1a), the span would be (0.558 − 0.5545)/(0.8 − 0) = 0.44 per cent meaning that the width of the distribution of the 1000 best misfits span only 0.44 per cent of the width of the distribution of all misfits (for the whole map). Similarly, for second point (Fig. A1b), the span is = 0.34 per cent, for third point (Fig. A1c), the span is = 15 per cent and for fourth point (Fig. A1d), the span is = 0.11 per cent. Therefore, the 1000 best misfits are very concentrated around their best value.

A P P E N D I X B : M I S F I T C O M P U TAT I O N
Fig . B1 shows minimum misfit maps that are constructed combining the misfits from all 1-D vertical profiles, computed for the inversion with three parameters (Fig. B1a) and with four parameters (Fig. B1b). It can be seen that the misfit in the inversion is improved using four parameters in comparison of using three. Figure B1. The maps of the minimum misfits from the inversion with the power-law parametrization at each cell; the inversion with (a) three parameters and (b) four parameters.

A P P E N D I X C : F -T E S T
We then check if this improvement in misfit is significant. The Fvalue is distributed with q − p and N − q degree of freedom; here, N = 18 (10 phase-velocity measurement of Scholte waves and 8 phase-velocity measurement of Love waves), first we fit N = 18 points dispersion curves with p = 3 parameters and second, with q = 4 parameters. The probability P f (F, v 1 ,v 2 ), for individual depth profile at all locations is shown in Fig. C1. One can see that the probability is less than 1 per cent for whole grid except some points in the northern part. Therefore, we have 99 per cent chances that the inversion with the four parameters bring more information about the subsurface at Valhall than the inversion with the three parameters.
On the basis of this statistical analysis, it is considered that the two layers of anisotropy significantly improved the inversion results, hence the inversion with four parameters is considered for final model. Figure C1. The p-value map from F-test to compare the two parametrization and check that the improvement in the fit with four parameters is significant or not. The probabilities are around 0 per cent for whole region except some points in north part of the network, it indicates that we need to include the two layers of anisotropy in the inversion, while in the north part of the network the negative anisotropy should be weak.

A P P E N D I X D : M I S F I T C O M P U TAT I O N F O R H Y B R I D I N V E R S I O N
Fig. D1 shows misfit map of the hybrid model, where misfit anomaly due to palaeochannel disappears. Figure D1. The map of minimum misfits for the final model that is constructed using the Gaussian layer parametrization for large palaeochannel region and power-law parametrization elsewhere. It could be noted that the misfit anomaly in the southeast part is disappeared with the modified parametrization.