Reprocessing and interpretation of legacy seismic data using machine learning from the Granada Basin, Spain

The Granada Basin (Spain) is a Neogene sedimentary depression with irregular geomorphology and deep depocenters. It is located in the most seismically hazardous part of the Iberian Peninsula with an historically experienced extremely destructive earthquakes, followed by periods of low to moderate seismicity. In 1980s the Chevron Oil Company collected a set of 30 deep seismic reflection sections in this Basin of which only the results on paper are kept. Due to the fact that many of these seismic profiles are currently located in urban areas and the economic cost of carrying out a similar exploration, it was decided to recover these old data and apply a post-stack treatment to improve their quality. The purpose of this study is to show the applied reprocessing flow and, with the new sections, to present a spatial model of the basin. The first stage of recovery and enhacement of seismic sections has consisted in three phases: first, high-resolution scanning of paper copies to TIFF images followed by the transformation of TIFF images to SEG-Y format; second, poststack processing workflow to increasing resolution and lateral coherence of these seismic lines; and third, it has been used a machine learning algorithm, among others, increasing the spatial resolution, signal-to-noise ratio, and coherence of the seismic signals. In addition, basement horizons, as well as three sedimentary sequences, were identified in all seismic sections and interpolated to create a three-dimensional basement model composed by normal faults, horst and grabens related to the seismotectonic behavior of the basin. As an overall assessment, this work is an example of the usefulness of ‘recycling ’ legacy seismic data, which nowadays are usually in archived boxes, but at the time required a great economic and acquisition effort.


Geological Background
The Granada Basin (GB) is a Neogene-Quaternary intramountain basin located in the south of the Iberian Peninsula (Spain); within the Betic Cordillera and bordered by the External Betic Zone to the northnorthwest and the Internal Betic Zone to the south-southeast.This basin (Fig. 1) has been accumulating sediments since the Neogene and is divided by normal fault systems that create distinct depocenters such as Granada, Cubillas, Pinos Puente-Fuente Vaqueros, Chimeneas, and Huetor Tajar (Morales et al., 1990;Sanz de Galdeano et al., 2003;Rodríguez-Fernández and Sanz de Galdeano, 2006;Agosta et al., 2012;Sanz de Galdeano et al., 2012;Azañón et al., 2013).In addition, the presence of the Mesozoic and Tertiary sediments in the External Betic Zone and the Paleozoic-PermoTriassic metamorphic alpine complex in the Internal Betic Zone contribute to the geological complexity of this region.
Tectonically, the GB is perched atop a low-angle extensional system, with ongoing extension influenced by sublithospheric processes such as slab retreat, partial detachment, and tearing of the Gibraltar slab (Bezada et al., 2013;Fichtner and Villaseñor, 2015;Gütscher et al., 2012;Heit et al., 2017;Spakman et al., 2018;Capella et al., 2020).This dynamic setting contributes to a predominantly normal faulting regime that shapes the basin's landscape, evident in the old seismic reflection lines and surface geology (Morales et al., 1990).Notably, fault activity during the Upper Pleistocene to Holocene period has linking specific structures such as the Atarfe, Pinos Puente, Padul, and Ventas de Zafarraya faults to seismic events, including the significant 1884 Andalusian earthquake (Reicherter et al., 2003), or most recently the seismic swarm in the Atarfe-Santa Fe sector between January and May 2021 (Lozano et al., 2022;Madarieta-Txurruka et al., 2022), as shown their location with red dots in Fig. 1.The Iberian Peninsula's highest seismic activity is recorded in the GB, mostly comprising microseismic events, and earthquakes rarely exceeding Mw 5.5 with depths around 20 km.To the the south-southwest there is an increase of the hypocentral depth up to 40 km (Morales et al., 1997;Galindo-Zaldívar et al., 1999;Carmona Rodríguez, 2009;Stich et al., 2020;Madarieta-Txurruka et al., 2021;Madarieta-Txurruka et al., 2022).

