Integrating streamer and ocean-bottom seismic data for sub-basalt imaging on the Atlantic Margin

ABSTRACT Stacked basalt flows cover much of the NW European continental margin, including potentially prospective sediments of the Faroes shelf. Such flows attenuate seismic energy, hindering sub-basalt structural imaging which is critical for both exploration and tectonic studies. Low-frequency, long-offset reflection surveys have yielded improved images below top basalt, while coincident ocean-bottom seismometer (OBS) data have mapped low velocity zones (LVZ) from tomographic inversion. Image correlation in a common depth domain is challenged by low spatial resolution of the tomographic image, the absence of turning rays in a LVZ and the lack of wide-angle arrivals in the reflection data. We integrate densely sampled reflection data with deep velocities from OBS data to give a common velocity model in a new, iterated, pre-stack depth-migration workflow using complementary constraints from the two datasets. The matched velocity model and depth image enable interpretation of (i) a uniform flood basalt sequence, 2–4 km thick beneath the Fugloy Ridge, with velocities correlating with those in the Lopra 1/1A borehole and (ii) a sub-basalt LVZ with chaotic reflections corresponding to sill-intruded, probably syn-rift sediments. Such a workflow could be extended to targets beneath high velocity salt or basalt and could provide constraints for 3D datasets. Supplementary material: Overlays of modelled arrival times on all the OBS gathers, similar to Figure 4, can be found at http://www.geolsoc.org.uk/SUP18433


