Reprocessing 2-D airgun seismic reflection data SALTFLU (salt deformation and sub-salt fluid circulation in the Algero-Balearic abyssal plain) in the Balearic promontory and the Algerian basin

In an ever more challenging context for the acquisition of seismic data in the Mediterranean Sea, reprocessing to improve the quality of legacy data has become increasingly important. This work presents the newly reprocessed, open access dataset SALTFLU acquired in the Algerian basin by the National Institute of Oceanography and Applied Geophysics (OGS) in 2012. We apply a ‘broadband’ reprocessing strategy adapted for offset-limited (3 km streamer for a target 4 km below the sea level) airgun reflection seismic data acquired in deep water settings. We then assess if the reprocessed images provide new geological insights on the Mediterranean sub-surface. The workflow relies on an integrated approach combining geophysics and geological interpretation to iteratively build the velocity model. In this way we aim to tackle some of the challenges linked to imaging deep complex geological structures containing high velocity contrasts with 2-D, offset-limited seismic data. We first broaden the bandwidth of the data through multi-domain de-noising, deghosting and a source designature using an operator derived from the seabed reflection. We then perform iterative migration velocity analysis, pre-stack time migration and multiple attenuation in the Radon domain to obtain time-migrated images. The initial velocity model is derived from the resulting time migration velocities, and geologically driven model updates are generated using a combination of travel-time tomography, seismic interpretation of the major salt horizons and velocity gradient flooding. The gradient flooding aims to reproduce the large scale first-order velocity variations, while the travel-time tomography aims to resolve the smaller second-order velocity variations. The results improve our deep geological knowledge of the under-explored Algerian basin down to the base salt and the pre-salt. Fluid indicators are imaged within the Plio-Quaternary of the Algerian basin, which we interpret as thermogenic or biogenic gas sourced from either the Messinian Upper Unit or from the pre-salt, migrating through a hydro-fractured salt. The reprocessed data image lateral and vertical seismic facies variation within the Messinian units that could shed new light on the tectono-stratigraphic processes acting during the Messinian Salinity Crisis. It also reveals numerous previously unresolved volcanic structures within the Formentera basin.