Study proposal
Between 1984 and 1986 the Chevron Oil and the General Geophysics S.A. companies carried out a geophysical survey in the GB (Fig. 2A), which included a gravimetric map of 250 ha, 30 seismic reflection lines with a total lenght of 412 km and two research boreholes, Granada-D1 (850 m depth) and .At present, this geophysical data set is accessible in hard-copy at the Andalusian Institute of Geophysics (IAG, University of Granada) and the Geological and Mining Institute of Spain (IGME, http://info.igme.es/SIGEOF/).The geophysical dataset may contain non-reproducible information from the study area due to the challenges posed by urban changes, high costs, and environmental restrictions, which currently make it difficult to replicate the experiment.
Within this framework, we have developed a methodology to convert these paper seismic sections into digital files that can be accessed and processed digitally.In this processing flow, we have incorporated the conventional post-stack processes, as well as the previous studies conducted by Cicala et al. (2021) and Brancatelli et al. (2022).Given this, the two main questions of our study are: How effectively can this methodology could be applied to recicled this legacy seismic data?; and Do these new seismic images provide new geological insights into the Granada Basin?.Finally, we hope that this approach can be extended to other seismically hazardous areas worldwide.

Old geophysical data set
Fig. 2 and Table 1 sumarise the main information form the geophisical set used in this study.The logs of two boreholes have also been scanned and digitized to be taken into account in the second phase of interpretation (Fig. 2B), even though our main objective is the recovery of the 30 seismic sections.Fig. 2C is a representative seismic reflection section-S-85154-that corresponds to the yellow line over the Sector 1 in the local map of Fig. 2A.The paper header in the left side of this section contains all the information we have been able to retrieve from the acquisition data and the old processing.Similar information goes for all paper sections.On this basis, Table 1 containt the most significant acqusition parameters and in the following section, we describe the historical processing.
Explosive charges located between 2 and 18 m depth were used as a seismic source, with charges (aluminium powder) varying between 1.45 g and 65.25 g.The sepread shooting was symmetrical with a lateral offset of 5 m.Eeach individual trace was made up of the sum of 18 geophones arranged in a circle or line, in order to eliminate random noise and the ground roll and amplify coherent signals.The field bandpass filter applied was, generally, 25 Hz (12 db/oct)-250 Hz (72 db/oct) and the original file format was SEGB.

Historical processing to obtain the seismic sections
The old seismic processing followed the classic multi channel seismic reflection flow (e.g.Yilmaz, 2001).Initially, the records underwent demultiplexing in order to disentangle the intertwined data streams.Next, during the pre-stack stage, a trace-amplitude correction recovery was implemented, which was then followed by the trace-editing step involving the kill of bad traces and the correction of polarities.Following a deconvolution process with a 120 ms operator length, a 200-2000 ms window, and 20% prewiting, the vertical (temporal) resolution is enhanced.The static corrections were applied to adjust elevations and mitigate the effects of weathering layer.These corrections were made using a velocity of 2500 m/s to minimize all traces at the datum plane (lowest elevation).The velocity analysis involves using constant velocity stacks (CVS) with 12 velocity ranges.Once the velocity laws are established for the entire line, dynamic corrections (NMO) are applied to all common depth point (CDP) gathers.Top-surgical muting was also implemented to remove non-reflective persistent events.The burte stack was achieved with 1200% to 1800% of CDP fold.During the post-stack stage, the wave equation migration (WEMIG) method was employed for accurate spatial relocation of seismic events.This was followed by the application of a variable time band-pass filter, with a range of 1-1200 ms for frequencies between 20 and 85 Hz, and a range of 1800-3000 ms for frequencies between 20 and 65 Hz.The ultimate improvements were made to the signal involved enhancing spatial coherence and implementing trace equalization to achieve consistent amplitude throughout the section.

Methodology used to recover the old seismic sections
The recovery of old seismic sections consist in three big steps which are described below.

Step I: Scanning and digitization of old seismic sections (paper to SEG-Y)
The seismic sections were originally printed in post-stack mode, in paper sheets, with a time window of 4 s (two-way travel time) with traces sorted by CDP.These paper-sheet seismic sections were approximately 2 × 1 m in size.The lateral composite of each seismic section shows general information about how the field data was collected, its geometry, and the processing settings.For example, the time sampling rate is spanning from 1 ms to 2 ms (Fig. 2C).Initially, the original sections included information about Lambert coordinates centered on the Madrid meridian, which were then transformed into ETRS89 zone 30 N geographic coordinates.
The high-resolution scanning procedure was performed using a plotter that was particularly adapted to the paper size.This resulted in full color TIFF images with a maximum resolution of 100 pixels per inch (ppi), which was chosen to maintain a reasonable file size managing (Fig. 2C).The TIFF images underwent a conversion to black-and-white (BW) scale, as well as cropping and resizing to meet the specifications of the SEG-Y Revision 2 transformation format code created by Farran (2012).This involved ensuring that the number of vertical pixels matched the number of samples per trace, and the number of horizontal pixels corresponded to the number of seismic traces that will be correlated with a navigation ASCII file to ensure geographic referencing in a further step.Ultimately, this technique resulted in the production of 30 BW images.
The process of converting TIFF images to SEG-Y format involved the use of two open codes.The initial action required the installation of the general SEG-Y libraries created by Hansen (2017) in the MATLAB environment.These libraries enabled the management of the resulting trace-header files.Then, for the direct conversion of each individual image, Farran (2012) developed the second code.
For the above conversions, the number of traces and samples per trace, as well as the ASCII file with the UTM coordinates for each trace in the ETRS89 zone 30 N, were considered as parameters.This enabled the georeferencing of every pixel along the profile, obtaning a SEG-Y file as shown Fig. 3A.After the data were successfully converted to the SEG-Y format (Araque-Pérez, 2023), it was observed that the first linear time marks were drawn every 100 ms, appearing as artificial reflections, or artifacts (as seen in Fig. 3B and C).

Step II: Post-stack seismic processing
The actual reprocessing flow was planned to be executed sequentially.Firstly, processing has been aimed at increasing the quality of the converted sections (Fig. 4A).To this end, GLOBE Claritas-Seismic Data Processing Software (petrosys.com.au) was employed.The following actions were taken (illustrative additional outcomes of each flow step are presented in the Supplementary Material).1.All seismic sections underwent a regularization process to mitigate the adverse effects of topographic corrections.Thus, the last 800 ms were removed to ensure that they did not influence subsequent interpretations.2. According with the resolution criteria, signal resampling was carried out from 1 ms to 2 ms to reduce the seismic section file size from 100 MB to approximately 50 MB.This was accomplished with the aim of an anti-aliasing filter, which was set up with a taper length equivalent to 20% of the cutoff frequency, established at 100% of the Nyquist frequency (250 Hz).Additionally, a zero phase adjustment was applied to prevent any changes in signal phase, thereby preserving the original waveform.3. Top muting was maked to eliminate spurious signals present above the reference level (topographic relief in this case).4. A predictive deconvolution step (Yilmaz, 2001) was used to minimize the remained timelines at 100 ms intervals.An input was a minimum phase Ricker wavelet (Ricker, 1953) with a lag size of 44 ms, which was applied to a time window of 3200 ms (with a 10% prewiting).This corresponds to 2.5 cycles of the original waveform, with a gap size of 100 ms (Gholamy and Kreinovich, 2014).This procedure successfully decreased the timelines while ensuring the preservation of data integrity.5. To improve the lateral coherence, a F-X deconvolution (5 traces of filter length, trace overlap of 20%, window length of 3200 ms and overlapped at 500 ms).
6.A time-variant band-pass filtering (Table 2) was utilized to reduce the remaining temporal noise.

Step III: Machine Learning implementation
The latest step has been aimed to improve the quality of the seismic sections as a general way.That is, to increase the signal-to-noise ratio and the coherence between traces, both in the vertical and horizontal direction.For this purpose, it has been considered to use the ML algorithms contained in the OpendTectdGB Earth Sciences software (V.7.0.4; dGB Earth Sciences -Tags -OpendTect, dgbes.com).As in the previos section, illustrative additional outcomes are shown on the Supplementary Material.
1.For all sections: a. UTM coordinates (ETRS89, zone 30 N) are encoded in trace headers.b 1 .A kill-trace prcedure, using GLOBE Claritas, over any poor signal data was make (e.g., Fig. 5A and B in wiggle mode and seismic color scale, respectivelly).The goal was the elimination of traces that had amplitude lower than half of the amplitude shown by adjacent traces.This sections set is defined as input data for ML implementation.
b 2 .Automatic Gain Control (AGC) was applied to to equalize the amplitudes of all seismic traces in the profiles obtained as outcome in the previous segment of the workflow (Fig. 6A and B in wiggle and seismic color, respectively).This sections set is defined as target data for ML implementation for the supervised deep learning tool.

Effective ML implentation tool (Groot et al., 2022).
Based on Unet architecture (Ronneberger et al., 2015), to interpolate traces with a low signal-to-noise ratio and coherency enhancing and replacing bad data sectors.The training ML model, as shown the workflow on Fig. 7, fills in missing seismic data from the input layer, which is trained and reconstructed as a new reconstructed image in the decoding part.For the application of the ML algorithm, it was necessary to establish the following parts of the neural network: a.The input layer (e.g., Fig. 5), which was segmented into 7584 blocks of 64 × 64 samples, each one overlapped 50% with its neighbors in the horizontal and vertical directions.The U-Net is a Keras-TensorFlow model with an input shape of 64 × 64 samples that are trained for 30 epochs (iterations) and a batch-size of 4, using a mean absolute error (MAE) as metric, a mean squeared error (MSE) as loss function, and "Adam" as optimizer.
The division of the input target into smaller sample blocks improves computation efficiency and makes data blocks more manageable in the encoder section of the U-Net.In this way, the encoder extracts a set of features such as the amplitude, frequency, and phase attributes from each fragmented block; making a seismic interpolation based on neural networks in the decoder section.So, the missing seismic traces are treated as a kind of "mask" in the reconstruction of the image where the amplitude of the loss signal are calculate updating the neural network trough the overlapping existing traces.
b.The target layer (Fig. 6) is convolved with the input segmented layer in labels, creating the training model, as shown the charts on Fig. 7.During this process, the MSE and MAE are represented in Fig. 8A  and B respectively.The low MAE value of 0.12 suggest that the signal was successfully reconstructed, and the decreasing trend in the learning curve (Fig. 8C) indicates that the deep learning algorithm was able to reconstruct the final seismic image with a low loss as validation (Fig. 8D).Furthermore, the above approach was applied to the remaining 29 seismic sections, which improved all previous seismic images overall.This has led to new interpretations, from which a Basement model has been obtained c.The results are reconstructed in an output layer (Fig. 9A for wiggle seismic image and Fig. 9B for a colored one), which serves as the final objective of ML procedure.The missing seismic values are computed using the input layer, which undergoes training and reconstruction   to generate a new seismic image in the decoding part of the ML algorithm.
The Groot et al. (2022) algorithm was selected because its convolutional neural network functions, working as an auto-encoder, enables an image-to-image regression that could lead to new seismic sections with higher data quality.Furthermore, the above approach was applied to the remaining 29 seismic sections, which improved all the legacy seismic images overall.This has led to new interpretations, from which a Basement model has been obtained.

Seismic interpretation
Three visible horizons as main seismic reflections are interpretated in these new sections.These horizons were then correlated with the sedimentary units found in borehole logs of Granada D1 (836 m depth) and Malaha 82 (482 m depth).The well-tie relation of the sedimentary units was made in the nearby seismic section S-85154, trougth the natural gamma-ray logs (Fig. 10A), and for the time-depth conversion the sonic logs (seen in Fig. 2B) were used to extrapolate a 2D-interval velocity model for the seismic section (Fig. 10B), which has been extended to the remaining profiles.
Preceding studies refer the existence of three sedimentary sequences (Morales et al., 1990;Rodríguez-Fernández and Sanz de Galdeano, 2006): (S3) Is a shallow a sequence characterized by conglomerates, sands, and clays dating back to the late Pliocene; (S2) is composed of evaporite sediments and micaceous sands from the late Miocene; and (S1) encompassing Tortonian sea flood plains within calcarenites.The base of S1 correspond to the basement of the GB.The basement top was interpolated trough all seismc profiles using the krigging method (Oliver and Webster, 1990) and subsequentlyit is related with the topographical relief.This process led to the creation of a three-dimensional basement model (2.5D interpolated surface), depicted in Fig. 12.
The simplified seismic processing workflow steps are illustrated in Fig. 11, specifically for the S-84151 seismic section located in the Cubillas depocenter.This figure sumarises all steps involved in the line restoration (conversion paper -SEGY, post-stack reprocessing, ML enhancement), and the final post-processing interpretation.

Reconstruction of gravimetric map
The geophysical campaign conducted by Chevron Oil included Bouguer anomaly maps and gravity field data for the GB region.The Bouguer anomaly maps, along with the original seismic sections, were reproduced on paper.An illustration of the original Residual Gravimetric Map can be found in the Supplementary Material.The maps were digitized utilizing the same scanner employed for acquiring TIFF images of the seismic sections.The TIFF map images were first georeferenced to the old Spanish Lambert coordinate system centered on the Madrid Meridian, and then transformed to the ETRS89 Zone 30 N system.Next, the isolines that corresponded to the residual Bouguer anomalies were converted into digital format, resulting in an ASCII file containing XYZ coordinates sorted in columns.In this file, the Z value represented the residual Bouguer anomaly, which could be accurately interpolated by kriging method (Oliver and Webster, 1990;Hansen, 1993) to generate a digital version suitable for an accurate interpretation.

Construction of the basement model
The basement model exhibits five known depressions, namely the Chimenas (CH), Pinos Puente-Fuente Vaqueros (PP-FV), Cubillas (CU), Granada (GR), and Huetor Tajar (HT), in Fig. 12.Among these, the Pinos Puente-Fuente Vaqueros and Cubillas depressions stands out as the deepest, reaching a depth of 1900 m and 1700 m respectively.Fig. 12A shows the general interpolation between the seismic profiles for the basement surface 3D view.The seismic profiles do not completely cover the Huetor Tajar depression, showing a portion of the depocenter that is not certain to have reached its deepest point.The newly proposed basement model (Fig. 12B) has been contrasted with the digitized Bouguer anomaly map depicted in Fig. 12C, see Table 4.This comparison was conducted in order to verify the reliability of the new basement depths, particularly in relation to the Cubillas' depression.
In relation to the Bouguer anomalies map, the depocenters of the GB exhibit residual Bouguer anomalies within the range of − 4 to − 10 mGals.These lower gravimetric values coincide with the depressions observed in the seismic basement (Fig. 12B and C).The Huetor Tajar depression, located on the western flank of the GB, presents the lowest gravimetric value of − 10 mGals, but this depression has not been completely covered by the seismic profiles, leaving uncertainty in the maximum depth reached.Following this, the Cubillas' depression, located in the northeastern region, displays a gravimetric value of − 8 mGals.Notably, the Cubillas' depression is separated from the Pinos Puente depocenter by the Sierra Elvira uplift, where the anomalies exceed 7 mGals.The Chimeneas and Granada depressions have a value of − 8 and − 5 mGals respectively (Table 4).

Cubillas depocenter analysis
The Cubillas basement is shown as an exemple of the detailled analysis that could be done in order to get new geological insigths.Fig. 13 summarizes graphically the entire procedure followed to study this sector.First (Fig. 13A), a detailed interpretation of the sedimentary sequences and velocities are considered for all implied seismic sections.Second, the results in 2D and 3D views are anlysed (Fig. 13B and C), obtaining a strong gradient similar to the Puente-Fuente Vaqueros depocenter.

New contributions
According to previous studies (Morales et al., 1990;Rodríguez-Fernández and Sanz de Galdeano, 2006), the Cubillas depocenter is situated at a depth of 1400 m, but in our new basement model this depocenter is placed at 1700 m; indicating a greater depth of 300 m.Furthermore, in the other depocenters examined, there is also a discrepancy of around 100 m in depth compared to the previous hypotheses.As a result, the increased sediment accumulation within the Cubillas and the others depocenters, the volume of sediment and sedimentation rate in this basin could surpass previous estimates.This surmise aligns with the isostatic compensation elevations of the Sierra Nevada and Alpujarride complexes, which are known to directly influence the subsidence of the depocenters along the eastern margin of the GB (Hürtgen et al., 2013;Madarieta-Txurruka et al., 2023).Nonetheless, sedimentation rate calculations has not explicitly performed in this study, but in future studies could address this gap to understanding of sediment dynamics in this region.
Rodríguez-Fernández and Sanz de Galdeano ( 2006) used an average velocity of 1.9 km/s for converting the time-section to depth-section, whereas we have used a variable velocity model with the stratification.This also implies differences in the geological characteristics assessed.This new contribution does not imply that previous procedures were incorrect; rather, it highlights that they were based on different approaches.In other hand, the velocities of the basement layers in our model correspond with the high-velocity anomalies identified in the Mesozoic basement by Serrano et al. (2002), indicating a consolidated structural framework for these sedimentary strata with highger velocities than typically observed in Neogene sedimentary rocks.The high basement velocities observed in our model also imply significant seismic impedance contrasts according with Gil-Zepeda (2002) observations.The velocity model for the entire GB is summarized in Table 3.
Considering the sedimentary composition of the GB, which includes detrital materials, sandy sediments, silts, calcareous marls, and limestone from top to bottom distributed in three sedimentary sequences from the Late Miocene, and the velocities ranging from 2400 m/s to 3300 m/s.If we establish a theoretical densities between 2.3 g/cm 3 and 2.5 g/cm 3 for these kind of sediments, it is feasible to calculate a basement density contrast within these values.Assuming that the sequence with the densest material (S1) is found at greater depth, and has greater compaction, is suggested a density variation of approximately 0.1 to 0.2 g/cm 3 .Using the equation g = 2πGΔρh, where g is the observed gravity in mGals, G is the universal gravitational constant (6.674 × 10 − 11 m 3 kg − 1 s − 2 ), and Δρ is the density contrast, the depth  4 for the maximum anomaly values at each depocenter.
From the previous observations, it can be deduced that for the Chimeneas (CH) and Huetor Tajar (HT) depocenters, the density contrast values that yield depths similar to those observed in this study, and even to those previously presented, are approximately 0.2 g/cm 3 .Conversely, for the Cubillas (CU) and Granada (GR) depocenters, these density contrasts are closer to 0.1 g/cm 3 .This variation may be attributed to the influence of isostatic uplift of the Nevado-Filabride and Alpujarride Complexes interacting with the Granada, El Fargue, and Alfacar faults (Madarieta-Txurruka et al., 2023).Such geological interactions, with the density contrast assumptions, validate the basement depth data obtained in this work.
The only depth basement depocenter data that lacks self-explanation within these framework is the basement depth in the Pinos Puente-Fuente Vaqueros (PP-FV) depocenter, where none of the density contrasts are enough to reach the depth of 1900 m, as shown the Fig. 12C for the residual Bouguer gravity anomaly of − 4 mGals.This suggests an unusual scenario where the overlying sediments are denser than the underlying basement.The resulting density contrast is approximately − 50.2 kg/m 3 or − 0.0502 g/cm 3 , which deviates from typical geological expectations where the basement rocks are denser than the overlying sedimentary layers.
The above could be explained by several geological factors occurred during the Messinian Salinity Crisis and the unique hydrogeological settings of the basin described as follows: The combination of these factors in the Pinos Puente-Fuente Vaqueros depocenter could contribute to the observed negative density contrast, specially when this basin was affected by a Salinity Crisis during the Later Miocene, as well other basins around the world (Deniaud et al., 1999;Carpentier et al., 2020;Duggen et al., 2003), forming evaporitic deposits in the GB (salt, gypsum).This scenario underscores the complexity of interpreting gravimetric data in geologically active and structurally complex regions like the GB, that are validated by the seismic profiles in the area which are clarifying a depth basement in this sector.
Furthermore, a hypothesis posits that the initial formation of the GB may have occurred within the Pinos Puente-Fuente Vaqueros depression, then progressing along a northeastern trajectory towards Cubillas (Morales et al., 1990).Based on the findings of this study, it is possible to  2012), it is important to note that the uplift of the Sierra Elvira Block served as a obstruction in this process.Hence, further research is required in order to address this inquiry.