INTRODUCTION
The ability to image through basalt flows is a major requirement for studying the structure and hydrocarbon potential of many continental shelves, including, inter alia, those of NW Europe, Brazil and the west coast of India. Here, we focus on techniques for imaging through thick basalt flows using seismic data from the Faroes shelf on the European continental margin (Fig. 1).
Continental break-up of the northern North Atlantic at c. 55 Ma was accompanied by the eruption of huge volumes of basaltic lavas. The basalts flowed across the continental hinterlands on both sides of the new ocean basin. In the North Atlantic region, the total volume of extrusive lavas reached more than 1 million km 3 (White & McKenzie 1989;Coffin & Eldholm 1994;Eldholm & Grue 1994). On the Faroes shelf the basalts are up to 5 km thick and extend across the Mesozoic basement (Sørensen 2003;White et al. 2003). In addition to the extrusion of thick basalt flows, extensive magmatism during 63-54 Ma was expressed through development of igneous centres, such as those on Skye, Mull and the Faroe Islands, and also included the intrusion of sill complexes in the Faroe-Shetland Basin (Naylor et al. 1999;Ritchie et al. 1999).
Flood basalts created at the time of break-up between the Faroe Islands and east Greenland extend over 250 000 km 2 (Andersen 1988;Waagstein 1988;Larsen et al. 1999), at least 40 000 km 2 of which lie in the Faroe-Shetland Basin (Naylor et al. 1999). The flood basalts flowed eastward and southward from their source in the main rift zone west of the present Faroe Islands, decreasing in total thickness from more than 7 km on the Faroe Islands (Richardson et al. 1998(Richardson et al. , 1999, to zero some 150 km away, mid-way across the Faroe-Shetland Basin (Fig. 1).
The main impediment to seismic imaging through basalts is that the layered basalt flows remove the high frequency components of reflections from beneath the basalts. This may be caused by one or more of the following: by scattering from irregular surfaces within the basalts (Maresh et al. 2006); by peg-leg multiples within the basalts caused by the many high impedance contrasts at flow tops and bottoms (Christie et al. 2006b); or by fluid flow within narrow cracks (Shaw et al. 2008). Whatever the physical cause, the effective quality factor, Q, of stacked basalt flows is very low: typically 20-35 in the Faroes region.
The remedy for lack of high frequencies in the signal returned from beneath the basalts is to ensure that the seismic source and the receivers are well tuned for low frequencies.
Typically, a centre frequency of c. 10 Hz will enable good penetration. The penalty of this is that lower frequencies provide reduced resolution, but reduced resolution from a low frequency source is better than no resolution because no energy is returned from beneath the basalts. The data discussed here were recorded with a source designed to produce low frequency energy and with the streamer towed deep (at 18 m below the sea surface) to enhance low frequency recording.
The approach we advocate here to improve the sub-basalt image is first to record good seismic reflection data using a powerful low frequency source and a long (here 12 000 m), deep-towed streamer (Hoare et al. 2005). By using ocean bottom seismometers we then produce a good compressional (P-) wave velocity model for the whole crust using a combination of the moveout of the shallower arrivals on the streamer data to constrain the sediment velocities, with deeper crustal diving waves and wide-angle reflections to constrain the deeper velocity structure through the basalts and underlying rocks (Parkin & White 2008;White et al. 2008;Roberts et al. 2009).
This was expanded to shear (S-) wave analysis, which assists in interpretation of rock types, especially in the region beneath the basalts and the underlying crust of the continent-ocean transition (Eccles et al. 2009(Eccles et al. , 2010. These separate studies of the conventional reflection data and the wide-angle data from ocean-bottom seismometer (OBS) were used as a starting point for further combined analysis, following favourable results from a pilot study of the streamer data (Spitzer et al. 2005). Finally, we show here that an iterated workflow, combining the tomographic velocity control with the streamer data, enables us to tune the seismic processing to remove multiples in the reflection profile data and to improve the pre-stack depthmigration image and interpretation of the structure of the crust beneath the basalt flows.

DATA ACQUISITION
The iSIMM02A seismic profile acquired in summer 2002 extends from the southeastern flank of the Faroe-Shetland Basin, crosses the Basin and the Fugloy Ridge to the east of the Faroe Islands and ends in the Norwegian Sea (Fig. 1). This 375 km profile is orientated (NNW-SSE) approximately in the seafloor spreading direction inferred from the oceanic magnetic anomalies. Multichannel seismic (MCS) and OBS data along the same profile were collected separately in two legs. Full details of the acquisition are found in White et al. (2002).
In the first leg, the NERC vessel RSS Discovery deployed 85 four-component OBSs along the profile (Fig. 1) at a spacing of 6 km, except over the southern slope of the Fugloy Ridge, our focus of interest, where the spacing was reduced to 2 km. The line was shot twice using a large-volume airgun array source (104 l), towed at about 20 m depth and synchronized using two different tuning modes, with the goal of maximizing the low-frequency bandwidth of the seismic signature (Lunnon et al. 2003). Low frequency energy has been shown to improve sub-basalt imaging by minimizing signal attenuation from scattering within the basalt (Ziolkowski et al. 2001;Maresh & White 2005;Maresh et al. 2006;Lau et al. 2007). Synchronization on the first bubble (Avedik et al. 1993) was compared with conventional synchronization on the peak pulse. While both signatures had similar amplitudes and dominant frequencies, the bubble-tuned signature was more compact and exhibited less reverberation (Lunnon et al. 2003). However, because source signatures were not available for every shot and there is more confidence in the shot instant for the peak-tuned data, these OBS data form the basis for the present analysis. The source was fired at approximately 100 m intervals (> 60 s), to avoid water-wave wrap-around, and the data time sampling was 4 ms. To monitor the speed of sound through the water column, Expendable Bathymetric Thermographs (XBTs) were also deployed.
The second leg was conducted by the WesternGeco vessel M/V Geco Topaz, which acquired coincident MCS data using a 12 km, Q-Marine streamer (WesternGeco 2002) along the same line. The streamer comprised 3840 single sensors, spacing at 3.125 m, with a near-offset of 105 m and towed at 18 m water depth. Data from each sensor were digitally combined to a group spacing of 12.5 m using a noise-adaptive filter (Özbek 2000) which preserved the low frequency signal response. Data were recorded up to 18 s at a sample interval of 2 ms. Similarly to the OBS leg, a large, low-frequency airgun array (167 l) was towed at 18 m depth and fired at a 50 m interval. Because near-gun source signatures were recorded for each shot, allowing far-field signature reconstruction and subsequent deconvolution (Ziolkowski et al. 1982), this line was acquired using bubble-tuning (Christie et al. 2006b).  Fig. 1. Bathymetry (GEBCO_08 Grid) and the location of iSIMM profile. This paper focuses on the red segment over the Fugloy Ridge. The small circles are OBSs and data from those in yellow are used in this paper. Black dashed line is the southeastern limit of the area covered by flood basalt (Stoker et al. 1993). Red circle is the Lopra-1/1A borehole. Inset map shows the location of the main map relative to UK.

WORKFLOW: COMBINED ANALYSIS
The combined analysis presented in this paper is a final imaging process that merges and refines results of previous analyses by Spitzer et al. (2005) and Roberts et al. (2009). The combined analysis workflow is essentially a normal pre-stack depthmigration (PreSDM) model-building process that incorporates ray-traced modelling to match the OBS travel-time data (Fig. 2). The purpose is to yield both a velocity model and a final image that are better-constrained by both datasets (Table 1).
The pre-migration processing of the streamer data is identical to that described in Spitzer et al. (2005). In summary this includes: a shot-by-shot deconvolution of data by the bubbletuned source signature; band-pass filtering; frequencywavenumber (F-K) filtering; geometrical spreading correction; multiple modelling and attenuation by wavefield inversion, and parabolic Radon Transform filtering. These steps enhance primary signals and attenuate multiples without the need for an accurate primary velocity model, which can be estimated by pre-stack migration later in the process after the strong, interfering multiples have been attenuated. The same premodelling processing of the OBS data, as described in Roberts et al. (2009), also applies in this analysis. This included offset and time corrections, due to the estimated drift of each OBS descending through the water column from the surface position at deployment, OBS clock drift and occasional time jumps. Peak-tuned vertical component OBS data were used for the modelling. Figure 2a shows a standard industry PreSDM workflow that normally uses only streamer data. This is the embedded workflow of our combined analysis workflow (Fig. 2b) and it is helpful to explain the latter workflow by introducing first the basic workflow and how it was adapted to our work. PreSDM model-building proceeds in a top-down manner.

Embedded basic pre-stack depth migration (PreSDM)
1. The initial velocity model for the first layer is built using a smoothed version of the interval velocity obtained from the best stacking, or root mean square (RMS), velocity model. For our water layer, velocities were constrained by the XBT data. 2. The data were pre-stack depth-migrated using this model to yield common image point (CIP) gathers. 3. Further post-migration multiple attenuation was performed with reference to the migration velocity used. This was important to allow more accurate picking of residual moveout (RMO). The CIP gathers were converted by vertical stretching to the time domain using the migration velocity model. The same model was used to convert the data back to depth. 4. The RMO was then picked on CIP gathers either automatically or by hand. In principle, when the migration velocities approach that of the true medium velocity, the correctlyimaged reflection will be focused at one depth (i.e. the reflection becomes flattened across the gather) and, so, the RMO should become zero. Due to the pervasive presence of peg-leg multiples, picking automatically below top basalt yielded poor results and so the RMO was hand-picked instead. For our processing package, hand-picking was performed by picking the RMS velocity in the time domain. 5. The RMS picks of RMO were then used to generate a set of linear tomographic equations relating changes in the input model velocities to changes in RMO. Subsequently, the changes in velocity model that are required to minimize the RMO were found by inversion using a linear assumption. 6. Such linearized inversion solutions may, however, be nonunique due to ray coverage being restricted to near-offsets.
To avoid over-fitting with unrealistic, short-wavelength structures, the procedure was first to obtain a spatially long wavelength solution and then to fine-tune this solution with progressively shorter-wavelength solutions. Due to the linear assumption, a damping factor was required to avoid large changes (>15%) in velocity in the solution. Therefore, the Core Model Update Cycle (Steps 2-6) was normally repeated several times until the solution converged. Due to real earth complexity, it was impossible to reduce the RMO to zero. 7. Once the migration interval velocity within the current deepest layer was finalized, a new horizon within this layer was picked from the migrated section to be the top of a new, deeper layer. Steps 1-7 were repeated for each successive new layer by holding the shallower layer velocities fixed until no more layers were needed. 8. Since the pre-stack depth migration algorithm was ray-based (Kirchhoff), a smooth model was necessary to minimize non-physical ray scattering. However, the model updating method usually created too much short-wavelength heterogeneity so the resulting model update was smoothed to obtain a final model. 9. This smooth velocity model was then used to migrate the data to produce a final migrated section.
The iterated velocity estimation and imaging sequence described above is normally applied to industry data using sub-critical arrivals to provide both the velocity model from the kinematics and the image from the waveforms. In this dataset, we have the benefit of wide-angle ray-paths from the OBS, which probe the deep subsurface along the same profile, and so we incorporated into our velocity model-building workflow a velocity inversion widely used in the academic community for OBS data, the RayInvr algorithm (Zelt & Smith 1992). This combined analysis approach (Fig. 2b) is explained below with respect to different model layers (Table 1).

Post-basalt (Steps 1-2)
The water velocity was estimated by using XBT measurements for a 1D model. This water velocity was not iterated. The water-bottom boundary was then digitized using the stacked PreSDM streamer data. The initial model for post-basalt sediment velocities was taken from that of Spitzer et al. (2005) and Roberts et al. (2009). These were both derived from interval velocities calculated by WesternGeco (2002; Fig. 3) from the streamer data-stacking velocities. Although updated by the Core Model Update Cycle (Fig. 2a), the final model was very similar to the starting model, showing that the initial stacking velocities were good. The base of this sediment layer, which is also the top-basalt horizon, was then digitized after migration. Note that OBS data do not provide as good a constraint on these two layers as do the streamer data, due to the sparse spatial sampling of OBS compared to that of the streamer (Table 1).
To model the basalt velocity, the RayInvr algorithm was applied to the OBS travel-time data. But, before that, we need to define the post-basalt layers in the RayInvr velocity model. For these layers, we used the same model as in Roberts et al. (2009): a constant velocity of 1.47 km s -1 for the water layer and a gradient layer for the sediment. Note that this model is not exactly the same as the PreSDM iterated model but the two are very similar, as mentioned above, as they were both derived  for fitting the sea-bottom and uppermost basalt arrivals. Since the two models should have similar overburden travel times for deeper arrivals, we decided we could safely omit the tedious step of adapting the PreSDM iterated model to RayInvr format.

Basalt (Steps 3-4 & 9)
We forward-modelled the basalt layer in RayInvr by fitting the basalt refracted phase (P B ; Figs 4-6a; Table 1). This phase was observed clearly as a first arrival in all the OBSs and provided excellent constraints on the basalt velocities throughout the model. In order to provide a good fit to the data, two layers with different vertical velocity gradients were required to model the arrivals from the basalt sequence. Velocity was continuous between these two layers as there was no evidence of a discontinuity in the data. The same P B phase could also be observed in the raw streamer data (Fig. 6c). In fact, Figure 6a shows that coincident streamer and OBS data are very similar as the OBS data merge almost seamlessly with the streamer data in the overlapping region. This formed the basis of a tight integration of both datasets. Unfortunately, the limited length of streamer (12 km) enabled refracted arrivals to be observed from only the topmost 1 km of this layer on the streamer data.
Using conventional semblance processing, the reflected streamer data also provided weak constraints on basalt velocities (Table 1). Much of the larger offsets were muted due to the wavelet stretch associated with normal-moveout correction (NMO) and pre-Radon processing (see section on Multiple removal). Figure 7 shows an example CIP gather which was used to pick RMO. Within the relatively shallow basalt layer, the streamer data were muted to 1 km offset at the top and 4 km at the bottom. Semblances were poorly defined within the basalt sequence and only slightly better at the base. Therefore, RMO picks from the streamer data within the basalt layer were guided heavily by the velocities obtained from RayInvr modelling in Step 9 (black curve; Fig. 7) and the PreSDM was mainly helpful in updating the migrated image of the basalt. Note that although the RMS picks in Figure 7 reflect residual moveout, they also give good insight into the interval velocities needed for focusing the reflections in the CIP domain.  To model the depth of the base basalt horizon, the boundary was constrained by both streamer and OBS data according to the following prioritized criteria: (i) PreSDM image. The base basalt horizon, being a major lithostratigraphic boundary, is likely to correlate with a change in reflection pattern in the PreSDM section ( Fig. 8).
Step-back in the time-offset domain. This is observed at many OBS and is strong evidence for the existence of a velocity reduction beneath a higher-velocity basalt layer, but   The corresponding ray-path diagram with the same colour code as above. Only top basalt and deeper phases are shown for clarity. The seismic section was pre-stack depth migrated using the same final P-wave velocity model (colour map) as is used in the ray-tracing. See text for modelling details. Similar plots for the OBS analysed can be found in SUP18433. the depth of the velocity reduction is not tightly determined due to uncertainty in picking the precise offset at which the travel-time step-back occurs (Fig. 4). The step-back is the result of a low velocity layer beneath the basalt layer . Rays entering the LVZ bend toward the normal (i.e. become more steeply dipping) and are not refracted back to the surface in the low velocity layer. Instead, they continue to dive until they penetrate a layer whose velocity exceeds that of the basalt, assumed to be the basement. This produces a step-back in time between the basalt and the basement-refracted phases (Fig. 4). (iii) Base basalt wide-angle reflection. Ideally, the base-basalt reflection time could be picked by correlating the wideangle reflection from the step-back to zero-offset. In reality, because of strong multiples at short offsets, heavy processing was required to bring out this reflection and hence adds uncertainty to its observed travel-time ( Fig. 6b; see also section on Multiple removal). Also, the presence of a wide-angle reflection does not always correspond to a velocity layer boundary. Where it could be picked, the wide-angle reflection provided good constraint on the approximate depth of the velocity reversal. For extra constraint, both OBS common receiver gathers and their corresponding shot gathers from raw streamer data were also used to model the wide-angle reflection (Fig. 6c).

Low-velocity zone (LVZ; Steps 5-11)
As mentioned above, there is no refracted wave from within the LVZ to constrain its velocity on long-offset OBS data. Furthermore, it was also difficult to estimate short-offset moveout of the reflections within this layer due to sub-basalt peg-leg multiples. Unlike the streamer data, multiple attenuation of short-offset OBS data usually yields poor results due to sparse trace spacing (100 m). Therefore, we assumed a velocity of 4.5 km s -1 (Roberts et al. 2009) as the best estimate for a starting model ( Fig. 2b; Step 5) to be updated later by the streamer side of the workflow. We then completed the RayInvr model by adding the deeper layers. Since the base of the LVZ was difficult to identify confidently from the OBS data at short offsets, we modelled the turning wave through the layer beneath, assumed to be basement (Figs 4 & 5; Step 6). The initial basement velocities were taken mainly from the model of Roberts et al. (2009). Reflections observed in the streamer data provided a better constraint on the velocities within the LVZ (Fig. 7). However, the RayInvr velocity model built in Steps 5-6 is not a suitable format for PreSDM. Therefore, we updated the PreSDM velocity model resulting from Steps 1-2 by merging it with the RayInvr model velocities for the basalt and sub-basalt layers (Steps 7-8). For simplicity, these layers became a single gridded layer in the PreSDM velocity model. Figure 7b shows a series of reflections with high semblance values that well define the LVZ. Unlike reflections within the shallower basalt layer, these reflections were observable over a wider offset distance range of traces, adding further confidence. They also helped to constrain the base LVZ boundary which could be defined by an increase in interval velocity compared to that within the overlying LVZ. By going through the Core Model Update Cycle again (Step 9), migration velocities were updated to produce a better focused depth image. These model updates were most significant for the LVZ layer as the basalt velocity updates were similar to the starting model and semblance picks for sub-LVZ reflections were less well defined (Fig. 7). At the end of the workflow, a final PreSDM velocity model was obtained. Figure  7c shows the same CIP gather converted to the time domain after depth migration using the final PreSDM model. Subsequently, velocities from the final PreSDM model were tested and adjusted by the RayInvr modelling process using OBS data (Steps 10-11). While their integration is conceptually straightforward, the PreSDM and RayInvr models define velocities differently. The PreSDM model comprised different layers, bounded by digitized horizons and, within each layer, the velocity field was defined by a rectangular grid. On the other hand, the RayInvr model, although also layer-based, did not use a grid within each layer but a linear velocity gradient defined by point velocity values at both top and bottom bounding horizons. While it was quite straightforward to transfer PreSDM boundaries to the RayInvr model format, we also needed to estimate linear profiles from the gridded velocity.
To do this, velocity profiles were obtained from the PreSDM model at every 10 km distance along the profile. Linear velocity gradients were estimated to update the lower basalt and the LVZ RayInvr model layers. Base-basalt and base-LVZ depths in the RayInvr model were adjusted to match the depths of the strongest velocity changes from the final PreSDM model. This model can be viewed as a smoothed version of the final PreSDM model. The LVZ velocities were unchanged but the base-LVZ boundary depths were forward modelled to fit the timing of the basement phase after the step-back (Fig. 4). Figure 5 shows data from OBS 58 which constrain the pinch-out of the LVZ in the north-northwest, based on the absence of a step-back.

Basement (Steps 11-13)
In a similar manner, the velocity for the upper basement layer beneath the LVZ was fine-tuned for a better fit with OBS travel-time data and reflectivity on the migrated MCS section. The modelling of deeper crustal layers is detailed in Roberts et al. (2009) and is beyond the scope of this paper. The resulting velocity model is the final RayInvr model (Fig. 8). This model was then used to replace the deepest layer of the PreSDM velocity model which, in turn, produced the final PreSDM imaging of the streamer data over the Fugloy Ridge ( Fig. 8b; Step 13).

Multiple removal
The post-migration multiple attenuation procedure (Fig. 2a) on the streamer data was crucial in revealing primary events in the CIP gathers before RMO values can be picked with confident. The migrated gathers were bandpass-filtered between 5 Hz and 30 Hz in the time domain, followed by high-resolution parabolic weighted-least-squares Radon filtering (Moore & Kostov 2002) and FK-filtering. The Radon filtering discriminates against multiples on the basis of their slowness and intercept time in the parabolic Radon domain, while the FK-filtering discriminates purely on the basis of apparent velocity in the frequency-wavenumber domain. Pre-water-bottom arrivals and signals overstretched by NMO correction at large offsets (so-called 'hockey sticks') were muted before Radon filtering. Due to the limited data extent in offset after muting, we used a 20% overlapping time window of 2000 ms and a 50% overlapping offset window of 2000 m to model coherent events. Filters were designed in both Radon and FK domains to reject modelled multiple events. A more forgiving filter was initially used to avoid rejecting primary events lying close to multiples in the Radon domain due to inaccurate migration velocities, but the filter was eventually tightened to retain only flattened to near-flattened events, as confidence in the velocity model improved, for stacking into the final migrated section.
To constrain the depth of the base of the basalts, an atypical, gather-by-gather multiple attenuation approach was used which enhanced base-basalt wide-angle reflection on the OBS data (Fig. 6). This involved bandpass filtering (5-30 Hz), predictive deconvolution, block moveout correction (Rupert & Chun 1975) to distinguish this wide-angle reflection from other reflections based on moveout differences, and FK-filtering. The block moveout correction, which lined up targeted reflections by time-shifting the whole trace by the required moveout, avoided stretch muting in the far offsets. The subsequent FK filter, however, removed all shallower reflections whose velocities were lower than that of the base-basalt wide-angle reflection and each gather had to be moveout-corrected individually to enhance the wide-angle reflection. While this approach was unsuitable for PreSDM analysis, which involves many more gathers than OBS analysis, it provided additional constraints for RayInvr velocity modelling (Table 1, Fig. 2b). Figure 8a shows the final P-wave velocity model resulting from the combined analysis as described above. The velocity structures delineate post-basalt sediment, basalt layers, a sub-basalt low velocity layer and an underlying high velocity basement layer. The post-basalt sediment increases in thickness from 70 m over the shallowest part of the Fugloy Ridge to 1.5 km towards the flanks of the ridge. The sediment increases in velocity with depth from 1.6 to 2.6 km s -1 . Beneath the sediment, the seismic velocity increases markedly to 3.5 km s -1 Under the high velocity-gradient layer in the upper section of the basalts, we observe differences in the velocity structures between the basalts overlaying the continental zone and those of the transitional zone (Fig. 8a). In the continental zone, the basalt velocities range from 5.2-5.7 km s -1 with very low vertical velocity-gradients (0.12-0.20 s 1 ). The total basalt thickness is 2-4 km over the continental section.

Velocity model
A clear sub-basalt low velocity layer of 4.5-5.5 km s -1 is delineated by the modelling. This layer is thickest (c. 2 km) beneath the axis of the ridge and thins toward the ridge flanks pinching out at the seaward end of the continental zone (at distance 157 km). The average velocity of the sub-basalt layer increases as the thickness of the layer decreases. The velocitygradient within the sub-basalt layer is very low (0.12 s -1 ) in the central region (186-212 km), but increases towards the two ends.
Underneath the low velocity zone, the velocity of the basement increases to 5.6 km s -1 at c. 6 km depth and 6.2 km s -1 at c. 9 km depth with a velocity-gradient similar to that of the lower basalt sequence. In the transition zone, the velocity-gradient at the same depth increases considerably to 0.49 s -1 beyond the pinch-out of the low velocity zone. However, a low gradient but high velocity layer (6.5-6.8 km s -1 ) is observed beneath 6-8 km depth, in contrast with the continental zone at the same depths. A more detailed discussion on this transitional layer, interpreted as intermixed continental crust and igneous intrusions, can be found in Roberts et al. (2009). Figure 8b shows the final pre-stack depth-migrated section of the streamer data using the velocity model shown in Figure 8a. When plotting both together (Fig. 8c), it becomes clear that the reflection pattern correlates closely with trends in the velocity model. The reflection from the top of the basalt is identical to the boundary identified in the velocity model which marks the top of the high velocity-gradient basalt layer. The strong vertical change in velocity-gradient into the main basalt layer is not manifested by any noticeable changes in the reflectivity but individual reflection events become sparser with depth.

Reflection pattern
Consistent sub-horizontal layering within the basalt sequences is observed. This is in contrast with the more chaotic pattern within the LVZ where some saucer-like structures can be observed (e.g. 182 km distance at 3.5 km depth; 199 km distance at 4.5 km depth). In the LVZ, the reflection amplitude is also higher with some reflections potentially defining the base-LVZ boundary. Underneath the LVZ, the reflection pattern changes markedly to less coherent reflections with relatively lower amplitude, except for a band of high amplitude coherent reflections at 8.8-9.0 km depth and 182-208 km distance. A deeper crustal reflection at 9.5 km depth is also observed near 160 km distance.
We also observed changes in the reflection pattern towards the transition zone. At 2-6 km depth, the basalt layers are characterized by a series of upward convex seaward-dipping reflectors (SDRs). SDRs are widely observed at the continentocean transition of volcanic margins (Menzies et al. 2002). They represent lava that flowed landward from fissure vents at, or near, the incipient spreading axis but subsequently tilted seaward due to thermal subsidence and loading as the new spreading axis receded, giving rise to their present shape. Underneath the SDRs, the reflectivity become very low and is coincident with the high velocity transitional crustal layer.

Constraints
The constraints on the final RayInvr velocity model by the OBS and of the final PreSDM model by the streamer data are shown in Figure 9. For the RayInvr model, Figures 9a-c show the ray coverage of the iSIMM OBSs by ray-tracing different traveltime phases through the final model. Within the basalt layer (Fig. 9a), both the refracted and reflected phases provide complete coverage of the layer with landward and seaward travelling rays crossing each other (except for beyond the most seaward OBS). Crossing ray-paths mean that both the depth and the velocity of the layer are well resolved. In contrast, the velocity of the LVZ is poorly constrained by the OBS data (Fig.  9b), because no refracted phases are returned to the surface from within this layer. The velocity within the LVZ is constrained only by reflections from the bottom of the LVZ. There is limited ray coverage at the seaward end of the LVZ. For the basement layer beneath the LVZ, there is very good ray coverage from the data recorded by the OBSs with turning rays from both directions crossing each other (Fig. 9c). This again provides good velocity and depth control in the tomographic inversion. The seaward end of the transitional crust is, however, less well covered because it is near the edge of the model. There are also gaps near 205 km distance at 6.5 km depth. The area without ray coverage at the base of the model corresponds to the lower crustal layer of Roberts et al. (2009) and is beyond the scope of this paper. Figure 9d is an illumination diagram showing the constraints of the PreSDM model velocities by residual moveout (RMO) picks using streamer data. The basalt layer is represented by a high number of RMO picks which constrain that part of the model. We improved confidence in picking the basalt layer RMO after the basalt velocities had been constrained by OBS data. The LVZ also has a high number of hit counts/cell (200-600) from the streamer data, which means that velocities here are well constrained. This is in contrast with the much poorer level of constraint from the OBS data alone (Fig. 9b). The deeper parts of the model are less well illuminated because reflections are relatively less easily observed at these depths except for a few high amplitude events.

Time section
When interpreting a post-stack migrated section, seismologists routinely compare it with its stacked section counterpart for any inconsistency which may have been introduced by the migration process. It is also helpful to examine the migrated section in the time domain to pick out any multiple events that may have been mistaken as primaries. Since our seismic section is a PreSDM section, we prepared our stacked-unmigrated time section (Fig. 10a) in a separate workflow using pre-migrated pre-stack data. To facilitate a fair comparison between the two sections, the pre-stack unmigrated data had the same multiple suppression processing applied before stack as the migrated section (converted to two-way time; Fig. 10b). They are also plotted to show similar gain within the LVZ.
It can be seen that the chaotic nature and the elevated amplitude of the reflections within the LVZ are already apparent in the stacked time section. These reflections are better defined after steeply dipping diffractions have been suppressed by migration. The difference between these reflections and the high-amplitude top-basalt reflection suggest that they are not likely to be either sea-bottom or peg-leg multiples. Note, however, that the difference in reflectivity between the LVZ and the basement underneath is not as obvious in the time domain as it appears in depth (Fig. 10c) and may result from wavelet stretching lowering the apparent frequency in depth.

DISCUSSION
The discrimination between the basalt sequence and the underlying geology remains the most challenging and yet critical problem of the region. Although sub-top-basalt reflection images have improved greatly due to recent advancement in MCS data acquisition (e.g. Hill et al. 2006) and processing (e.g. Gallagher & Dromgoole 2007), interpretation of base basalt remains a challenge in many areas. As basalt layering can easily be mistaken as sediment layering from the reflectivity alone, our results show that velocity modelling consistent with the reflection image plays an important role in resolving this problem.
Basalt Figure 11 shows the velocity profiles across our final model gathered into three different regions. Our modelled velocities within the interpreted sub-aerial basalt flows (Fig. 12) show striking similarity throughout the model with those of the Lopra-1/1A borehole sonic log (0.6-2.8 km below top basalt), even though the latter shows high resolution variations and two high velocity peaks near the top due to dolerite intrusions (Christie et al. 2006a). Note that our modelling technique and the limitation in OBS coverage have smoothed out these high resolution features seen in well logs. Therefore, our basalt velocities show good agreement with those from the Lopra borehole, suggesting that they have a similar origin. According to the PreSDM section (Fig. 12), these basalt layers are observed to be relatively smooth and continuous for long distances, which strongly suggests sub-aerial emplacement, as at Lopra. Therefore, based on both velocities and the reflection pattern, we conclude that the Fugloy Ridge basalt we imaged is similar to the stacked flow units drilled at Lopra. The high vertical velocity-gradient, low velocity layer in the upper part of the extruded basalt sequence is likely to be caused by weathering due to exposure of the final flooding events and reduction in flow thickness when the Iceland plume moved away as seafloor spreading proceeded. Velocities are lower at the top of each basalt individual flow (2.5-5.5 km s -1 ) compared with their centres and bottoms (5-6 km s -1 ) in the SE Greenland Margin (Planke & Cambray 1998). Therefore, thinner flows tend to decrease the velocity averaged over the length scale detectable by seismic waves. The older and much thicker lower unit of the basalt seismic sequence (Fig. 11) has a low vertical velocity-gradient but high velocities. It is likely to comprise thick, massive flow units that characterize flows from frequent basalt flooding near to the source. Dolerite intrusion may possibly be widespread, as in Lopra. Such similarity may provide the basis for extrapolation in basalt velocities over the area between Lopra and the iSIMM profile, i.e. the eastern Faroe platform (Fig. 1).
As described in the 'Method' section, the constraint on basalt velocity by the streamer data alone is poor (Table 1) due to a combination of high basalt velocities at shallow depths, as under the ridge. Therefore, it is likely that constraints are required from other sources for proper depth imaging of the basalt layer in similar environments. Furthermore, accurate sub-basalt velocities from travel-time inversion are possible only if the basalt velocities are also known well. Therefore, even though the basalt velocities themselves are not of primary interest to petroleum exploration, they are essential for imaging and interpreting the sub-basalt structures. In our case where a LVZ is present beneath the basalt, the OBS data are particularly useful in identifying the base basalt horizon from the measured step-backs in diving wave arrivals. Unfortunately, coincident OBS data are not normally acquired along industry MCS profiles. Furthermore, well-ties are sparse in basalt terranes such as the Faroes shelf and are often a long distance from the seismic profiles. However, extrapolations may be appropriate if basalt flows are shown to be uniform over a long distance, as argued above.

Base basalt
An accurate interpretation of depth of the base of the basalt flows is critical for reducing risk in exploring sub-basalt reservoirs. Our results show that to reach such targets requires drilling through over 2 km of basalt within our area of interest. According to its velocity structure, the basalt sequence to be drilled is likely to comprise massive flows with very little interbedded sediments. This would provide useful input into project planning.
Our base basalt boundary is simply the top boundary of the LVZ, which suggests that a different lithological unit is encountered. Such velocity inversion conveniently provides constraints on modelling the boundary through the use of step-backs and base-basalt wide-angle reflections in the OBS data. Since step-backs are often observed at offsets beyond the normal length of streamers, we infer that top-LVZ interpretation from streamer data alone may be prone to errors. We find that PreSDM velocity modelling using sub-critical reflections from near the base basalt did not provide enough constraint in its identification because of uncertainty in deriving interval velocities from the moveout velocity. Fortunately, the fact that the observed changes in reflection patterns in the PreSDM section correlate so well with our modelled base basalt provided good support for its interpreted depth. However, had the reflection pattern of the sub-basalt unit not been distinguishable from that of the basalt sequence, confirmation could come only from step-backs on the OBS data. This highlights the importance of the use of step-backs and wide-angle reflections from OBS data on identifying and imaging of base basalt.

Low Velocity Zone (LVZ)
Our results show good correspondence between the velocities and the reflection fabric within the LVZ. Here, we propose slightly different interpretations for the three sections of this sub-basalt unit (Figs 11, 12). The main component of the LVZ everywhere is sediments, but the velocities and reflection fabric differ along the LVZ according to the proportion of hyaloclastite and intruded sills.
We interpret the middle section (188-210 km) to be formed by syn-rift sediment up to 2.5 km thick, intruded by sills. Velocities within the LVZ (4.6-4.8 km s -1 ) are lower than those of the hyaloclastites in the lower section of the Lopra well (Fig. 11b). Although the LVZ velocities may be slightly too high for normal sediment, the presence of sill intrusions may be responsible for the elevated velocity in the sediment. This is also supported by the reflections that are similar to the saucer-shape reflections commonly observed from sills within sediments. The chaotic reflectivity pattern of this unit is consistent with perturbation by intrusion, mass-wasting events or deposition under high energy environments. Therefore, this unit may potentially contain a reservoir. Without a well-tie, it is unclear what formation this is equivalent to on the Shetland side of the Faroe-Shetland Trough, although it must be of lower Paleocene age or older.
For the section (212-220 km) where the LVZ thins landward to c. 1 km toward the Faroe-Shetland Trough, the velocity profiles closely resemble those from Lopra (Fig. 11c). Therefore, we infer that the LVZ here is perhaps hyaloclastite prograding to deeper water to the right of Figure 12 into the Faroe-Shetland Trough. Such an interpretation implies a palaeo-coastline near the point where the top LVZ starts to dip steeply (c. 207 km). The palaeotopography on the left (north) of the figure would have been higher than on the right (south), but subsequently rotated seaward due to withdrawal of the spreading axis, with subsequent cooling and subsidence.
The part of the LVZ (160-186 km) where it pinches out in a seaward direction toward the transition zone has velocities intermediate between those of the middle section and the landward section (Fig. 11a). We also find that the reflection pattern is more similar to the flood basalt above (Fig. 8). Therefore, we interpret this region as the sediment-to-basalt transition with the intensity of basaltic intrusion increasing seaward. As this part of the profile was close to the location of the continental break-up, we expect the syn-rift sediment to have the most interaction with the basalt, resulting in a higher percentage of igneous intrusion than the landward section.
Despite these differences between the three sections of the LVZ, the velocities in general concentrate more toward the lower end of the spectrum (Fig. 11d). From this we infer that the LVZ is mainly sedimentary in origin.

Basement
Underneath the LVZ, we modelled high velocities (c. 5.6 km s -1 ) which we interpret as caused by continental basement. Although the velocity model suggests a near-uniform unit from top of the basement to the bottom of the model, we see hints of layering within the upper basement in the migrated section (Fig. 8). Furthermore, there are strong reflections at c. 8 km depth that resemble fault-block structures (Fig. 12). Therefore, we interpret the upper part of the basement as syn-rift basement, which may be highly consolidated early syn-rift sediment, and the layer beneath it as pre-rift basement, which may be Lewisian gneiss (Earle et al. 1989;Dean et al. 1999).
Since the velocities of the interpreted syn-rift basement are high, the syn-rift basement is likely to be composed of clastics with a high degree of metamorphism and possibly uplift and subsequent erosion. The high velocity may also imply a high density of sill intrusions that would increase the sediment velocity. However, unlike within the LVZ, we do not observe individual bright reflections within this layer such as would be produced by sills. It may be that much of the seismic energy has been attenuated within the basalt before being reflected by the sills, if they exist, or that the higher velocity and density of the compacted syn-rift sediments produce less of an impedance contrast with the putative igneous sills, thus making them less visible.
A reflection at 8.5-6.5 km depth and 175-220 km distance is interpreted as the top of pre-rift basement (Fig. 12). Its geometry suggests a normal fault which may be related to the formation of the Fugloy Ridge continental fault block and the crustal thinning toward the Faroe-Shetland Trough. However, a velocity discontinuity is not found across this boundary in our model (Fig. 8), suggesting that the velocity contrast which has caused a reflection from this boundary may be from a layer thinner than is resolvable by OBS travel-time inversion. As this boundary appears continuous with a group of strong, deep reflections (160-170 km distance; > 9 km depth), it may represent another sill intrusion. The coherency of this reflection may, however, have been degraded by over-migrated residual multiples. Below this boundary, the basement is relatively unreflective and the reflection pattern is typical of crystalline continental crust. With velocities of c. 6 km s -1 and a low vertical velocity-gradient, the basement is reasonably interpreted as of continental origin (Fig. 8).
Transition zone Seaward of model distance 156 km, we observe below top basalt an abrupt change in both the velocity structure and the reflection pattern from that of the continental section. This is in agreement with previous results (Eccles et al. 2009;Roberts et al. 2009). Such crust, which can neither be classified as continental nor oceanic, is termed transitional here. The reflection pattern within this crust also correlates well with the velocity model (Fig. 8c). SDRs are observed in the transitional upper crust where the velocities are lower than those of the continental basalt layer. The SDRs suggest thinner flows when compared with the continental section and may explain the relatively lower basalt velocities. SDRs are a result of a progressively subsiding surface as crustal thinning proceeds at the margin and are commonly observed at volcanic rifted margins. A strong vertical velocity-gradient results in high velocities in the lower transitional crust where reflections are hardly observable. However, strong reflections are found deeper in the lower crust (> 10 km) and seaward of the current model where high velocity lower crust is also observed (WesternGeco 2002; Roberts et al. 2009).

Sub-basalt imaging
Although sub-basalt imaging, beneath the Faroe-Shetland Trough in particular, is challenging for reasons mentioned earlier, the potential prize justifies the effort. The result is development of insights into best practice for this geological setting in data acquisition, processing and model-based interpretation. Our combined-analysis depth-imaging approach has allowed us to make a confident determination of the location of the base of the basalts, and has revealed new detail below them and is, therefore, effective. However, the effectiveness of the approach described here depends on two main factors: (1) two-way transmission of seismic energy through the basalt layer with sufficient signal-to-noise ratio to reveal structure; (2) acquisition of coincident streamer and wide-angle OBS data. We have also assumed, with later justification, that vertical-horizontal anisotropic effects are small, although they may be included in the workflow provided the anisotropic parameters can be estimated.
For basalt penetration, the results demonstrate the benefit of using acquisition parameters designed to enhance lowfrequency generation and recording for both the OBS and MCS surveys. Not only are turning wave arrivals through and from beneath the basalt layer more easily observable with low frequencies, but the image of the LVZ using reflected data from the streamer also achieved a resolution finer than has ever been achieved before over the Fugloy Ridge (< 1 km horizontally; < 0.3 km vertically; Fig. 8b). It is the streamer data that provide constraints on the velocities of the LVZ. Therefore, we achieved our goals of both deep penetration and reasonable resolution.
Regarding the need for coincident streamer and OBS data, it is obviously important to ensure that both datasets sample the same transect, otherwise errors may be introduced by projecting velocity structures and images from two different profiles. Furthermore, the application of mutual constraints between OBS and MCS datasets would not be valid, rendering featureto-feature correlation from one dataset to another inappropriate. However, acquisition of coincident streamer and OBS data is not yet common practice, especially in the commercial sector, despite the increased commercial availability of autonomous and cable-mounted seabed systems (e.g. Beaudoin & Ross 2008;Foster et al. 2008). However, our approach, resulting from an academic-industry collaboration, may pave the way for a fresh approach on how seismic data are collected in basalt terranes.
When using both sub-critical and wide-angle data together as constraints, we need to consider anisotropy. Although basalt layers are thought to be potentially anisotropic, borehole seismic results from the Faroe Islands at both the Lopra-1/1A well, and the Vestmanna and Glyvursnes wells, found both vertical polar and azimuthal anisotropy to be very low (Christie et al. 2006a;Bais et al. 2009). Therefore, for simplicity, we assumed any anisotropy on our profile to be negligible.
In addition to the generic difficulty of sub-basalt imaging discussed above, MCS imaging beneath the Fugloy Ridge presents specific challenges: (1) relatively shallow-water depths produce pervasive water column multiples; (2) high velocities just below seabed contribute to many orders of multiples and cause large NMO stretch, limiting the outer mute offset; (3) large impedance contrasts within the top basalt create high amplitude peg-leg multiples. These challenges prompted a significant effort on multiple attenuation, which has been at least partly successful, although leaving some residual multiple energy. Since both multiple attenuation (including stacking) and migration contribute to the observed clarity in sub-basalt structures, it would be interesting to separate their effects to investigate which of them is more important. As we are able to observe detailed sub-basalt structure even in the pre-migration stacked section (Fig. 10a), the improvement in the image quality of the sub-basalt layer is determined more primarily by a better multiple attenuation whereas the migration algorithm provides better focusing of the primary events. Since our approach in multiple attenuation is rather conventional, the key element was the new velocity model produced by our combined workflow which was used to discriminate primary events from the multiples. We, therefore, support the conclusion by Gallagher & Dromgoole (2007) that a critical element of imaging in basalt terranes is iterated velocity analysis which takes the geological model into account.
Our combined workflow can be adapted to other areas which may benefit from its flexibility in using both OBS and streamer data interchangeably as constraints. It should work equally well for imaging below tabular salt bodies as with sub-basalt imaging, where velocity inversion beneath a high velocity layer is a common problem. The OBS data would provide a good constraint on the velocity and thickness of a salt layer, while the streamer data would give better constraint on sub-salt velocities. Even for structures that are without a high velocity layer and velocity inversion, our combined workflow promotes velocity models that are consistent with both OBS and streamer datasets. We were confined to layer-based models by the model representation of the RayInvr modelling approach adopted here, but we encourage future development of grid based velocity models which may be better for more complicated structures. Likewise, this technique should not be limited to just 2D datasets. Where 3D MCS and OBS data are available, the construction of a final 3D velocity model and pre-stack depth-migration results should be possible by adapting our workflow accordingly.

A new workflow with mutual constraints between OBS and
MCS data is reported. The advantage of such an approach over conventional processing techniques is that it is able to produce a subsurface image and velocity model that are consistent with each other and with both datasets. 2. Pre-stack depth migration (PreSDM) with a constrained velocity model improves the primary image. However, the improvement in the quality of the sub-basalt image results primarily from the better multiple attenuation made possible with an iteratively constrained velocity model. 3. Low frequency sources, such as the airgun source we designed (with peak amplitude at c. 9 Hz), help to image below top basalt by providing broad bandwidth at low frequencies to combat scattering and attenuation losses. 4. Our results show excellent correlation of flood basalt flow velocities with Lopra borehole data (Fig. 11). We suggest, therefore, that the basalt flows are laterally uniform in morphology across much of the eastern Faroe platform. 5. Accurate interpretations of the base basalt (or top LVZ) boundary, which is of key interest to hydrocarbon exploration, relies heavily on OBS constrains from observations of step-backs. Based on this, the total basalt thickness was mapped and found to be 2-4 km over the continental zone beneath the Fugloy Ridge. 6. The comparison of LVZ velocities with those of the Lopra borehole guides an interpretation of the LVZ and its lateral variation (Figs 11,12). It is suggested that the LVZ comprises mostly sediment with sill intrusions but becomes progressively dominated by intrusions where it pinches out toward the transitional crust on the continent-ocean boundary. Its landward end towards the Faroe-Shetland Trough is interpreted as part of the prograding hyaloclastite succession. 7. The same workflow could be applied to tabular salt bodies and 3D datasets.