Submerged Marine Terraces Identification and an Approach for Numerical Modeling the Sequence Formation in the Bay of Biscay (Northeastern Iberian Peninsula)

Submerged sequences of marine terraces potentially provide crucial information of past sea-level positions. However, the distribution and characteristics of drowned marine terrace sequences are poorly known at a global scale. Using bathymetric data and novel mapping and modeling techniques, we studied a submerged sequence of marine terraces in the Bay of Biscay with the objective to identify the distribution and morphologies of submerged marine terraces and the timing and conditions that allowed their formation and preservation. To accomplish the objectives a high-resolution bathymetry (5 m) was analyzed using Geographic Information Systems and TerraceM®. The successive submerged terraces were identified using a Surface Classification Model, which linearly combines the slope and the roughness of the surface to extract fossil sea-cliffs and fossil rocky shore platforms. For that purpose, contour and hillshaded maps were also analyzed. Then, shoreline angles, a geomorphic marker located at the intersection between the fossil sea-cliff and platform, were mapped analyzing swath profiles perpendicular to the isobaths. Most of the submerged strandlines are irregularly preserved throughout the continental shelf. In summary, 12 submerged terraces with their shoreline angles between approximately: −13 m (T1), −30 and −32 m (T2), −34 and 41 m (T3), −44 and −47 m (T4), −49 and 53 m (T5), −55 and 58 m (T6), −59 and 62 m (T7), −65 and 67 m (T8), −68 and 70 m (T9), −74 and −77 m (T10), −83 and −86 m (T11) and −89 and 92 m (T12). Nevertheless, the ones showing the best lateral continuity and preservation in the central part of the shelf are T3, T4, T5, T7, T8, and T10. The age of the terraces has been estimated using a landscape evolution model. To simulate the formation and preservation of submerged terraces three different scenarios: (i) 20-0 ka; (ii) 128-0 ka; and (iii) 128-20 ka, were compared. The best scenario for terrace generation was between 128 and 20 Ka, where T3, T5, and T7 could have been formed.