Discussions
This manuscript presents two actions.The initial approach has been characterized by a methodical process, demonstrating the sequential actions taken to restore and update old seismic sections through the utilization of current signal processing techniques.In the second case, geological findings have been achieved, such as the analysis of the Cubillas depocenter, which revealed a significant difference of 300 m in its depth compared to previous studies.In contrast, the other depocenters only differs in 100 m with the previous desccriptions of their maximum depth.
The most novel aspect of the methodological part has been the use of the ML algorithm proposed by Groot et al. (2022) to restore seismic traces from sections with low signal content.In our case, the low MAE and MSE values (0.12 and 0.04, respectively) obtained indicate the remarkable improvement of the seismic image reconstruction.
It should be mentioned that this ML altorithm differs from classical algorithms for improving the lateral and spatial coherence of seismic sections because the latter operate by modifying amplitudes and dips between neighbouring traces, while the ML is fixed on the whole blockpattern.This difference ensures less artifacts creation, while the presence of bad traces (or null traces) produces artifacts in classical algorithms.
The geological results obtained in the second part of this study have shown the success of the rescue seismic sections, because it it has been possible to identify sedimentary sequences and the underlying structures providing new geological information of the GB.This new obtained data, in conjunction with the previous studies, have been significant.The 3D reinterpretation of these sedimentary horizons facilitates leading to a deeper understanding of subsurface structures, and the model can also be used in seismic hazard modeling of the basin.
One limitation of seismic secctions employed was that the profiles were constrained to a two-dimensional (2D) mode.This constraint restricts model resolution, because the calculation of the extented basement and others horizons can only be done using interpolation methods leading to a 2.5D representation.