Introduction
The Mediterranean Salt Giant (MSG) is a vast and thick evaporitic unit that was deposited during the Messinian Salinity Crisis (MSC), an extreme environmental episode that occurred from ~ 5.97 to ~ 5.33 Ma Krijgsman et al. 1999;CIESM 2008;Manzi et al. 2013). Since its discovery in 1970, and despite intense multidisciplinary research, the MSG and the related MSC are still poorly understood and subject to many unresolved controversies (Camerlenghi and Aloisi 2019). This is in large part due to a lack of data in the deep-water offshore domain, which represents more than 80% of the volume of 1 3 13 Page 2 of 28 MSC-related sediments, emphasizing the need for new and improved offshore data to better constrain the MSC (Roveri et al. 2014). Multi-channel seismic reflection profiling is the most common geophysical method applied to image the architecture of offshore basins and to prospect for potential drilling sites. Acquiring new marine seismic data, however, is challenging due to high acquisition costs associated with marine operations. Numerous offshore seismic datasets have been acquired in the last decades in the Mediterranean Sea, providing countless 2-D profiles that could no longer be acquired today because they cover areas that are currently subject to restrictions related to obtaining exploration permits (Diviacco et al. 2015). Many of these datasets are currently poorly exploited due to lack of public data access or poor-quality seismic processing. Reprocessing legacy data, therefore, is a potential source of new geological information that could be extracted from these dormant datasets, directly contributing to a better understanding of the MSG and the MSC.
The SALTFLU ('Salt deformation and sub-salt fluid circulation') 2-D multi-channel seismic reflection dataset was acquired in June-July 2012 by the R/V OGS-Explora, Eurofleets Cruise No. E12 (acquisition parameters listed in Table 1). The survey was planned to study the influence of the MSG on pore fluid circulation during basin evolution since the post-Messinian. The legacy SALTFLU processing followed a 'narrowband' approach, without deghosting, using narrow bandpass filters coupled with a source designature based on statistical deconvolution and no zero-phasing of the target wavelet. The filtering eliminated much of the low frequency signal (below around 6 Hz), whilst the source designature boosted the high frequency noise and produced a mixed-phase wavelet that has inconsistent phase across the survey. The wavelet also contains strong residual energy (likely from the bubble pulse) that overprints and obscures the primary signal, particularly the shallow geology close to the water bottom (Jovanovich et al. 1983;Sheriff and Geldart 1995;Yilmaz 2001).
Common goals of reprocessing are to improve the data bandwidth, spatial resolution, signal-to-noise ratio, reflector continuity and, where relevant, seismic-borehole correlations (Sadhu et al. 2008;Lille et al. 2017). In this study we aim to better image the salt structures, particularly the base salt, previously interpreted from the legacy processing as a flat surface lying around 2.7 km below the seafloor, at a water depth of 3.2 km (Dal Cin et al. 2016). Due to the complex overburden geology and the short far-offset (3.1 km) with respect to the target depth, the main imaging challenges include accurately resolving velocity variations, eliminating multiples and improving the signal-to-noise ratio at depth. We confront these challenges by outlining three key stages that should be systematically included in modern processing flows for similar marine seismic datasets: i) 'broadband' processing, ii) multi-domain denoising and demultiple, and iii) geologically guided velocity model building using iterative pre-stack migration and travel-time tomography.
Compared to traditional 'narrowband' seismic processing, 'broadband' processing aims to improve the spectral accuracy by expanding the data bandwidth and restoring frequency content attenuated by the source-and receiverside 'ghost' effect and by seismic absorption Lille et al. 2017). High frequencies are most strongly affected by seismic absorption, which preferentially attenuates the highest frequency parts of the spectrum (Futterman 1962;Sams et al. 1997). Recovering this information improves the resolution and results in more accurate true amplitudes, improving the performance of other data processing steps such as velocity model estimation and migration, and quantitative interpretation such as amplitude variation with offset (AVO) analysis, impedance inversion, and attribute analysis (Chopra and Marfurt 2007;Mavko et al. 2009;Amundsen and Zhou 2013;Lille et al. 2017). Conversely, lack of low frequency signal commonly results in poor focusing in the deep part of the section, as low frequencies suffer less from scattering and absorption, so they penetrate deeper and display a better trace-to-trace moveout coherence, allowing us to build a more accurate velocity model (ten Kroode et al. 2013). In marine seismic data acquired with an airgun source, the source signature is a combination of a relatively broad impulsive signal (approximately a minimum-phase wavelet), periodic oscillations caused by the so-called 'bubble pulse' and an inverted polarity 'ghost' multiple, caused by the time-delayed reflection of the signal from the sea surface (Ziolkowski 1970;Sheriff and Geldart 1995;Hegna and Parkes 2011;Watson et al. 2019). GI-guns are designed to increase the primary-to-bubble ratio by injecting air inside the bubble generated by the primary pulse so that it collapses quickly (Pascouet 1989). Here, the GI-guns were operated in harmonic mode, where generator and injector volumes are equal, thereby reducing the bubble oscillation, without completely suppressing it (Landrø 1992;Dondurur 2018). The data bandwidth is principally widened by performing a source designature, whereby the primary airgun impulse and the bubble pulse are collapsed into a sharp, zero-phase wavelet (Sheriff and Geldart 1995;Amundsen and Zhou 2013;ten Kroode et al. 2013;Baldock et al. 2013;O'Driscoll et al. 2013). Deterministic deconvolution, where the operator is designed using an estimated source wavelet, can yield geological imaging superior to traditional statistical deconvolution methods, particularly for recovery of low frequencies and the preservation of amplitude information (Yilmaz 2001;Sargent et al. 2011;Scholtz et al. 2015;Davison and Poole 2015). Deghosting, instead, aims to deconvolve both the source-and receiver-side ghosts from the wavefield, further sharpening the wavelet and removing the frequency 'notches' associated with the ghost effect (e.g., Chuan et al. 2014;Davison and Poole 2015;Tyagi et al. 2016;Willis et al. 2018).
The quality of the bandwidth enhancement depends on the signal-to-noise ratio of the input data (Amundsen and Zhou 2013). It is therefore essential to attenuate as much as possible the low frequency noise (e.g., reverberation from the direct arrival, 'swell' noise caused by pressure fluctuation near the sea-surface) beforehand, because any remaining low frequency noise not correlated with the source pulse may be artificially boosted by the source deconvolution and deghosting filters (Yilmaz and Baysal 2015). Thanks to lower computational costs in recent decades, multi-channel filtering and analysis in transform domains has become routine for noise reduction (Schultz 1985). For example, 'swell' noise can often be better attenuated by predictive filters in the frequency-space (F-X) domain than by a simple time domain low-cut filter that results in loss of the low frequency signal along with the attenuated noise (Liu and Goulty 1999;Schonewille et al. 2008). Multiple attenuation is also better tackled by move-out discrimination techniques in the parabolic Radon domain rather than by traditional statistical deconvolution based methods, for example (Basak et al. 2012;Verschuur 2013). A 'broadband' processing flow combined with an efficient multi-domain noise separation can improve the signal-to-noise ratio in the deep part of the section (in our case below the MSG), allowing considerable improvements in velocity model building (Chuan et al. 2014).
Seismic imaging restores the correct geometry of seismic reflectors and requires an accurate velocity model of the subsurface (Jones and Davison 2014;Jones 2015). The legacy SALTFLU data were imaged using a Kirchhoff prestack time migration. Time domain migrations are relatively robust to errors in the velocity model but are only wellsuited to imaging geology containing weak lateral velocity variation (i.e., approximately 'layer cake' geology), as they do not properly account for ray path refraction. This can lead to degradation in image quality in time domain images of complex geology such as salt diapirs. Depth domain migrations, instead, can more accurately reproduce the ray paths of reflections in the subsurface, but the image quality is more sensitive to velocity errors (Sheriff and Geldart 1995;Yilmaz 2001;Jones and Davison 2014). Migration velocities are generally estimated based on the flatness of reflections on common midpoint gathers after migration (Tsvankin and Thomsen 1994;Jones 2015). In depth domain, the process of picking reflectors is often automated and used as input to travel-time tomography. Iterative rounds of analysis of the residual curvature of reflectors on depth migrated gathers followed by travel-time tomography to calculate model updates (Jones 2015).
In this study, we aim to showcase a 'broadband' reprocessing strategy designed to improve imaging of the MSG for the SALTFLU dataset. We demonstrate multi-domain de-noising, deghosting and a source designature using a seabed reflection derived operator. We then perform multiple attenuation and geologically driven iterative migration velocity analysis. Our results include pre-stack time and depth migrated images, which we compare to the legacy 'narrowband' processing. These reprocessed images highlight some new geological insights that these new results provide on the salt system of the Algerian basin, the seismic expression of the MSC and the basement structure of the study area.

Geological setting
The study area is located south of the Balearic Islands (Spain) with water depths between 1000 and 2800 m (Fig. 1a). The area covers the transition from the Formentera basin, on the Balearic Promontory, to the deep Algerian basin, which is marked by the steep NE-SW Emile-Baudot Escarpment, a NW-SE volcanic transform fault system (Acosta et al. 2001). The Algerian basin has previously been described as a Neogene back-arc oceanic basin that opened in response to the roll-back of the subduction of the Alpine-Maghrebian Tethys (Rehault et al. 1984;Verges and Sabat 1999;Mauffret et al. 2004).
There are no boreholes within the study area to tie to the seismic profiles or to constrain the velocity model. The nearest borehole, Alger-1, is located 60 km to the south-east (Fig. 1b), in a perched basin, at only 100 m water depth, in a different geological setting (Burollet et al. 1978). The closest borehole located in the deepwater Algerian basin is the ODP Site 975 (Fig. 1b), more than 200 km away from the line, that penetrated only the uppermost Messinian sediments (Comas et al. 1996).
Previous studies based their geological interpretation on the comparison with the analogous western Mediterranean margins, describing an Oligo-Miocene to Plio-Quaternary sediment cover, with the presence of thick (locally up to 2 km) Messinian evaporitic units both on the Balearic Promontory, and in the deep offshore Algero-Balearic basin (Fig. 1b;Wardell et al. 2014;Driussi et al. 2015;Dal Cin et al. 2016;Camerlenghi et al. 2018;Raad et al. 2020). These used the nomenclature from Lofi et al. 2011Lofi et al. , 2018, with a Mobile salt Unit (MU) identified in both the deep Algerian basin and the intermediate depth Formentera basin. In the Formentera basin, the MU separates the Bedded Unit 2 (BU2) from the Bedded Unit 3 (BU3; Raad et al. 2021). Laterally, when the salt pinches out, the BU2 turns into the Bedded Unit 1 (BU1) and/or the BU3. These bedded Units b Map of the same seismic lines acquired in the area superimposed to the present-day spatial extent of the MSC marker, with DSDP and ODP drillsites (in red the ones that sampled MSC evaporites) and industrial well Alger-1 (in black) (modified from Lofi 2018). Topo-graphic highs on the sea-bottom illustrate the presence of piercing diapiric structures. SF refers to SALTFLU lines presented in this paper, SB05 refers to survey SBAL-DEEP 2005, MS046 to profile 46 of the Mediterranean survey. c Schematic cross-sections across the Algerian basin showing the presence of complex salt structures below the surface (modified from  have been interpreted as Evaporites (primarily gypsum) mixed with clastic sediments (Ochoa et al. 2015;Raad et al. 2021). In the deep basin, the MU is highly deformed, with complex and steep salt structures locally deforming the seafloor, and separates the Messinian Upper Unit (UU) from the pre-salt unit (Camerlenghi et al. 2009;Dal Cin et al. 2016;Camerlenghi et al. 2018).

Methods
Reprocessing of the SALTFLU dataset is based on the following major processing stages ( Fig. 2): 1. Noise attenuation 2. Bandwidth enhancement, including source designature and source-and receiver-side deghosting. 3. Multiple attenuation, pre-and post-migration. 4. Iterative pre-stack migration and velocity model building in time and depth domains. 5. Post-migration processing, including seismic attenuation compensation.
Processing is performed using the REVEAL software from Shearwater Geoservices. After loading, raw data are checked for integrity using header values (Field File Identification Number (FFID), shot point number, channel number) and cross-checked with the observer logs to assess FFIDshot point integrity. Near-trace plots are visually analysed to identify any problems with the source array (e.g., misfires, missing shots) and receiver channels (e.g., noise bursts, electrical interference). The acquisition setup did not record information on the positioning of the streamer, preventing accurate location of the common mid-points (CMPs). Instead, we use the ship positioning GNSS data and project the streamer along the smoothed ship track, resulting in CMP locations that roughly honour the ship track and likely feathering of the streamer. The CMP binning is done using 6.25 m bins along the profile, giving a nominal fold of 60. After merging the navigation geometry into the trace headers, we subtract the recording time shift (40 ms), calculated as the average zero-offset time delay of the direct arrival (linearly extrapolated from the recorded data with a water velocity of 1511 m/s corresponding to the slope of the direct arrival in the data).

Noise attenuation
This stage aims to attenuate noise that could affect the deghosting and the estimation of the source signature. Firstly, data are resampled from 2 to 4 ms. A 3 Hz/30 dB to 110 Hz/96 dB Butterworth bandpass filter is applied in the frequency domain (Fig. 3b). The low-and high-cut values are chosen after visual analysis of the octave frequency panels to ensure no elimination of signal (Fig. 3a). The 3 Hz low-cut attenuates the very low frequency (1-3 Hz) part of the swell noise caused by pressure variations along the streamer (Bedenbender et al. 1970;Dondurur 2018). The 110 Hz high-cut filter is an anti-aliasing filter for later resampling to 4 ms (Fig. 3d, e).
The remaining low frequency part of the 'swell' noise is addressed by several iterations of frequency-space (F-X) domain prediction filtering (Hornbostel 1991;Schonewille et al. 2008). We perform up to 5 iterations on commonreceiver gathers, followed by 5 iterations on common-shot gathers. The first two iterations target broad noise bursts at all frequencies, while later iterations target local noise bursts with frequencies between 0-3 Hz, 0-5 Hz, and 6-12 Hz respectively. The total number of iterations depends on the amplitude of the 'swell' noise for each individual profile. The F-X prediction filtering has the secondary effect of attenuating amplitude spikes present in the dataset (e.g., from electrical interference).
After this initial filtering and resampling, missing shots and traces removed during shot and channel edits (Sect. 2) are interpolated on common-receiver gathers. This is performed by frequency-wavenumber (F-K) domain interpolation in two steps: 21 traces × 500 ms patches to target small gaps, then 125 traces × 300 ms patches to target larger gaps. Finally, forward and inverse linear Radon (τ-p) transforms are performed on shot gathers to further reduce the remaining 'swell' and random noise (Fig. 3d).

Deghosting and source designature
Deghosting is performed on shot gathers in the F-K domain to attenuate the source-and receiver-side 'ghost' effect resulting from sea surface reflections. Shot gathers are transformed into the F-K domain, where the notch location varies as a function of frequency and horizontal wavenumber, illustrating that the ghost time-delay is dependent on the emergent angle (Amundsen 1993;Day et al. 2013). The optimal deghosting operators for each frequency and wavenumber are estimated through a least-squares inversion scheme, minimizing the error between the model and the recorded data, parameterized based on the acquisition and medium parameters (i.e. receiver depth, water velocity, sea surface roughness). An inverse two-dimensional (2D) Fourier transform is then 13 Page 6 of 28 performed to return the deghosted data in the time-space domain. The deghosting is window-based to account for the non-stationary nature of the ghost in time and space. Prior to deghosting, shot gathers are padded by inserting 30 near-offset traces (362.5 m) and 10 far-offset traces (112.5 m) to prevent amplitude leakage from the forward and inverse transforms (Yilmaz 2001). The receiver-side ghost is estimated in local patches (500 ms × 11 channels), and the source-side ghost is estimated in broader time patches (1000 ms × 240 channels). After deghosting, a further F-K filter is applied to remove artefacts that may be introduced during the deghosting, particularly at the edge of the gathers. A regularization parameter is applied (0.01 for the receiver side, 0.001 for the shot side) to prevent the amplification of low frequency noise (Wang et al. 2017;Denisov et al. 2018).
After deghosting, we apply a deterministic source designature to collapse the bubble pulse, sharpen the wavelet and zero-phase the data (Fig. 4a, b). We follow the workflow of Sargent et al. (2011), where the far-field source signature is derived from the flattened and stacked waterbottom reflection (Fig. 4c, d). A matching filter is then designed to convert the estimated far-field source signature to a zero-phase band-limited Ormsby wavelet (4-8-70-120 Hz, Fig. 4e). Using a single deterministic operator yields better sub-seabed imaging and low frequency recovery superior to probabilistic methods (e.g., surface consistent, spiking and predictive deconvolution) and modelled source methods (Sargent et al. 2011;Davison and Poole 2015;Maunde et al. 2017). Converting the seismic data to zero-phase tends to improve temporal resolution and lowers interference between closely spaced reflectors, allowing for better delineation of the stratigraphy, and is a necessary pre-condition for Kirchhoff-type migrations (Sheriff and Geldart 1995).
Finally, an inverse-Q filter (phase only) is applied immediately before migration to compensate for dispersion effects caused by seismic attenuation (Sams et al. 1997). Dispersion introduces a phase-shift with increasing propagation depth. Compensating for this is necessary to ensure that the wavelet remains zero-phase with increasing propagation depth (Wang 2006). A two-layer Q model is used, with Q = 1000 in the water layer (essentially non-attenuating). We choose a value of Q = 70 in the subsurface, based on a constant single value that best flattens the frequency spectrum and balances amplitudes with depth across all the profiles in the survey.

Multiple attenuation
The water depth exceeds 1500 m in the survey area and the free surface 'water bottom' multiples generally arrive later than most of the reflected primary signal at the target depths. High amplitude short period internal multiples are, however, showing bandpass filtered frequency intervals (each one octave) of the raw field data; b Input near-trace-plot (left), outputs after low-cut filter (middle) and high-cut filter (right); c Frequency spectra cor-responding to b; d input shot gather after resampling to 4 ms (left) versus after linear low-frequency noise removal steps (middle), and difference panel between the two (right); e frequency spectra corresponding to d observed throughout the survey, generally linked to the strong impedance contrasts associated with the salt (Fig. 5). We therefore attempt to attenuate both long and short period multiples using moveout discrimination techniques.
Before migration, we perform multiple attenuation in Radon domain (Fig. 5). Normal moveout (NMO) correction is performed using the RMS velocity model (Sect. 3.4.1) to flatten primary energy, and the CMP gathers are transformed into the Radon domain using a least squares parabolic Radon transform (Sacchi and Porsani 1999). The Radon domain multiple model is defined by a time-varying inner moveout mute (Fig. 2). A top mute is applied at 200 ms below the water bottom to avoid attenuating primary energy in the shallow Plio-Quaternary (largely free of internal multiples). We then perform the inverse Radon transform on the multiple model and subtract from the original data.
We perform a second pass of Radon multiple attenuation post-migration, with a similar workflow to the pre-migration demultiple. For the depth domain processing (Sect. 3.4.2), the depth migrated Common Depth Point (CDP) gathers are first converted to time domain using the time migration velocities, the Radon demultiple is applied, and the CDP gathers are converted back to depth domain using the same velocity field. The multiple model contains a weak component of the primary signal (Fig. 5c), but this 'primary leakage' is considered acceptable as it does not significantly alter the amplitude of the reflectors after stacking.

Time domain
Our approach to velocity model building in the time domain involves picking RMS velocities on semblance panels on migrated CMP gathers. We perform several iterations of migration velocity analysis, gradually refining the velocity field.
A sparse initial velocity model is picked from semblance analysis of unmigrated CMP gathers every 500 CDP (~ 3 km). This initial velocity model is primarily used as a Fig. 4 Results from the bandwidth enhancement stage. a Near-trace gather before deghosting (left), after deghosting (middle) and after source designature (right); b frequency spectra corresponding to a; c frequency spectra of the far-field source signature estimated from the stacked seabed and the target Ormsby wavelet displayed on d and used for designing the source designature (debubble and zero-phasing) operator; e the resulting source designature operator stacking velocity to quality control intermediate processing stages, and for the initial pre-stack time migration.
We use a Kirchhoff pre-stack time migration to image the data. This choice is dictated by a limited computing resources and the flexibility of this algorithm to velocity errors. Before migration, the data are regularised into evenly spaced offset bins of 50 m using a partial NMO correction. The migration aperture is parameterised by analysing the migration impulse response (Appendix A). Aperture is broadened from 1500 m at 1500 ms to 3500 m at 5000 ms. After the migration, we perform Radon domain multiple attenuation (Sect. 3.3) with broad parameters (minimum moveout − 100 ms, maximum moveout + 300 ms) to avoid attenuating primary reflections for the following velocity analysis.
We perform migration velocity analysis on the migrated CMP gathers every 250 CMP (~ 1.5 km). The migration followed by velocity analysis is repeated three times until the primary reflections are flattened on the migrated gathers.

Depth domain
Depth domain processing aims to estimate an interval velocity model and produce depth migrated seismic images. Similar to the time domain model building, we use an iterative approach by repeatedly performing velocity analysis on migrated gathers, gradually refining the velocity field. Our criteria for a successful velocity model are that the primary reflections are flat after depth migration.
We use a Kirchhoff pre-stack depth migration to migrate the data in depth domain. The input data to the migration are the regularised CMP gathers (after Radon demultiple and phase-only inverse-Q filtering) as were used in the pre-stack time migration (Sect. 3.4.1). The parameters used for traveltime computation are detailed in Appendix B.
The initial velocity model is built using a water velocity (which we assume to be constant) and sediment velocities derived from the time domain velocity analysis (Sect. 3.4.1). The water velocity is estimated by performing a range of constant-velocity depth migrations (1450-1550 m/s). We choose a constant water velocity of 1525 m/s for the whole SALTFLU survey, which best flattens the water bottom reflection across the different profiles and at different water depths. To derive the interval sediment velocities, we first flatten the time domain velocities along the water bottom. The RMS velocity models in time domain are then converted to interval velocities in depth domain using a Dix conversion (Dix 1955). The resulting interval velocities are Gaussian smoothed (100 m × 500 CDP) before shifting back to the original water bottom depth. An initial depth migration is performed using this initial velocity model and the input gathers to the final pre-stack time migration.
Following this initial pre-stack depth migration, we obtain velocity updates using a global tomography procedure comprising automated offset-dependent residual moveout picking, ray tracing and travel-time tomography (proprietary to REVEAL software). The velocity update in each iteration is limited to -10% and + 20% of the input velocity, respectively. The automated picks are quality controlled by visually checking that the picks correlate with major reflectors, eliminating picks with low trace-to-trace correlation and eliminating picks with anomalously high moveout correction. Grid-based travel-time tomography is then used to compute model updates from these edited picks. The parameters used for tomographic updates for each velocity model building iteration are detailed in Appendix C. We consider that the velocity model is good enough when most reflectors appear flat on migrated gathers. The first two iterations are performed only for the post-salt sediments, updating only between the interpreted water bottom and the top salt horizons. Initially, the top salt horizon is interpreted from the time domain data and converted to depth using the initial velocity model. The picking of the top salt horizon is refined for each velocity updates to fit the new depth image and limit the tomographic iteration to the post-salt.
The sediment model (Fig. 6a) is then flooded with salt velocity from the top salt horizon to the base of the profiles (Fig. 6b) The transition from the sediment velocity to the salt velocity is Gaussian smoothed (50 × 50 m) to enable projected rays to pass through the sharp velocity contrast during travel-time computation for the migration and tomography. The base salt horizon is interpreted on the depth migrated stacks after salt flooding, and a new velocity model is built with a salt velocity flooding that stops at this newly picked base salt. The pre-salt is also flooded with a velocity gradient starting from the base salt horizon. We test the velocity with a range of starting velocities and velocity gradients and assess the overall flatness of reflectors in the pre-salt. After testing, the best results were obtained with a gradient of 2700 m/s + 1.2 (m/s)/m. Finally, the final interval velocity model is obtained after two more iterations of tomographic velocity updates following the pre-salt flooding (Fig. 6c). After the final depth migration, Radon domain demultiple is performed as previously described for the time domain processing (Sect. 3.3).

Post-migration processing
After the final migration, residual moveouts are automatically picked every 4 CDP (~ 20 m) using semblance velocity spectra. The residual moveout analysis is limited to 5% change, and the resulting picks are smoothed spatially (25 ms × 200 CMP). The residual moveout field is applied to the migrated CMP gathers before stacking. Additionally, we apply an inverse-Q filter (amplitude only) to compensate for the frequency-dependent amplitude attenuation caused by seismic absorption. The Q model is the same as was Page 11 of 28 13 previously for the phase-only inverse-Q filter (Sect. 3.2). We then apply an inner angle mute (at 1.2°) to target nearoffset residual multiples and an outer angle mute (at 48°) to remove far-offset refractions and data affected by wavelet stretching. The CMP gathers are stacked, and a low-cut frequency filter is applied (Ormsby, 500 ms, 4-12 Hz) in a short 400 ms window below the seabed to target any residual low-frequency noise post-stack. Resulting outputs are displayed on Figs. 7 and 8.

Results
The resulting seismic sections can be plotted using free software packages such as Seismic Unix, 1 SeiSee 2 or Opend-Tect 3 etc. We assess the quality of the newly reprocessed images by comparing them with the pre-existing vintage images (Figs. 7 and 8). Figure 7a shows a section of SF06 on the Balearic promontory where the bubble noise and the multiples have been well attenuated compared to the vintage data (Fig. 7b). This allowed for better estimation of the velocity for migration by revealing new geological reflectors that were not well imaged before (between CMPs 20 000 and 21 000, around 2.3 s and around 2.6 s). Figure 7c shows a section in the deep Algerian basin where the pre-salt reflectors and lateral terminations of reflectors against the fault are better imaged than in the vintage data (Fig. 7d). The bandwidth enhancement appears to improve the low-frequency content of the data, boosting the amplitude of the base salt and the pre-salt reflectors (Fig. 8a, c). In the vintage processing images the base salt was very poorly imaged, and the pre-salt reflectors were nearly invisible (Fig. 8b, d). These improvements allow for better distinction between multiples and primary signal during velocity picking, improving the move-out based multiple elimination. Despite this, some residual multiple energy is still visible below the base salt (e.g., multiples overprinting the base of the salt and the pre-salt reflectors on Fig. 7c). In vintage processing sections, the steep faults are associated with high amplitude migration noise (migration "smiles" artefacts, Figs. 7f, 8c, d) that is significantly lower amplitude on the reprocessed data (Fig. 7c, 8a, b). This is achieved thanks to the use of better velocity model and a better migration algorithm compared to the vintage data. Improved imaging of steep reflectors is also demonstrated in Fig. 7e, where faults are better imaged than in Fig. 7f. The bandwidth enhancement also better balances the amplitude variations, emphasizing better the contrasts in acoustic impedance, particularly within the evaporites, and emphasizing geometries that were not well imaged before (e.g., terminations against the unconformity in the Plio-Quaternary on the left side of Fig. 7e). This is also well illustrated on Fig. 8a, b, where steep reflectors connected to complex salt geometries (e.g., salt-rooted faults or diapirs crests) are much better imaged than on vintage images (Fig. 8c, d), where migration artefacts mask the steeply dipping reflectors. The vertical resolution of the reprocessed data can be estimated using the Rayleigh criterion of λ/4 (where λ is the seismic wavelength) (Kallweit and Wood 1982). RMS velocities and dominant frequencies range from 1511 m/s and 60 Hz at the seabed, to 5500 m/s and 10 Hz at the presalt. The resulting nominal vertical resolution varies from approximately 6.3 m at the seabed to 137 m below the salt. The legacy time migrated stacks display a slightly higher dominant frequency of 65 Hz at the water bottom, hence a very slightly higher vertical resolution, notably in the first 500 ms below the water bottom (5.8 m using the aforementioned velocity at the seabed). However, the reprocessed data shows a much better imaging of the deeper MSG thanks to a broader bandwidth, a better signal-to-noise ratio, and an improved migration compared to the vintage data.
The pre-stack depth migration corrects for geometrical distortions linked to the strong lateral velocity variation between the evaporites and the surrounding sediments (Fig. 6). The pre-salt reflectors are not always clearly visible on the final sections, thereby limiting the reliability of the velocity model beneath the salt (see Sect. 5.2). The 'pullups' of the pre-salt beneath the salt structures seen in the time migrated section are flattened after depth migration with an average velocity of 4300 m/s.

Limits of the pre-processing strategy
The chosen pre-processing strategy is determined by the computing costs and our ability to separate noise from primary signal. The time resampling from 2 to 4 ms eliminates the frequencies above 125 Hz and attenuates the frequencies from 90 Hz. This is due to the smooth roll-off of the applied filter that is necessary to avoid ringing. By doing so, the highest resolution part of the data is lost but this allows a significant reduction of computing times. Resampling divides by two the size of the dataset and lowering the Nyquist frequency to 125 Hz filters out frequencies that would otherwise become aliased after resampling. That way, spatial transforms (e.g. F-K, Tau-P and Radon filters) can be performed without further trace interpolation to unwrap the aliased energy. Resampling allowed us to divide the number of samples (hence the data size) by 2. It also removes the need to interpolate between receiver channels done to avoid spatial aliasing of high frequencies (due to the acquisition receiver under-sampling in space). Therefore, resampling from 2 to 4 ms reduces the overall data size by 4, considerably improving computing times. The aim of this reprocessing is to improve the image of the evaporites and the salt generally lying below a thick Plio-Quaternary cover, therefore we consider this resolution loss an acceptable trade-off.
A future study focusing on the shallower sub-surface could benefit from the preservation of these higher frequencies by processing the dataset at 2 ms sample rate.
During the acquisition of SALTFLU, the source array and streamer were towed at approximately 3 and 4 m below the sea surface, respectively. Deep streamer towing increases the operational weather window, reduces noise, and increases signal penetration (Carlson et al. 2007). However, a timedelayed reflection from the sea surface (a so-called 'ghost') can interfere destructively with the recorded signal and limit the resolution of the data (Schneider et al. 1964). These ghost reflections introduce notches in the spectrum of the recorded wavefield which particularly attenuate the low frequency component (Yilmaz and Baysal 2015). Assuming a vertical ray path, the first frequency notch should appear at f = v/(2z), where z is the source or receiver depth and v is the speed of sound in water (Amundsen and Zhou 2013). A water velocity of 1511 m/s at the surface (estimated from the slope of the direct arrival), and source and receiver depths of 3 m and 4 m, respectively, result in a notch around ~ 252 Hz for the source-side 'ghost', ~ 189 Hz for the receiver-side 'ghost', and ~ 108 Hz for the combined source and receiver 'ghost'.
In the data, the 'ghost notch' varies between approximately 118-150 Hz on near-trace gathers (offset 100 m) depending on the lines, and it gradually shifts towards higher frequencies along channels until it disappears above filtered frequencies after the 10th receiver (offset 325 m). Such values broadly correspond to the combined source and receiver 'ghost'. The systematic variation in the dominant frequency of the ghost notches likely reflects the varying depth of the streamer with offset, despite the use of birds to attempt to maintain the whole streamer at the target depth of 4 m (Table 1). Poor streamer depth control can result from a slow vessel speed, rough weather, or strong currents perpendicular to the streamer. With a Nyquist frequency of 125 Hz after resampling, the source-and receiver-side 'ghost' notches are mostly outside the data bandwidth (Fig. 3c), but deghosting still helps to extend the lower part of the bandwidth and removes the interference of the 'ghosts' with the primary wavelet in time domain (Raj et al. 2016). A deeper towing depth would have been desirable to improve the recording of the lower frequencies and to attenuate the 'swell' noise, which is particularly strong and dominates the low frequency part of the data. It could not be fully eliminated pre-stack without attenuating some primary signal. The best results were obtained with iterative F-X prediction filtering, which better preserved the low frequency primary signal compared to the low-cut filter applied to the vintage dataset, but it required many iterations with different sliding windows on both common-shot gathers and commonreceiver gathers. Other approaches tested include F-K filtering and time-frequency filtering, but these were either unsuccessful at separating the noise from the primary signal, or too computationally expensive. Tau-P domain filtering is also an efficient way to separate noise from primary signal (Basak et al. 2012). Here, we apply a shot domain forward and inverse Tau-P transform, with a limited slowness range, without further muting in the Tau-P domain. The 'swell' noise is incoherent enough to not sum constructively during the transform and is therefore well attenuated. The remaining noise embedded in the data is low amplitude enough to ensure a stable deghosting operator estimation, but it is boosted after the deghosting stage and an additional denoising step is required after bandwidth enhancement (Fig. 2).

Limits of the multiple elimination and accuracy of the velocity models
As the water depth exceeds 1500 m in the survey area, long period surface-related multiples are not superimposed on the target primary signal. Therefore, surface-related multiple elimination (SRME; Verschuur 2013) is not applied, and all multiples (long and short period) are attenuated solely by move-out based methods. However, short period internal multiples are present and not fully eliminated (Fig. 5a-c). They are generated by strong acoustic impedance contrasts either within the Upper Evaporites, or between the water bottom and the Top Evaporites (Fig. 5d). They prevent the resolution of velocity variations below and within the evaporites and can be misinterpreted as primary signal. During the processing they are only partially attenuated by moveout discrimination on CMP gathers, likely because they are generated by a high velocity layer (the upper Messinian unit above the salt). With a maximum offset of 3100 m, they only yield a small move-out difference compared to the primary reflections. Other approaches for attenuating multiples include boundary-related internal multiple method (Verschuur and Berkhout 1996), data-driven internal multiple removal method (Jakubowicz 1998) and inverse scattering internal multiple removal (e.g. Araújo et al. 1994). The second method has been tested but it was not applied in the final flow due to significantly increased computational times (16 h for 5000 shots, while the survey consists of more than 170 000 shots).
The presence of remnant multiple can lower the accuracy of the time and depth domain velocity analyses. The time domain velocity model is built through several iterations of semblance-based migration velocity analysis (with a pre-stack Kirchhoff time migration, and Radon multiple attenuation). The remaining multiples and the out-of-plane reflections can generate high semblance values that can be wrongly picked. Consequently, automatic picks had to be carefully quality controlled during the residual move-out analysis stage after both the time and depth migrations. Internal multiples generated within the Upper Evaporites can also be wrongly picked as the top salt horizon used for the velocity flooding. Two reflectors R1 and R2 could be interpreted as the top of the salt (Figs. 9b and 10b).
The top of the salt is difficult to follow laterally because of the complex salt structures and the presence of several reflectors within the salt that could represent residual internal multiples, out-of-plane reflections or the presence of internal clastic beds within the salt. The transition from the Upper Evaporites seismic facies to the 'transparent' salt seismic facies is gradual. The impedance contrast between the base of the Upper Evaporites and the salt is not clearly marked. Based on the strength of the amplitude and the lateral continuity of the reflectors, R2 seems to be the most likely horizon marking the top of the ductile salt. However, despite the likely presence of internal multiples below R2, some reflectors are believed to represent geology because they display different geometries and seismic facies compared with overlying evaporites. Furthermore, intra-salt reflectors have already been observed in the deep Mediterranean salt (e.g. Bertoni and Cartwright 2015;Netzeband et al. 2006;Hubscher et al. 2007;Feng et al. 2016). We interpret these events as intra-salt reflectors generated by the presence of clastic beds with the mobile salt. Intra-salt reflectors have already been observed in the eastern Mediterranean, where they lower considerably the average salt velocity to ~ 4260 m/s (Feng and Reshef 2016). On common reflection-point gathers, these reflectors align at a velocity of ~ 3250 m/s (Fig. 10) and would be over migrated if salt flooding was performed from R2. To preserve them from the post-migration Radon demultiple, the salt flooding is performed starting from reflector R1. This choice has a strong control on the velocity used for the flooding to flatten the base salt (a lower top salt implies a higher flooding velocity to flatten the base salt), the resulting geometry of the salt structures on the stacked sections, and the depth of the base salt and pre-salt units. A better identification of the top salt could be guided by drilling through the UU until the top of the salt in the Algerian basin.
The final interval velocity models used for Kirchhoff depth migration are sampled every 12.5 m laterally and 5 m vertically. Tomographic updates were considerably smoothed and regularized to avoid velocity 'bulleyes'. In the absence of nearby wells, the only quality check on the interval velocity model is the flatness of reflectors in common reflection point gathers after migration. However, considering the limited offset of 3.1 km compared with the depth of the target (the base salt is > 4 km below sea level for most of the survey area) and the strong velocity variations, the reflections are flattened under a substantial range of velocities. This implies that there are large uncertainties in the velocity model, particularly for the evaporites and the pre-salt, where the residual internal multiples and the lack of strong pre-salt reflections make  it difficult to constrain the velocity model. The ray-based Kirchhoff migration is not the best suited strategy for salt imaging due to (Leveille et al. 2011): (i) a limited ability to handle steep and large velocity variations which requires a smoothing velocity model and limits the quality of the output images and (ii), the apparition of migration "swing" artifacts in poorly illuminated zones under or near salt bodies due to limitations of the assumption of the ray-based scheme. More advanced migration techniques (e.g., anisotropic migration or reverse time migration) have not been tested. They are more costly and are accurate with sufficiently detailed velocity models, thus they require corresponding higher resolution velocity estimation techniques (e.g., waveform inversion or anisotropic model building; Jones 2015). These techniques require dense and long-offset recording of the reflected wave field, a good knowledge of the geology and petrophysics distribution within the medium, with ideally nearby well control, and good computational power (Virieux and Operto 2009). SALTFLU data is sparse (~ 6 km spacing between the lines), with no nearby wells, and the maximum offset recorded is ~ 3.1 km, while the local subsurface is known to include a highly deformed and deep (> 4 km) salt . For these reasons, the Kirchhoff migration is preferred (Yilmaz 2001;Vardy and Henstock 2010;Dondurur 2018). In practice, it remains a flexible and robust approach where the velocity model is not well constrained and the objective is primarily the imaging of salt structures, not the subsalt imaging, but reverse time migration could potentially improve the imaging of the salt structures.

Area of interest within the central Algerian basin observed on the reprocessed data
The newly reprocessed data better image the pre-salt reflectors, allowing an improved understanding of the salt tectonic system from late Miocene to today in the central Algerian basin , with a peak of contractional salt deformation at early to mid Plio-Quaternary. This peak could be linked to an important episode of shortening of the Tell-Atlas fold-and-thrust belt, along the Algerian Margin. A contemporaneous tectonic event has been also identified along the Balearic margin, with a regional unconformity, and the reactivation of pre-existing faults leading to the uplifting and the tilting of the Balearic slope .
SALTFLU is one of the highest seismic resolution datasets available in the Algerian basin. It allows detailed study of the late Miocene seismic facies, the acoustic basement in the Formentera basin, and faulting within the Plio-Quaternary that could provide new insights on the geological evolution of the region. Combined with other neighbouring datasets, regional sections extending from the central Mallorca depression to the Algerian margin through the Formentera basin and the deep Algerian basin could be drawn.

Fluid circulation and mud volcanoes
We observe several amplitude anomalies that we interpret as fluid indicators in the south-central Algerian foredeep (Fig. 11a, b) and along the Balearic margin (Fig. 12a,  b). In the Algerian foredeep we interpret pushdowns and amplitude dimming below high amplitude anomalies as gas chimneys (along black arrows on Fig. 11). They are systematically located above piercing salt diapirs. This suggest that they need a migration path through the evaporites to escape vertically into the Plio-Quaternary. Somehow, they do not migrate further up dip along the Plio-Quaternary strata (Fig. 11b), suggesting the potential presence of stratigraphic traps. It is not known if this gas is thermogenic or biogenic in origin. Previous studies in the Mediterranean sea have shown that cross-evaporite fluid flow through an hydro-fractured MU is likely, implying that the gas observed in the SALTFLU profiles could potentially be sourced from a pre-salt source rocks (Dale et al. 2021;Oppo et al. 2021).
A thermal modelling study in the eastern Algerian basin from Arab et al. (2016) favours the pre-salt Oligo-Miocene sourcing, with the Messinian shale as a possible biogenic gas source. In the Algerian foredeep, seismic fluid indicators are only observed in the western part of the study area, which is characterized by relatively undeformed Plio-Quaternary deposits, and less deformed salt compared to the eastern part . The absence of gas chimneys in the eastern part could be due to the absence of traps, and/or migration pathways, and/or the immaturity or absence of a source rock. If there were fluids migration but no traps, we would expect to see evidence of fluid escape features, such as pockmarks, at the seabed. Arab et al. (2016) conclude that discharge and accumulation of the pre-salt thermogenic fluids is limited to the Algerian margin or beyond the slope toe (60 km from the coastline), with some structural traps associated with the Quaternary north-verging thrust ramps. Further eastward, the basin widens and the thrust ramps are located further away from the seismic data .We speculate that the fluids may not be able to migrate far enough within the basin to be observed locally of the eastern lines. Since the western part of the Algerian basin also displays a higher heat flow than the eastern part (Poort et al. 2020), this spatial variation could have had an impact on the production of fluids, whether biogenic or thermogenic.
Along the Balearic margin, shallow high amplitude variations could suggest the presence of gas pockets (Fig. 12a). Locally, disrupted and dimmed reflectors below a locally arched overburden could record the presence of fluid migration paths (Fig. 12b). The nature and the source of these fluids are also unknown. Previous studies link these seismic amplitude anomalies and the normal faulting on the Balearic margin to hydrofracturing induced by the fluid circulation (Urgeles et al. 2013;Wardell et al. 2014;Del Ben et al. 2018;Dale et al. 2021).
The presence of mud volcanoes in the Algerian basin has previously been suggested by Camerlenghi et al. (2009). In salt basins, shale diapirs are common in areas of shortening, which can pressurize and mobilize shale, as can vertical loading (Jackson and Hudec 2017). Mud volcano systems have been recognized in the neighbouring Alboran sea, where they are sourced by pre-Messinian sediments (Aquitanian-Burdigalian in age; Sautkin et al. 2003;Blinova et al. 2011;Medialdea et al. 2012;Somoza et al. 2012). We observe several weakly reflective diapir contacts, with no underlying velocity pull-up but with a narrow dimming zone instead (Fig. 11b). Migrating this structure with a high salt velocity result in over-migration artefacts, suggesting that either (i) they are not made of high velocity material, but could be made of low velocity mud; (ii) the structures are out-of-plane with respect to the 2-D seismic profiles; or (iii) the salt is allochthonous and too thin to generate a pull-up. They could testify to an active mud volcano system, but the evidence is too poorly constrained to verify this hypothesis.

Seismic facies variations within the Messinian seismic units
Along the Balearic margin, the Upper Unit is interpreted as an evaporite-rich layer, with two subunits UU1 and UU2 separated by an intra-UU unconformity, and possibly, an intra-UU salt layer (Dal Cin et al. 2016;Camerlenghi Fig. 12 Sections in depth of lines SF03 (a) and SF08 (b) showing amplitude anomalies and signal disruption that may indicate the presence of fluid migration along the Balearic margin. Arrows indicate areas of blanking and/or disturbed bedding, representing possible fluid migration pathways. Vertical axes indicate depth below the sea level. Vertical exaggeration × 3. Positioning map from Fig. 1 et al. 2018). On the newly processed data, UU1 in bounded by R4 and R3, and UU2 is bounded by R3 and R2. The intra-UU unconformity would correspond to R3.
We question the presence of an intra-UU salt layer and of an intra-UU erosion surface along R3, as reflector truncations along R3 are not observed regionally in the dataset. The Reflector R3 is the strongest amplitude reflector observed within the whole survey (Figs. 9b and 11). It is the reflector that would separate the UU1 from the UU2 described by Camerlenghi et al. (2018). This regional strong impedance contrast could mark the part of the UU containing the largest proportion of evaporites and could record a climax for the UU stage. Velocity profiles do not show evidence of a high velocity layer in the UU (Fig. 10). On the contrary, the average velocity of the UU is low: it increases downward, from ~ 2500 m/s at R4 to ~ 3250 m/s at R1. These values cannot be associated to high velocity evaporites such as gypsum and/or anhydrite (~ 5000 to 6000 m/s; Schreiber et al. 1973;Mavko et al. 2009) or halite. This could indicate a clastic-rich layer, with an increase in the proportion of low-velocity material, where the last unit of UU is predominantly made of clay, such as in the Unit 7 of the Levant Basin (Gvirtzman et al. 2017). The scientific wells that penetrate the topmost UU in the Algerian basin confirm the presence of a sequence of mud and marls, interbedded with evaporites (dolomites, gypsum and anhydrite) that are possibly rich in organic matter (ODP-975, DSDP-124 and DSDP-371;Ryan et al. 1973;Hsü et al. 1978;Comas et al. 1996).
The shallowest Upper Unit (UU1), between R4 and R3, exhibits lateral amplitude variation, both toward the Algerian and the Balearic margin (LAV on Figs. 9 and 11). The reflectors of this unit lose their high amplitude expression toward the nearest Balearic margin and toward the deepest Algerian basin. This loss of amplitude could result from a variation in the acoustic impedance contrast, that could be related to a change in lithology or pore fluid content. This uppermost unit of the UU also thins down toward the margin. It nearly falls below the seismic resolution (Sect. 3) and is then expressed as a medium to strong amplitude reflection with negative polarity compared to the water bottom. It is not clear where it pinches out, but it appears to do so before the underlying unit (UU2), between R3 and R2. This last underlying unit also displays a sudden lateral variation (Figs. 9b,13), pinching out towards the margin. The point where most reflectors onlap the R3 reflector could coincide with the point where the overlying unit falls below the seismic resolution. These sudden changes could indicate a different depositional environment, either due to local variation in accommodation space, sediment supply and/or salinity.
The lowermost seismic unit of the UU, interpreted in this study between reflectors R1 and R2, is characterised by relatively discontinuous medium amplitude reflectors embedded in a low amplitude seismic background 'matrix', with a velocity comparable with the overlying units (3000 ± 250 m/s; Fig. 9b). This unit could represent a transition facies from the ductile salt (dominantly halite with poor clastic content) to more brittle Upper Evaporites, richer in clastic and gypsum evaporites. Many faults can be tracked down to the reflector R1 which attests to brittle behaviour (Fig. 9b). Whether it is part of the UU or the MU, this unit pinches out before the overlying units.

Seismic expression of the Messinian Lower Unit in the central Algerian basin
The absence of a high amplitude Lower Unit (LU) below the salt, as it is expressed in the Provencal basin (Lofi et al. 2011), could indicate the absence of the LU in the Algerian basin; that the LU is much thinner than the vertical resolution of the data; or that there is an absence of strong impedance contrasts within this unit and within the pre-Messinian units. In the Levant basin, in the eastern Mediterranean Sea, the LU is also absent and the onset of the Messinian Salinity Crisis is marked by the Foraminifers Barren Interval, a 10 s-of-meters thick, evaporite-free, shale unit that records the entire duration of the first stage of the crisis (Manzi et al. 2021).
It could be that the record of the onset of the MSC along the Balearic margin of the deep-water Algerian basin is similar to the record of the MSC in the Levant basin. Toward the southern Algerian margin, however, the expression of the pre-salt units changes laterally and displays a set of continuous low amplitude and low frequency reflectors that could record the LU (Fig. 13). Such seismic facies have been associated with the LU in previous studies, where it would be made of 100-200 m of plastic grey marls and gypsum based on nearby industry wells (Burollet et al. 1978;Medaouri et al. 2014). However, they seem difficult to differentiate from the Tortonian underlying units (Leprêtre 2012). This could confirm a lack of impedance contrast between the two (nearby wells showed that the underlying Tortonian and Serravalian units was also made of gray marls, with intervals of pyroclastics sands or limestones; Medaouri et al. 2014) suggesting that in the Algerian basin, the LU, when present, is not as rich in evaporites that in the Provencal basin. In that case, the LU could be a lateral equivalent of the Foraminifers Barren Interval, deposited preferentially in the deepest part of the Algerian basin. Alternatively, contemporaneously with the deposition of the LU in the Provencal basin, MU is already being deposited in the Algerian basin, and LU is a lateral equivalent of the lower MU (Roveri et al. 2019).

Insights on the tectonic setting of the Balearic promontory
The new SALTFLU sections also provide improved imaging of the Formentera basin, which could allow better stratigraphic correlation with the Messinian units from the Mallorca depression and the deep-water Algerian basin and a better understanding of the tectonic history of the Balearic promontory (Figs. 13, 14a).
When aligning the seismic data to the top surface of the evaporites (Fig. 13), it seems straightforward to correlate the BU3 with the UU above the salt. Following the interpretation of Raad et al, (2021), if the BU1/BU2 (contemporaneous in age) are equivalent to the Stage 1 primary Lower Gypsum, according to the correlation suggested by Roveri et al, (2019), the equivalent of the BU1/BU2 in the deep offshore is the thin Foraminifers Barren Interval (likely below the seismic resolution). This could explain why a pre-salt evaporitic unit is not observed in the deep basin, while it is recorded by the BU1/BU2 in the intermediate depth basin.
In the reprocessed SALTFLU dataset, the pre-salt basement appears highly variable, containing several sharp and highrelief structures inferred to be of volcanic origin (Fig. 14b). These structures match the geometry and the seismic expression of igneous bodies in volcanic provinces, such as the Taranaki basin (Infante-Paez and Marfurt 2017). We interpret several magma complex and intrusive sills beneath the Miocene to the current sedimentary cover of the Formentera basin (Fig. 14b). They are overlain by a highly deformed unit, that is onlapped by pre-Messinian reflectors wedging towards the top of the magma conduits. These wedges indicate that the volcanoes were already present during the deposition of these pre-Messinian units. Extensive faults within the Plio-Quaternary cover suggest that the Formentera basin is still actively deformed (Fig. 14b).
Collapse structures are observed at the top of the BUs, in the deepest part of the Formentera basin (Fig. 14c). The top surface of the BU3 unit is generally conformable with the Plio-Quaternary in the Balearic promontory (Raad et al. 2021) and is erosive only on topographic highs, not in the depressions where salt is preferentially accumulated. This suggests that the depressions were never aerially exposed to erosion or dissolution. Instead, the collapse features could be due to the circulation of undersaturated fluids or diagenesis, similar to features observed in the Messinian evaporites in the eastern Mediterranean (Bertoni and Cartwright 2015). Heat flow measurements on the Balearic promontory record anomalously low values that support the presence of groundwater fluid circulation (Poort et al. 2020).

Conclusions
We present the reprocessing of 2-D airgun seismic reflection data acquired offshore the Balearic Islands in the Algerian Basin in 2012. The goal is to improve the imaging of the Messinian evaporites and the overall basin architecture. The reprocessing strategy is designed to better image the pre-salt reflectors in an amplitude-preserving manner, attempting to overcome the challenges of imaging complex and deep subsurface geology using offset-limited seismic data. This has been done using a 'broadband' processing strategy, multiple attenuation and an imaging approach that integrates geophysics and geological interpretation to iteratively build the velocity model. The resulting pre-stack migrated images, in depth and time, display improved reflector continuity and amplitude preservation, particularly for the pre-salt. The processing flow presented here could be applied to other offset-limited legacy airgun reflection seismic datasets in the deep Mediterranean Sea. The efficiency and efficacy of the workflow strongly depends on the original acquisition parameters: the limited offset of most vintage data (often < 3 km) compared with the depth of the salt structures (the base salt lies in between 3.5 and 5.5 km) inhibits our ability to separate signal from noise and to accurately resolve the subsurface velocity distribution using moveout-based processing techniques. The reprocessed data image several fluid indicators, lateral amplitude variations, salt structures and volcanic structures in greater detail and more precision. These new results provide insights into the evolution of the under-explored Algerian basin and the Messinian Salinity Crisis. These outcomes highlight the value of reprocessing legacy academic seismic data, particularly when considering how challenging acquiring new seismic and borehole data has become in the western Mediterranean Sea.
Appendix A See Fig. 15. c Zoom 2 along line SF08 shows the Messinian evaporites in the Formentera sub-basin, with several depressions at its top that suggest the presence of collapse structures, chaotic seismic facies in the Plio-Quaternary that suggests the presence of mass transport deposits (within dotted black circles), and disturbed signal below the Evaporites that could indicate fluid circulation. Vertical exaggeration × 5

Appendix B
See Table 2.

Competing interests
The authors have no relevant financial or nonfinancial interests to disclose.
Ethical approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit https:// creat iveco mmons. org/ licen ses/ by/4. 0/.