INTRODUCTION
Marine terraces are geomorphic features of ancient shorelines, which constitute the geological record of past sea-level positions (Lajoie, 1986;Burbank and Anderson, 2011). About 1020 sequences of coastal indicators including the one related to the last interglacial maximum (Marine Isotopic Stage, MIS5e) have been reported in literature (compilations in Pedoja et al., 2014Pedoja et al., , 2018. Due to the lack of high resolution data, difficulties in carrying out direct observations and dating of underwater landforms, the submerged parts of the coastal sequences are totally neglected and studies dealing with submerged marine terraces are rather scarce (e.g. Passaro et al., 2011;Johnson et al., 2014;Jara-Muñoz et al., 2017;Ricchi et al., 2018). The analysis of submerged marine terraces and sea stacks, combined with sea level variations and tectonic settings, helps us to elucidate the Quaternary coastal evolution of the SE Bay of Biscay. Galparsoro et al. (2010) were the first to attempt the characterization of submerged features and identified eight submerged marine terraces in the SE Bay of Biscay at approximately: −37, −52, −56, −70, −73, −75, −87, and −92 m below sea level.
In this study we use a high resolution bathymetrical model of the Bay of Biscay down to depths of −116 m, to map submerged shoreline angles using morphometric analysis. Alongside, the detection of sea-stacks was carried out to better define the number and distribution of the submerged terraces. Due to the lack of dating, to constraint the chronology of the successive submerged marine terraces we used numerical modeling.

BACKGROUND
The study area is located on the SE Bay of Biscay of the continental shelf, N of Iberian Peninsula. It represents a coastal stretch estimated to be of ∼ 150 km (Galparsoro et al., 2010; Figure 1).

Geology of the Study Area
The substratum of the studied coastal stretch is mainly composed of Cretaceous and Tertiary flysch, and at some sites, urgonian limestones are also outcropping (Ábalos, 2016;Euskal Kostaldeko Geologia Sinplifikatua, 2017). The submerged sequence of marine terraces is located on a passive margin that experienced uplift during Upper Cretaceous -Miocene times, in the Alpine orogeny (García-Mondéjar et al., 1985). Along the Cantabric Margin, various uplift rates have been reported: 0.06 -0.15 mm/yr (Jiménez-Sánchez et al., 2006;Jimenez-Sánchez et al., 2011), 0.07 -0.15 mm/yr (Alvarez-Marrón et al., 2008) or 0.08 mm/yr (Aranburu et al., 2015). Regarding neotectonics, the Quaternary Active Faults Database of Iberia (QAFIv3) (IGME, 2015) does not report any fault in the area. Nevertheless, taking into account the neotectonic map (IGME and ENRESA, 1998), to the east of the study zone a fault is mapped and proposed to have been active during Pliocene or possibly Quaternary times (Figure 1).
The modern coastline is roughly E-W oriented. More than 70% of it is characterized by moderate to high sea-cliffs (20-150 m). Sandy beaches (predominantly of "pocket" type) make around 21 km (14%) (Pascual et al., 2004). Cape Matxitxako constitutes the main headland of this part of the Basque coast, where the continental shelf is the narrowest 7 km (Uriarte, 1998; Figure 1).

Oceanography of the SE Bay of Biscay
The southern coast of the Bay of Biscay is exposed to large storms from the NW, due to the North Atlantic low pressure systems (Galparsoro et al., 2010). Strong NW swell waves dominate and are the most common sea state within the study area (Galparsoro et al., 2010). Using the information of the Bilbao offshore buoy Liria et al. (2009) summarized that in summer the swell waves period is <10 s more than 75% of the time with representative heights of 1.5 m. In winter the swell waves have higher period (i.e., 13 s), and heights >2 m more than% 50 of the time. Several times a year, significant swell wave heights can exceed 5 m, and with return periods of 20 years the waves can reach 10 m.

Marine Terraces and Past Sea-Level Positions
Marine terraces are geomorphic markers that constitute the geological record of past sea level positions (Lajoie, 1986;Burbank and Anderson, 2011). These features are composed by a sub-planar surface that dips seaward, and are limited landward by a steep fossil sea-cliff (Lajoie, 1986; Figure 2). The intersection between the shore platform and the sea-cliff is called the shoreline angle, which represents the maximum height that reached the sea level in the marine terrace formation (Lajoie, 1986). Staircase sequences of coastal terraces, either marine, coral reef or sedimentary are present along several coasts, in any geodynamical context, due to the combination of sea-level variations and landmass uplift (Pedoja et al., 2011(Pedoja et al., , 2014. Rocky shores are dominated by erosive processes which leave morphological remnants such as arches, sea-stacks or active seacliff (Sunamura, 1992). More specifically, sea-stacks, are rocky columns that were originally part of the cliff, but that nowadays are separated from the mainland, creating an islet (Figure 2). The depth of the top of the sea-stack, if flat, can give an approximation of the depth of the immediately above terrace whereas its base lies on a former platform (Figure 2).

MATERIALS AND METHODS
We used a 5 m resolution bathymetric Digital Bathymetric Model (DBM) of the Basque continental shelf until a depth of −116 m (with respect to the mean sea-level in Alicante). The DBM was generated through a Multibeam Echosounder between 2005 and 2008, with a high-resolution SeaBat 8125 and SeaBat 7125 MBESs (RESON, 2002(RESON, , 2006Galparsoro et al., 2010). Most of the work was carried out using the SeaBat 7125 model; its operational frequency is 400 kHz, producing 256 beams in a 128 • angle swath and using up to 50 swaths per second. The beam width is 0.5 • along-track and 1 • across-track, producing very small  footprints; these, in turn, result in high horizontal resolution digital elevation models (DEM). The MBES was coupled with an Agp132 (TRIMBLE) global position system, receiving differential corrections. An OCTANS III (IXSEA) gyrocompass and motion sensor was utilized, to compensate for the movement of the vessel. Furthermore, a portable SVP 15 (RESON) was used, to measure sound velocity profiles throughout the entire water column (Ernstsen et al., 2006). The software package PDS2000 was used to integrate the MBES data, with the information from all the auxiliary sensors during the surveys -data acquisition and synchronization. This software was used in real-time, as well as in the post-processing of the integrated data. Tidal correction was applied using the nearest tide gauge and 1 m resolution seafloor DBM was produced in projected coordinate UTM, Zone 30 N (WGS84). The DBM was generalized into a 5 m grid, in order to increase the speed of computational processing; finally, it was exported into ESRI grid format and integrated into a Geographic Information System (Galparsoro et al., 2009).
The information relative to the seafloor type was obtained from the spatial data infrastructure of the Basque Government 1 . Different lithostratigraphic maps are available for the emerged part of our study area, with different scale and detail levels. The lithologies of the studied coastal stretch were provided

Morphometric Analysis
To detect the submerged fossil rocky-shore platforms and seacliffs, various methods were combined. First, we developed a Surface Classification Model (SCM) using TerraceM R (Jara-Muñoz et al., 2019; Figures 3A,B). The SCM linearly combines the slope and the roughness of the terrain to discriminate different geomorphic elements of submerged marine terraces. Using the seafloor type map from the spatial data infrastructure of the Basque Government, we selected from the DBM only the areas that laid on rocky and mixed (rocky/sedimentary) sea floors. Then, potential fossil sea-cliffs were isolated, by combining slope threshold values of 0.5 • , 1 • , and 2 • and roughness threshold values of 0.3, 0.6, and 0.9.
To identify the submerged sea-cliffs, we also extracted a contour map from the DBM ( Figure 3C) and the hillshaded map (angle of light incision set at 0 • , and a vertical exaggeration factor at the maximum available value: 99.9). Such hillshaded map, discriminates the sediment bodies from the rocky seafloor areas and highlights other features such as submerged channels ( Figure 3D).
To get the most accurate position and depth of the successive submerged shoreline angles, we used bathymetric swath profiles. We used a swath width of 250 m based on the resolution of the bathymetric datasets, orienting the profiles perpendicular to the trace of the fossil sea-cliffs ( Figure 3E). We used the maximum topography from swath profiles to determine the location of shoreline angles, considering that the maximum topography represent the original terrace morphology without the effects of erosion. For that purpose, we marked two points to define the fossil sea-cliff and another two points to define the fossil platform. Then lineal regressions are calculated upon the segments enclosed by the points, and extrapolated to find the intersection that marks the position of the shoreline angle ( Figure 3F). Vertical errors of shoreline angles are based on the extrapolation of the 2σ ranges of the linear regressions (Jara-Muñoz et al., 2019). The shoreline points lying on mainly sedimentary bodies were avoided, as they could underestimate their original height.
We used the original 5 m resolution DBM to map the shoreline angles of submarine terraces, to generate the hillshaded and contour maps. Instead we used a resampled version of 20 m/px in the detection and visualization of paleocliffs with the Surface Classification Model of TerraceM. The overall results are consistent after carrying out this resampling, as it can be observed comparing them with the hillshaded and contour maps.

Sea-Stack Analysis
Sea-stacks are relicts of eroded ancient marine terraces, and can be used to reconstruct the initial marine terrace morphology. The flat top of a given stack can indicate the depth of the terrace, above which the stack is emplaced. Such approach helps to detect the presence of submerged terraces in areas where nowadays their original morphology has been eroded and their shoreline angles are not detected. Thus, it could help to complete the sequence of terraces.
Sea-stacks are characterized by closed isobaths on the contour map, determinable with geographical information system (GIS). Then, we measured the range of the potential stacks from the bottom to the top, and we kept only the ones whose height is >2 m. Some of the closed isobaths were found on sedimentary zones, and consequently were not taken into account. Moreover, because some closed isobaths can be holes and depressions instead of sea stacks, the obtained close isobaths were also analyzed manually and the holes removed and discarded. Finally, the depth of the top of each stack was extracted from the bathymetry.
Such method to detect sea-stack through the analyze of closed isobaths has also its limitations. First, the real height of the stack could be underestimated in some occasions. Second, the detected feature can be a wide islet, without a columnar morphology. Third, the top of the sea-stack can be sharp instead of flat.

Modeling Marine Terraces Generation
Landscape evolution models (LEM) can be useful to estimate the age of marine terraces, especially when direct sampling and dating are difficult, as for submerged sequences of marine terraces. Herein, we used a landscape evolution model based on the wave energy dissipation theory proposed by Sunamura (1992). The model assumes wave erosion as a linear function of the rate of wave energy dissipation (Sunamura, 1992) that follows an exponential increase in relation to depth landward (Jara-Muñoz et al., 2019). The model allows simulating the generation and erosion of submerged terraces during a certain period of time, for different morphologic, climatic, and tectonic conditions.
For this study, a composite sea-level curve was generated combining various sea-level data. In this study we infer that marine terraces formed before MIS 5e have a lower probability to be preserved due to erosion from subsequent sea-level variations (e.g. Jara-Muñoz et al., 2016). Indeed, terraces before the MIS 5e may have been affected by strong wave erosion during the episode of sea-level rise that preceded the MIS 5e and former highstand periods before. For that reason, we used a composite sea-level curve truncated at 128 ka. Our composite sea-level curve for Upper Pleistocene (Figure 4) was obtained combining the sea level curves used by Arz et al. (2007), Stanford et al. (2011), and Grant et al. (2012). The sea-level curve of Stanford et al. (2011) is itself a combination of previously published sealevel curves (Lighty et al., 1982;Shinn et al., 1982;Robbin, 1984;Macintyre et al., 1985Macintyre et al., , 1995Macintyre et al., , 2004Digerfeldt and Hendry, 1987;Fairbanks, 1989;Chappell and Polach, 1991;Bard et al., 1996Bard et al., , 2010Hanebuth et al., 2000Hanebuth et al., , 2009Yokoyama et al., 2000Yokoyama et al., , 2001Cutler et al., 2003;Toscano and Lundberg, 2003;Peltier and Fairbanks, 2006).
We performed several sensitivity tests to find the best scenarios that may favor the development and preservation of submerged marine terraces (Figure 4): • Scenario 1: The submerged terraces were formed between 20 ka and present (the yellow line in Figure 4). • Scenario 2: The submerged terraces were formed between 128 and 20 ka (the green line in Figure 4). • Scenario 3: The submerged terraces were formed between 128 ka and present (the red line in Figure 4).
The wave height of the study area was estimated from 9 years of radar altimetry measurements (produced and distributed by the AVISO website of the French Spatial Agency, CNES, at www.aviso.com as part of the Ssalto/Duacs ground-processing segment). Different latitude/longitude positions were tested to extract the wave height. We selected wave height values between 3.7 and 3.4 m which are the 90% percentile of the data obtained in our study area (Supplementary Figure S1), assuming that these waves would be the most erosive ones related to storm events (Jara-Muñoz et al., 2016).
To constraint uplift rates in this area we used the database of Pedoja et al. (2018) based on the MIS5e benchmark widely distributed along the western coast of Europe. Estimates in the proximities to our study area, show uplift rates of 0.049 m/ka (Supplementary Figure S1), in Castro Urdiales (Guilcher, 1972;Flor and Flor-Blanco, 2014) and in Pointe Saint Martin (Biarritz) (Deler, 1932).
The initial slope of the model was constrained by measurements along predefined profiles in the continental shelf. For that, two points were selected on the profile, and using the distance and elevation differences, the inclination was calculated. The slopes from the bottom to the top of the profile were lower than 2 • in most cases. For modeling values of 0.5 • , 1 • , 1.5 • , and 2 • were tested.
We constraint the initial erosion rates based on previous studies close to our study area. To the West of our study area, in the coast of Lapurdi, mean coastal retreat ranges from 0.05 to 0.5 m/yr, but are typically of 0.15 m/yr (Aubié et al., 2011). To the East, in Galicia, Perez-Alberti et al. (2013) reported rates between 0.09 and 1.18 m/yr (typically between 0.21-0.38 m/yr) in granites. Even though in the model the input data is not directly the cliff retreat, for the initial erosion rates, values of 0.25, 0.5, 0.75, and 1 m/yr were tested. In scenario 1 (i.e., sequence formed during 0-20 ka), the initial erosion rates tested for 2 • of slope, were only 0.25 and 0.5 m/yr. Frontiers in Earth Science | www.frontiersin.org FIGURE 4 | The generated sea-level curve, and the tested three scenarios. MIS, marine isotopic stages from Railsback et al. (2015).

The Submerged Marine Terrace Sequence of the SE Bay of Biscay
The Surface Classification Model (SCM), the hillshaded map and the isobath map reveal a staircase morphology of the continental shelf, with irregularly distributed and preserved submerged marine terraces. In advance, the submerged marine terraces mapped are the ones that present a better preserved and recognizable fossil sea-cliffs. Submerged terraces are recognized from the eastern part of the study area until Barrika to the West, but the best preserved sequence is observed in the central part, especially in front of Ogeia (Figures 5A,B). The low general slope of the continental shelf has to be taken in mind, and note, that even the fossil sea-cliffs are highlighted clearly in some areas (e.g. Figure 3C), the slope of some of those cliffs can be of very low degree.
The detected submerged terraces distribution along the continental shelf and the main characteristics of the stretches that we defined are summarized on Tables  There are two differentiated groups: one between ∼ −32 and ∼ −17 m with a higher amount of stacks (more than 50 in some cases) and another between ∼ −32 and ∼ −95 m, for which there are less stacks (mainly less than 30). In this last group, the number of peaks is especially high at three depths: −54, −75/−76 and −90 m (Figure 6).
Qualitatively, in some zones there is a concentration of seastacks whose top coincide at a similar depth, and that seem to indicate eroded ancient platforms or stands in sea-level. Some of those potential levels are highlighted in Figures 7, 8. The rest of the areas are available in Supplementary Figures S10-S15.

Modeling Submerged Marine Terrace Sequence Generation
The best time span for submerged marine terrace creation is between the Last Interglacial Maximum and the Last Glacial Maximum (LGM). In this scenario, well developed submerged terraces are generated at ∼ −6, ∼ −36, ∼ −52, ∼ −68 and ∼ −118 m. Other lower magnitude shoreline angles are detected at ∼ −80 and ∼ −97 m ( Figure 9A). The profile generated between Frontiers in Earth Science | www.frontiersin.org TABLE 1 | Submerged marine terrace distribution from Jaizkibel to Ogeia in segments of the continental shelf with the most remarkable characteristics.

Continental shelf location
Location in Figure 5A Detected submerged terraces and shoreline angles depths MIS5e until nowadays, shows a smoothed and curved profile, with very low magnitude cliffs. Along this profile, the inflection points at ∼ −36 and ∼ −68 m are remarkable ( Figure 9B). The large morphological differences between both scenarios are due to the effect of the last period of sea-level rise after the LGM. Thus, the last transgression would have reworked the former coastal landscape and eroded older geomorphic evidences of past sea-level positions, giving as a result an smoothed surface.
The scenario related to the last transgression, is characterized by a profile with poorly marked cliffs where the main shorelines are located at ∼ −92 and ∼ −104 m ( Figure 9C). In this scenario the maximum initial erosion rate considered for modeling was of 0.5 m/yr.

Processes Implied in the Submerged Marine Terraces Distribution
In the attempt to elucidate the factors that control the location and spatial variability of the submerged marine terraces, we analyzed the effects of wave climatology, lithology and stratification direction. Nevertheless, didn't observe clear control from a single factor but what we suppose is the combination of these variables.

Wave Climatology
Wave direction and energy, can affect to the formation and development of rocky shore platforms (e.g. Trenhaile, 1999). The data available from the Bilbao offshore buoy (Boya de Bilbao-Vizcaya, 2019), record main wave directions of N292.5 • E and N315 • E ( Figure 10A). Considering the present coastline direction, the marine terrace sequence is best preserved in zones where the coastline is located parallel to wave direction. In this sense, the presence of Matxitxako cape, could act as a barrier, reducing the energy of the waves eastward of that point. The study carried out by Galparsoro et al., 2008 regarding the annual power distribution along the coast, could be correlated with waves energy along our study area ( Figure 10B). Their figure  (Figure 10B) confirms the low energy in the central segment of the analyzed coastal stretch. In the best preserved submerged marine terrace sequence (Figures 3C, 5B), it is also visible the effect of wave direction in terrace preservation. To the east of the headland (which would be protected from the waves) the levels are clearly preserved, whereas the western flank, exposed 2 | Submerged marine terrace distribution from Ea to Barrika in segments of the continental shelf with the most remarkable characteristics.

Continental shelf location
Location in Figure 5A Detected submerged terraces and shoreline angles depths to the main wave direction, exhibit submerged terraces that are not equally developed.

Continental shelf in front of Ea
From Barrika to the west no drowned shorelines were detected along most of the coast (Figure 10C). One fact that could explain the lack of shorelines could be that this coastal stretch is orthogonal to main wave direction and waves energy reaches its highest values at that sector. Such strong hydrodinamics conditions certainly would not favor the preservation of previously formed submerged terraces. Such observation would fit with observations made onshore as the coastal stretch Barrika to Punta Galea exhibits one of the best preserved and recognized rasa in our study area.
Nevertheless, the relation between high wave energy conditions and the lack of submerged marine terraces preservation, does not occur along all the studied shore. From Orio to the east and on the eastern flank of the Matxitxako cape, the modern coastline is perpendicular to wave direction (high hydodinamic), but some drowned marine terraces are well preserved. On the other hand, in front of Zumaia, Ondarroa and around Ogeia, with low hydrodinamics, the preserved submerged marine terraces distribution is different in each zone (Figure 10). Such observations suggests that other factors aside wave climatology influence the distribution of submerged marine terraces in this area.

Lithology
Several authors have pointed out on the effect of lithology and rock resistance in cliff recession rates, and hence in rocky shore platform generation (e.g. Sunamura, 1992;Woodroffe, 2002;Prémaillon et al., 2018). The main part of the studied coast is composed by flysch type rock successions from the Cretaceous period, and the Paleocene and the Eocene epochs. Secondly, at some sites, the Urgonian limestone (Cretaceous) also outcrops, with variable extent depending on the considered map (Figures 11A,B).
The stratification of flysch type rocks is notorious on the seafloor and can be followed in the hillshaded map, from Ondarroa to the east, with well developed submerged marine terraces ( Figure 11D); and between Barrika and Punta Galea where no submerged marine terraces are detected ( Figure 11C). Among flysh type sequences, some parts are more siliciclastic, more carbonaceous or the number and thickness of turbidites may vary, which affect to the generation and characteristics of rocky shore platforms (Baceta et al., 2012). Nevertheless, not such a big difference in submerged marine terrace preservation would be expected between those areas with quite similar lithological characteristics, suggesting that other factors are affecting the submerged marine terraces distribution offshore those coastal stretches. The differential preservation potential of the submerged terraces seems somehow controlled by lithology in the eastern part of Matxitxako cape ( Figure 11E). There, a submerged marine terrace is detected, in spite of being within a high wave energy area. We infer that there the seafloor could be carbonated, as onshore, there is a carbonated unit close to this area. The limestones, are comparatively more resistant to erosion than flysch, increasing the preservation potential of marine terraces. This carbonated unit could also act as a barrier, protecting the submerged terraces carved in presumably flysch unit substrate.

Stratification Direction
As stated before by Everard et al. (1964) or Trenhaile (1987) the dip, strike and thickness of the bedding influences the generation and development of modern rocky-shore platform. In our study area, the disposition of the stratification relative to the current shoreline has been highlighted in the littoral zone of the Basque Coast Geopark to explain the different morphologies of present coastline (Baceta et al., 2012).
On the continental shelf, submerged marine terraces are detected with variable bedding strikes. Parallel, oblique or perpendicular to the shoreline (e.g. Figures 11C,D, 12A-C). Note that around Punta Galea (Figure 11C), the stratification is perpendicular to the coastline and are not well preserved submerged marine terraces, whereas in other areas of the shelf under the same perpendicular strikes, submerged terraces are detected (e.g. Figure 11D).
The effect of the stratification is suggested for the area from Orio to the east, with high wave energy conditions and flysch type substrate, but where submerged terraces are preserved (Figures 12B,C). Along this coastal stretch, the bedding is oblique or sub-parallel along most of the continental shelf. According to Trenhaile (1987), the beddings which are parallel and with a high dip angle relative to the coast, are the most resistant to erosion. Such bedding disposition could explain why some submerged terraces are preserved in these areas of high hydrodynamics and not at other sites of similar wave energy conditions such as in front of Punta Galea, where the bedding is perpendicular to the wave direction ( Figure 11C).

Sea-Stacks Implications in Submerged Marine Terraces Sequence
The fact that the number of stacks is higher at certain depths, is probably associated to the submerged terrace levels, as some   Figure 7A. Dotted lines highlight stacks that could be associated with previously detected submerged terraces. White, T2; Orange, T5; Light blue, T7; Brown, T8; Dark blue, T10; Light green, T11.   The areas shallower than 30 m can give valuable information, as it is the area with the higher amount of sea-stacks and a shoreline was only detected at ∼ −13 m. Nevertheless, the number of sea-stacks at −19 and −26 m is considerably higher than at other depths. This can be because sea level has stayed for longer time at those depths, even though the submerged marine terrace morphology is not observed at present.
Aside from it, the presence of these stacks grouped at certain sites and depths, suggests the presence of former rocky shore platforms, which helps to complete the terrace sequence in places where nowadays no shorelines are present (Figures 7, 8 and Supplementary Figures S10-S15).

Submerged Marine Terraces That Fit With LEM Results
The submerged marine terraces generated on the scenario between MIS 5e and the Last Glacial Maximum, could fit with T3, T5, and T8. The main uncertainties come from two variables used to obtain those results. First, the best initial erosion rate for the generation of the staircase morphology was set at 1 m/yr, which is probably too high, taking into account the calculated cliff recession rates of 0.05 and 0.5 m/yr (mainly 0.15 m/yr) for the rocky shore platform in the rest of the Basque coast (Aubié et al., 2011). Second, the best staircase morphology is generated with a platform gradient set at 2 • . In profiles generated by the LEM with such initial slope, the width of the continental shelf between 0 and −120 m is ∼ 6 km ( Figure 9B), whereas the width of the continental shelf off the Basque coast most generally exceed 6 km (Figure 1).

Submerged Marine Terraces That Do Not Fit With LEM Results
T4, T7, and T10 were not detected in the scenarios tested in the LEM. In the case of T4 and T7, these are generally narrower submerged marine terraces, so we could infer shorter morphogenesis times, that are not reproduced in the models generated by the LEM. One hypothesis is, that they were generated between the Last Glacial Maximum (LGM) and present. We correlate the depth of those submerged terraces with the curves derived by Stanford et al. (2011) (Figure 13), that present episodes for which sea level seems to be rather stable and characterized by lower sea-level change rates ( Figure 13A). In other word, we compared the depth of these periods of sealevel stability (Figure 13B) with the depth of terraces identified from the bathymetry.
Sea-level change minima that could coincide with the described marine terraces are highlighted by the green and orange ellipses on Figure 13A, Those two periods would range between approximately 11.5 and 12.8 ka and between 9.8 and 10.4 ka. Then, using the sea-level curve of Stanford et al. (2011) (Figure 13B), the depths in those periods were calculated. Between 11.5 and 12.8 ka the sea level would be between −57 and −73 m in the 99% confidence curves, and between −60 and −68 m in the 67% confidence curves; whereas between 9.8 and 10.4 ka sea level would be between −35 m and −47 m in the 99% confidence curve. These depths could coincide with terrace levels T7 and T4, respectively.
The second hypothesis would be that these generally narrower submerged marine terraces were generated prior to the last transgression. In front of Ogeia, T-7 is barely distinguishable in the western part, while in the eastern part the sequence is best preserved (Figures 3C, 5B). Such morphological disymmetry for the submerged terraces could mean that the terraces were eroded during the last transgression, and hence, they were previously generated. Aside from this, considering that the sea stacks can be remnants of ancient rocky shore platforms, the ones grouped at depths of −61 m (Figures 8A,B), could indicate that previously a terrace was carved into that substrate. The occurrence of these stacks at that depth could support the idea that some terraces were generated prior to the last transgression.
Regarding T-10, due to its extension, planar character in some places, and all its associated stacks that are close to that submerged marine terrace depth, it could have been generated previous to MIS5e.
The terrace generation and erosion process could be as follows. T-10 could have been generated before MIS5e. Then, prior to the Last Glacial Maximum, T-3, T-4, T-5, T-7, and T-8 could be formed. With the last transgression, the ancient sea-cliffs and rocky shore platforms may have been smoothed or totally eroded.
Another factor that should be considered, is the inheritance or reoccupation. Some of the present rocky-shore platforms could have been generated not only during the present sea level, but also in past moments when sea level was similar to the recent one like at MIS5e Trenhaile, 2001Trenhaile, , 2002Trenhaile, , 2018Blanco Chao et al., 2003;Pedoja et al., 2018). A similar scenario may have occurred for submerged marine terraces which could be reoccupied at certain depths in several times during past sea level oscillations. In other word, the generation of the submerged terraces might be in part polygenic, and their generation could have started before MIS5e.
The reason for their preservation, could be related to a combination of wave climatology (energy and direction), lithology, and bedding strike at different segments of the continental shelf. Another point could be the general low slope of the shelf. A low slope platform would dissipate more wave energy, and waves wouldn't be able to erode completely the previous terraces, enhancing their partial preservation (Trenhaile, 2014;Jara-Muñoz et al., 2016).

CONCLUSION
The analysis of a high-resolution submerged DEM of the Bay of Biscay, reveals a staircase morphology, with irregularly preserved and distributed submerged marine terraces along the continental shelf. In total, 12 submerged terraces were detected. Six of them show more lateral continuity and are better preserved along the continental shelf; theirs shoreline angles are found between ∼ −34 and 41 m (T3), ∼ −44 and −47 m (T4), ∼ −49 and 53 m (T5), ∼ −59 and 62 m (T7), ∼ −65 and 67 m (T8) and ∼ −74 and −77 m (T10). We propose a novel analyze of drowned seastack analysis to detect marine terraces with poorly detectable shoreline angles.
Within the current stage of knowledge, we were not able to discriminate or to connect the submerged terrace and shoreline setting with a single variable. Rather, the degree of generation and preservation of fossil shorelines and submerged marine terraces, seems to depend on an interplay of variables such as wave climatology, lithology and bedding direction that would act simultaneously.
In absence of dating, we propose, through numerical modeling that submerged terraces T3, T5, and T8 could have been generated between MIS 5e and the Last Glacial Maximum. Nevertheless, the submerged terraces T4, T6, and T10, do not fit with that scenario. T4 and T6, could have been created during the last transgression, but the fact that at some places they seem to be partially eroded at those depths suggests a previous formation. Regarding T10, due to its extension in some places and the number of sea stacks that are associated to that submerged terrace, its generation probably took place before MIS5e.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
PB-L conducted the primary data analysis and prepared the manuscript. JJ-M generated the stacked sea level curve, the landscape evolution model, and helped in the data analysis, discussion, and in the preparation of the manuscript. KP, EI, IÁ, AA, and IG helped in the data analysis, discussion, and in the preparation of the manuscript.

FUNDING
PB-L is beneficiary of a predoctoral grant from the Basque Government (PRE_2018_2_0290). Financial support was also provided by the research project IT 1029-16-GBV6 from the Basque Government and US18/14 from the University of the Basque Country.