Conclusions
This study was designed to assess the feasibility of retrieving and utilizing old geophysical data.Ultimately, the validation of this work would be deemed positive if these data, appropriately processed using current techniques, yield new insights into the area where they were collected compared to their initial processing stage.If this "recycling" proves to be successful, this manuscript provides an illustrative example of this course of action.
Through this manuscript we have shown that the recovery of legacy seismic sections of the Granada Basin and their subsequent use and interpretation has been positive in obtaining new geological information about the complex Granada Basin.
It should be noted that the methodology described can be replicated using any seismic reflection processing software, such as free available Seismic Un*x code.Furthermore, it could be integrated and executed by     any machine learning tool utilizing the U-Net architecture coded in Python, as well as the one developed by Groot et al. (2022), to overall enhance the seismic data quality.
The authors of this study are affiliated with the Andalusian Institute of Geophysics of the University of Granada (https://iagpds.ugr.es/instituto/presentacion/historia; IAG-UGR), which was established in 1902 and comprised three divisions: astronomical, seismological, and meteorological.Over time, additional areas of study such as volcanology and geophysical exploration were included.The extensive duration of these studies results in a significant accumulation of data in the Earth

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Geologic map of the Granada Basin modified from Araque-Perez (2024).The ellipses highlight the five GB's depocenters: Huetor Tajar, Chimeneas, Pinos Puente -Fuente Vaqueros, Cubillas, and Granada.Blue lines: Main active faults (QAFI, 2021).Red dots: historical earthquakes.The bottom right black box depicts the location of the Granada Basin on the Iberian Peninsula.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. (A) Location of 30 seismic profiles and 2 boreholes on a digital terrain model (UTM ETRS89 zone 30 N System).The cyan rectangles indicate the three zones with the densest seismic profiles distribution: 1-Pinos Puente and Chimeneas, 2-Granada, 3-Cubillas.The green filled polygon represents the location of the Granada city.(B) Lithology and digitized logs from Borehole Granada D1 with the markers identified for the sedimentary sequences.(C) Original paper seismic section (S-85154) and its location in Zone 1 through a yellow arrow.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. (A) Seismic Section S-85154 converted to SEG-Y format.The highlighted rectangular region denotes an enlarged area of interest.(B) Within the enlarged portion, certain flat artifacts, resulting from the original division of timelines at 100 ms intervals.(C) Same image as (B), changing the color scale to emphasize timelines that produce flat artifacts.The arrows indicates the existence of timeline marks (artifacts) within the seismic section.

Fig. 4 .
Fig. 4. (A) Example of a seismic section (S-85154) after applying the first processing flow stage.(B) Zoomed-in area showing the suppressed biased horizons, as is compared with the same counterpart of unprocessed profile of Fig. 3B.

Fig. 5 .
Fig. 5. (A) Seismic section (S-85154) in wiggle mode with blanked traces (input image) before applying the ML algorithm.(B) The same image as (A) but with a seismic color scale.

Fig. 6 .
Fig. 6. (A) Seismic section (S-85154) in wiggle mode after the first stage workflow, adding an AGC procedure to be set as a target image for ML.(B) The same image as (A) but with a seismic color scale.

Fig. 7 .
Fig. 7. Machine Learning implementation with a U-net architecture.The encoder divides the input image (on the up-left side of the U-net scheme) into smaller 64 × 64 sample blocks, which convolve the symmetric units of the target image (on the down-left side of the U-net scheme).The decoding component reconstructs a new image (right side of the U-net scheme).The dotted red frame highlights the architecture of the U-net ML architecture schematic depicted at the bottom of the ML workflow.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

1.
Hydrothermal Alterations: The presence of CaCl 2 -rich hydrothermal waters in the GB (García-Veigas et al., 2013), might have chemically altered the basement's density.Such hydrothermal activities can lead to the leaching of denser minerals or the deposition of lighter materials, thereby decreasing the density of the basement relative to the overlying sediments.2. Evaporite Recycling: The recycling of evaporites, influenced by non-marine water inputs, is another process discussed by García-Veigas et al. (2013).This can change the mineral composition and density of the basement.The introduction of lighter minerals or the removal of denser components during the recycling of previously precipitated evaporites could result in a less dense basement.3. Tectonic and Structural Modifications: Significant tectonic activities during the Messinian Salinity Crisis could have disrupted the stratigraphic layers, possibly bringing less dense materials into the basement or altering the compaction state of the sediments.This is suggested by Riding et al. (1998) and could lead to a reversal in the expected density.

Fig. 8 .
Fig. 8. Training and validation resume of the ML tool in the S-85154 seismic section.(A) Epoch advance versus mean square error (loss, or typically MSE).(B) Epoch advance versus mean absolute error (MAE).(C) Graphic of epoch advance versus learning between the input and target data.(D) Number of iterations made it (by all bach sizes) versus loss.The blue line represents the learning model (trained model), and the orange line is the validated final model as reconstructed image.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. (A) Seismic section (S-85154) in wiggle mode after applying ML tools as a new reconstructed image.(B) The same image as (A) but with a seismic color scale.

Fig. 10 .
Fig. 10.(A) Seismic section S-85154 tied to gamma-ray data from Granada D-1 and Malaha 82 boreholes in order to convert seismic time to depth.The analysis identified three principal horizons in the seismic section.The uppermost horizon (S3) represents a near-surface sedimentary sequence from the late Pliocene.The middle horizon (S2) consists of sedimentary sequences of marine and alluvial facies from the late Miocene.The lowermost horizon (S1) indicates sea flood plain facies that date back to the Tortonian identified by (Rodríguez-Fernández and Sanz de Galdeano, 2006).The colored dashed lines in the borehole data indicate marker units that are correlated with seismic horizons in the seismic profile.(B) Velocity model derived from the well-tie procedure on S-85154 image.CH: Chimeneas depocenter; PP-FV: Pinos Puente-Fuente Vaqueros decpocenter.

Fig. 11 .
Fig. 11.(A) Digitized Raw Seismic line S-84151.The arrows shown the biased reflectors from original timelines marks.(B) Seismic line S-84151 after the first processing cycle.(C) Seismic line S-84151 after the second processing cycle with ML algorithm application.(D) Horizon selection, stratigraphic sequences, and fault interpretation after fully implemented data processing.

Fig. 12 .
Fig. 12. (A) A general spatial view of seismic profiles with the calculated basement model.(B) Basement model depth.The color scale represents the altitude above mean sea level, and the isolines reflect the depth below the topographic relief.(C) Residual Bouguer anomalies, overlayed on the digital elevation model (UTM ERTS89 zone 30 N coordinate system).The white square contains the study case of Cubilla depression.SE: Sierra Elvira.

Fig. 13 .
Fig. 13.(A) Seismic section S-84151 with normal faulting in the Cubillas depression.The stratigraphic sequences are on the velocity model color scale.(B) Map of the Cubillas depocenter with elevations mean sea level (UTM ETRS89 30 N) and active faults (QAFI, 2021).(C) 3D map of the Cubillas' basement.Stratigraphic sequences are illustrated in the seismic sections.

Table 1
Summary of acquisition parameters.

Table 3
Summary of velocities per sedimentary sequence within the GB.

Table 4
Description of Basement observations.The observed depths in the HT's depocenter are inconclusive due to the limited seismic coverage in the region, it is traversed by only one seismic section from 1984's survey (section S-84150 in Fig.2A).Sciences field, which requires a substantial investment of financial, human, and environmental resources.Consequently, by updating them, we not only acquire new insigths, but also acknowledge the contributions of our predecessors.