New insights into tunnel valley locations and Cenozoic exhumation in the southwestern Baltic Sea based on Machine Learning aided seismic refraction tomography

The Cretaceous and Cenozoic evolution of the North German Basin is shaped by complex processes involving basin inversion, uplift and erosion, extension and several periods of Quaternary glaciations. Based on a densely spaced long-offset 2D seismic profile network covering the Bays of Kiel and Mecklenburg, we employ a Machine Learning algorithm to pick refracted first-arrival travel-times. These travel-times are used in a travel-time tomography to derive velocity models for the approximately upper 800 m depth of the subsurface. Investigating velocity-depth relations within the Upper Cretaceous strata and analyzing lateral velocity anomalies within shallow depths provide new insights into the magnitude of the Cenozoic basin exhumation and the locations of glacial tunnel valleys. Our findings suggest that previously observed bent-up structures in seismic images are caused by heterogeneous velocities in the overburden and do not represent actual reflectors. We provide strong indications that these misinterpretations of imaging artifacts are related to tunnel valleys even though these valleys might not always be resolvable in seismic reflection or sediment sub-bottom images. Comparing Upper Cretaceous velocity-depth trends to reference trends reveals significantly higher velocities in our study area. We interpret these differences as overcompaction and estimate the apparent Cenozoic exhumation in the Bay of Mecklenburg to be about 475 m. Within the Bay of Kiel, we observe an increase of the apparent exhumation from about 385 m (south) to about 480 m (north). Our study demonstrates the importance of near surface velocity analysis for the investigation of geological processes in shallow marine settings.


Introduction
The North German Basin (NGB) is a sub-basin of the Southern Permian Basin and part of the Central European Basin System (Fig. 1). Its development since the Cretaceous consists of multifaceted processes, including basin inversion and several intermittent phases of extension and uplift accompanied by salt structure developments. Quaternary glaciations and tunnel valley evolution further overprinted the sedimentary record (e.g. Kley 2018;Ahlrichs et al. 2020Ahlrichs et al. , 2021Huuse and Lykke-Andersen 2000;Hübscher 2014, 2017). Understanding and quantifying these processes provides crucial insights for many considerations regarding the concepts of basin development as well as the economic role of the subsurface within the NGB. In this study, we focus on Cenozoic basin evolution in the Baltic sector of the NGB to improve the understanding of Paleogene and Neogene uplift events and the effects of Quaternary glaciations. This contributes to recent discussions regarding the nature, timing and magnitude of different uplift phases as well as the representation of glacially related structures in seismic images (e.g. Ahlrichs et al. 2021;Schnabel et al. 2021;Nielsen et al. 2002;Wenau and Alves 2020;Frahm et al. 2020). Therefore, we investigate the locations of glacial tunnel valleys and estimate the magnitude of the basin exhumation during the Cenozoic.
Glacial tunnel valleys are channel shaped features that incise the underlying strata and result from glacial meltwater transport (e.g. Huuse and Lykke-Andersen 2000; Kirkham et al. 2021). Investigating their locations, orientations and dimensions offers a basis for further work regarding their genesis and the extension of glacial ice advances and retreats (Kehew et al. 2012;van der Vegt et al. 2012). Additionally, the investigation of tunnel valley locations is important with respect to identifying imaging artifacts in seismic images. Al Hübscher (2014, 2017) observed and mapped several of these valleys in the Bay of Kiel by means of seismic reflection data (left inset in Fig. 1). The authors noted anticline-like upward bent strata and near-vertical faults underneath the identified valleys. Similar bent-up structures and faults were also observed without recognizable valleys above them and a geometric correlation between valleys, valley infill deposits and bent-up structures could not be identified. The authors also assumed soft, unconsolidated and gassy sediments with low velocities as infill deposits and therefore did not postulate a causal relationship between valleys and the bent-up structures and faults beneath (for details see Al Hseinat and Hübscher 2014). These observations were recently revised by Frahm et al. (2020), who showed that some bent-up structures are imaging artifacts resulting from high velocity infill within the tunnel valleys above. They suggested that this might also hold for bent-up structures not accompanied by visible valleys, as some tunnel valleys may be too small to be resolvable in seismic images. The investigation of bent-up structures without observable valleys above, their potential for misinterpretation and their relation to glacial tunnel valleys represents a major objective of our work.
The investigation of the magnitude of the basin exhumation during the Cenozoic constitutes a second main objective. It represents a challenge, as uplift related erosion and non-deposition result in the lack of synkinematic units within the sedimentary record. However, the influence of exhumation is retained in the compaction status of the preserved sediments. The porosity of sedimentary rocks decreases during the deeper burial of the strata due to increasing pressure and resulting compaction (Bulat and Stoker 1987). Consequently, this implies an increase in the seismic velocity. This process is thought to be irreversible and thus enables the use of velocities as proxy for the maximum burial depth (Corcoran and Doré 2005). Post-depositional uplift and erosion decreases the burial depth and results in apparent overcompaction relative to the present depth of burial and thus, higher than expected velocities at this depth. The resulting velocity anomaly can be used as a measure of the post-depositional reduction of the thickness of the overburden (exhumation) relative to the unit in question (Japsen 1998(Japsen , 2000. We aim to derive exhumation estimates for the Cenozoic for the Bay of Kiel and Bay of Mecklenburg and to investigate the impact of the northern basin margin of the NGB onto this exhumation. Other authors have already investigated the Cenozoic uplift within Fig. 1 Structural overview of the Baltic sector of the North German Basin and its adjacent sedimentary basins, modified after Ahlrichs et al. (2021). The red box marks the study area and the inset on the upper left displays the locations of the used marine seismic profiles as well as previously observed bent up structures. The orange marked profile was used for training the Machine Learning algorithm. The inset on the upper right shows the approximate location of the study area within the Southern Permian Basin. Purple shading indicates the Central European Basin System. Bent up structures after Al Hübscher (2014, 2017) the Danish areas north of our study area and within the Bay of Mecklenburg (Japsen 1998;Rasmussen 2009;Schnabel et al. 2021). This enables a detailed comparison and allows conclusions about the regional variations of the exhumation magnitude within and adjacent to our study area.
We base our investigation of tunnel valley locations and Cenozoic exhumation on a dense grid of high resolution and long-offset seismic data. Seismic imaging offers a powerful tool to investigate the subsurface and derive the main effects that shaped its evolution. One of the fundamental prerequisites for seismic data analyses are reliable seismic velocities. In this study, we present a data-driven scheme based on seismic refracted first-arrival travel-times to calculate velocity models for several marine seismic profiles located in the Baltic sector of the NGB (Fig. 1). We employ a machine learning (ML) based Convolutional Network to pick refracted first-arrivals and subsequently use the picks in a travel-time tomography to derive velocity models. We then analyze the velocity distributions and variations above the observed bent-up structures and within the Upper Cretaceous unit in order to investigate the impact of Quaternary glaciations and Paleogene and Neogene exhumation events. Our study area in the Baltic sector of the NGB offers unique prerequisites for our study. With the available long-offset seismic data we can derive high-resolution velocity models in a purely data-driven fashion. Previous works from other authors in this area allow us to validate our applied methods and compare the refraction travel-time tomography approach with other methods, such as migration velocity analysis (Schnabel et al. 2021), offset-depending seismic imaging (Frahm et al. 2020) or analyzing well information (Japsen 1998, van Dalfsen et al. 2006).

Geological setting
Our study area is located in the Baltic sector in the northeast of the NGB and comprises the Bay of Kiel in the west and the Bay of Mecklenburg in the east (Fig. 1). It is part of the intracontinental Southern Permian Basin spreading across large parts of Central Europe. The western part of the Bay of Kiel comprises the Eastholstein Trough, which is part of the NNE-SSW trending Glückstadt Graben that marks the western boundary of the study area and the transition from the Baltic to the North Sea sector of the NGB (Maystrenko et al. 2005a, b). The northern basin margin is composed out of the Ringkøbing-Fyn High and the Møn High representing basement highs stretching in approximately E-W direction (Berthelsen 1992;Clausen and Pedersen 1999). The Grimmen High located east of the Bay of Mecklenburg marks the eastern margin of the study area and of the NGB (Fig. 1).
Basin formation of the NGB started in the Late Carboniferous to Early Permian with extensive volcanism, wrench faulting, lithospheric thinning and subsequent basinwide thermal subsidence lasting until Middle Triassic times (Maystrenko et al. 2008;Ziegler 1990;Kossow and Krawczyk 2002). During the ongoing subsidence, the basin was filled with terrestrial Rotliegend sediments, Zechstein evaporites and red-bed clastics (Fig. 2). During the Middle Triassic, a regional sea-level rise caused the development of Muschelkalk carbonate platforms until a decrease in sea level during the Late Triassic reestablished terrestrial conditions and resulted in the sedimentation of the Keuper sequence (Kossow and Krawczyk 2002). In the Late Triassic, large-scale E-W directed extension triggered intensive salt movement in the study area (Hansen et al. 2005;Hübscher et al. 2010;Maystrenko et al. 2005b). Following a major transgression in the Rhaetian, Early Jurassic sediments were deposited under shallow marine conditions (Fig. 2) (Maystrenko et al. 2008;Stollhofen et al. 2008).

Phases of inversion and uplift from Jurassic to Neogene
The long-term subsidence regime changed in the Middle Jurassic. Induced by the central North Sea Doming event, the study area was uplifted resulting in a period of nondeposition and erosion (Underhill and Partington 1993;Hansen et al. 2005Hansen et al. , 2007Hübscher et al. 2010). Almost the entire Jurassic sequence as well as parts of the Upper Triassic successions were eroded (Kossow and Krawczyk 2002;Hübscher et al. 2010;Ahlrichs et al. 2020). An eustatic sea level rise during the Early Cretaceous initiated resumed sedimentation in the Albian and marks the onset of a period of relative tectonic quiescence lasting until the Turonian (Fig. 2) (Kossow and Krawczyk 2002;Hübscher et al. 2010;Ahlrichs et al. 2021). Starting in the Turonian to Santonian, a major plate reorganization marked by the onset of the Africa-Iberia-Europe convergence led to NNE-SSW directed intraplate thrusting transmitting compressional stress into the European foreland basins leading to basin inversion, uplift and erosion (Kley and Voigt 2008;De Jager 2007;Kley 2018). The Grimmen High was uplifted by about 450 to 570 m and minor salt remobilization occurred in the Bay of Mecklenburg (Kossow and Krawczyk 2002;Ahlrichs et al. 2021). Due to the approximately parallel alignment of the Glückstadt Graben compared to the principle stress regime, the transmitted stress caused only minor uplift and salt movement in the Glückstadt Graben area (Maystrenko et al. 2005a, b;Ahlrichs et al. 2021).
A renewed period of large-scale uplift during the Paleocene, most likely caused by domal uplift in central Germany, eroded Selandian and parts of the Danian sediments in the study area von Eynatten et al. 2021). The subsequent period until the middle Eocene is marked by relative tectonic quiescence, comparable to the structural style following the Jurassic/Early Cretaceous uplift event ). This quiescence persisted until the late Eocene, when almost E-W directed extension, possibly related to the opening of the European Cenozoic Rift System, reactivated salt movement in the study area lasting until the Miocene (Maystrenko et al. 2005b;Ahlrichs et al. 2021). Contemporaneously to the renewed salt activity in the late Eocene, a period of increased tectonic activity driven by the Africa-Iberia-Europe convergence resulted in N-S directed shortening leading to inversion and uplift of some basins within the Southern Permian Basin system. Thereby, the exact impact onto the NGB remains unclear and aspect of recent discussions (Kley and Voigt 2008;Kley 2018). During the Miocene, the principle stress regime in the study area changed towards the currently persisting NW-SE extension Al Hseinat and Hübscher 2017). Consequent reactivation of preexisting fault zones due to stress transmitted from the Alpine orogeny resulted in reshaping the margins of the North Sea Basin and uplift of the hinterland, which resulted in an uplift of the Ringkøbing-Fyn High of about 300 m (Japsen 1998;Rasmussen 2009). This caused erosion of much of the Miocene, Oligocene and parts of the Eocene deposits in the study area. The overall uplift of the Ringkøbing-Fyn High and the Bay of Mecklenburg resulting from all observed uplift events during the Cenozoic is estimated at about 600 m and 450 m, respectively (Rasmussen 2009;Schnabel et al. 2021).

Quaternary development
Since the mid-Quaternary, the development of the NGB was affected by at least three major glaciations with ice sheets fully covering the study area (e.g. Huuse and Lykke-Andersen 2000; Ehlers et al. 2011;Houmark-Nielsen 2011). Accompanied by the ice advances, glaciers eroded large parts of the pre-glacial Quaternary and Neogene sedimentary units (e.g. Sirocko et al. 2008). Glacial landforms that can be observed throughout the southwestern Baltic Sea include moraines, ice-marginal valleys, tunnel valleys and the effects of fluvial erosion (Sirocko et al. 2008;Al Hseinat and Hübscher 2017). Tunnel valleys primarily occur as cutand-fill structures filled with Quaternary deposits such as till, glaciolacustrine clay and silt and/or meltwater gravel and sand. Several authors investigated the velocities of valley infill deposits and observed infills with higher velocities than the surrounding sediments (e.g. Jørgensen et al. 2003;  Kossow and Krawczyk (2002) and Kley (2018). Sea-level curves based on Haq et al. (1987) and Haq (2014Haq ( , 2017Haq ( , 2018 Gerok et al. 2017;Kristensen and Huuse 2012). Especially glacial till is thought to be responsible for this observation (Jørgensen et al. 2003). Kirkham et al. (2021) investigated the internal structures of tunnel valleys and discovered complex internal patterns that are suspected to result from several phases of deposition, erosion and overprinting during possibly more than one glaciation. Therefore, structures and infill deposits can vary considerably between but also within tunnel valleys.

Seismic database and stratigraphy
The seismic data we use in this study originates from the "BalTec" dataset collected during the MSM52 cruise led by the University of Hamburg in cooperation with the Federal Institute for Geosciences and Natural Resources (BGR), University of Greifswald, Polish Academy of Sciences, Uppsala University and German Research Centre for Geosciences Potsdam (Hübscher et al. 2016). The dataset consist of a network of high-resolution 2D multichannel seismic profiles located in the Baltic Sea between the Bay of Kiel and southeast Gotland. The acquisition layout included an active seismic streamer length of 2700 m and thus enabled the recording of refracted waves. The receiver spacing was 12.5 m and the source station interval amounted to approximately 25 m. From this dataset, we selected 19 seismic profiles covering the Bay of Kiel and Bay of Mecklenburg for the velocity model calculation and analysis in this study. Seismic images corresponding to these profiles were created following the processing scheme described by Ahlrichs et al. (2020).
The stratigraphic subdivision of the individual units in the seismic images and velocity models follows the extended stratigraphic framework of Ahlrichs et al. (2020Ahlrichs et al. ( , 2021. The stratigraphy is thereby based on and validated with several deep research and hydrocarbon exploration wells as well as results from previous studies in the area (for details see Ahlrichs et al. (2020Ahlrichs et al. ( , 2021 and references therein). In this study, we focus on six seismic horizons ranging from the Upper Cretaceous to the Quaternary, namely the base Quaternary, base Eocene, base Paleocene, base Maastrichtian, base Campanian and top Turonian (Fig. 2).

Theory of ML based first arrival picking and travel-time tomography
Calculating velocity models for the selected 19 seismic profiles requires to determine the exact first-arrival times of the refracted P-waves. In case of an increasing velocity with depth, refracted waves propagate through the subsurface at layer interfaces with the velocity of the lower layer (Kearey et al. 2002). Due to their relation to the layer velocities, they offer first hand information about the velocity structure in the subsurface. These information can be extracted by picking the first-arrival travel-times and then modeling ray paths and inverting for velocity structures that correspond to the picked travel-times (Hübscher and Gohl 2014).
At offsets larger than the crossing point with the direct wave (see e.g. Fig. 3.13b in Kearey et al. (2002)), refracted P-waves are always the first arrivals in seismic shot gathers. Picking the first-arrivals manually is time-consuming; therefore, a ML approach represents a powerful and efficient way to pick these travel-times. ML algorithms are able to automatically detect patterns in data and then use these patterns to classify the data in a desired way (Murphy 2012). In this study, we employ a ML algorithm that is inspired from and follows the general architecture of the U-Net from Ronneberger et al. (2015). The U-Net is a Convolutional Network in form of an autoencoder, which breaks down highlevel features into simpler representations and then restores the original format while retaining the new representations. This architecture helps the algorithm with the classification of complex input data (Goodfellow et al. 2016). The main aspects of the contracting path of the autoencoder consist of five sequences of two consecutive convolutions (3 × 3 kernel), each followed by a rectified linear unit (ReLu) activation function, and a max pooling operation (2 × 2 kernel). The last contracting sequence does not contain max pooling, instead the expansive path of the autoencoder starts with an up-convolution (2 × 2 kernel) followed by two consecutive convolutions (3 × 3 kernel), which are again each followed by a ReLu function. This expansive sequence is repeated four times. Subsequently, an additional convolution (1 × 1 kernel) and a softmax activation function results in the desired probability distribution, which acts as classification (for details see Fig. 3.6 in Ronneberger et al. 2015). In contrast to the setup from Ronneberger et al. (2015), we use the Adam optimizer as optimization procedure during the training of the algorithm. It is more time efficient, easy to implement and requires little memory space (Géron 2019;Kingma and Ba 2014). In addition, we incorporated batch normalization, which can further improve the optimization procedure as well as the regularization performance (Goodfellow et al. 2016).
We use the first-arrival refracted P-wave travel-times to calculate velocity models by performing travel-time tomography. Therefore, we apply the PSTomo_eq algorithm developed and described by Tryggvason (1998) for local earthquake tomography. It was later enhanced for the specific use of controlled source data with known source parameters (Bergman et al. 2004;Tryggvason et al. 2009). The forward problem is based on the wavefront tracking principle introduced by Vidale (1988) and is solved by finding a first order finite difference solution to the Eikonal equation, which results in a first-arrival travel-time field (Podvin and Lecomte 1991). Ray paths are found by tracing the rays backwards from the receivers to the source location following the paths perpendicular to the isochrones (Hole 1992). This approach is able to calculate unconditionally stable and accurate first-arrival travel-times even in structures with large velocity contrasts and arbitrary shapes (Podvin and Lecomte 1991;Buske and Kästner 2004). The inverse problem employs the iterative LSQR conjugated gradient solver by Paige and Saunders (1982) to minimize the objective function (Eq. 13 in Tryggvason 1998). The model smoothness can be controlled by the number of LSQR steps in the inversion and by the weighting factor for the smoothness constraints applied to the objective function (Eqs. 9 and 13 in Tryggvason 1998). We slightly modified the PSTomo_eq algorithm by adding an additional code fragment that allows to change the inversion grid cell size during the iterative calculation of the velocity models. This helps to avoid stability issues and allows to use a small inversion grid size (Vesnaver and Böhm 2000).

Velocity model calculation
Obtaining refracted first arrival travel-times for the entire seismic dataset by using the described ML algorithm requires to train the algorithm with manually picked firstarrival travel-times. We perform manual picking on 2301 seismic shot gathers from one seismic profile located in the Bay of Kiel (Fig. 1). The manual picking accuracy is estimated between 2 and 4 ms depending on the data quality.
Manually picked travel-times then serve as labels during the supervised training of the ML algorithm and are divided into three subsets. Labels from 400 shot gathers are assigned to a test subset, which is used to assess the performance of the fully trained algorithm. Out of the remaining labels, 80% are used to train the algorithm and 20% are used for validation during the training. Training is an iterative process, which is stopped after 40 training epochs as the training and the validation error starts to diverge. This indicates that more trainings epochs will not lead to a further improvement of the algorithm (e.g. Goodfellow et al. 2016). Subsequently, we apply the fully trained algorithm to the test subset to quantify the picking performance. ML based picks that are within the maximum estimated manual picking accuracy of 4 ms are regarded as accurate first arrival travel-time picks. This is the case for over 91% of the picked travel-times. Figure 3 illustrates an exemplary shot gather with manually picked labels and picks from the ML algorithm. Despite some missing picks and minor differences at large offsets (indicated by arrows and circle in Fig. 3), the overall performance of the ML based picking is satisfying.
For all seismic profiles and using all shot gathers, we use the picked refracted first-arrival travel-times to calculate velocity models by applying travel-time tomography. For simplicity, the initial velocity model is identical for all seismic profiles and varies only with depth. The initial velocity increases from 1500 m/s at the surface to 3500 m/s in 1000 m depth and is based on previous knowledge of the velocity structure in the study area (e.g. Ahlrichs et al. 2021;Schnabel et al. 2021). The tomography uses 16 iterations, whereby the last iteration consist of forward modelling only, therefore yielding the final root mean square (RMS) error of the velocity models. The forward modelling uses a uniform grid cell size that is set to 20 × 20 m in the first six iterations and 10 × 10 m in the last ten iterations. The inversion grid cell size changes after every third iteration; from 100 × 40 m in the first to 20 × 10 m in the last inversion iteration. Model regularization is controlled by a set of parameters that aim to find a compromise between minimizing the RMS error and calculating smooth velocity models. With these settings, we obtain velocity models with final RMS errors varying between 3 ms and 5 ms. An exemplary velocity model from the Bay of Mecklenburg is illustrated in Fig. 4a. All velocity models are calculated in the depth domain, Fig. 3 Comparison between manually picked travel-times from the test subset (green) and ML based picks (blue) for an exemplary shot gather. Circle indicates traces with slightly differing travel-time picks between offsets of around 2400 and 2500 m. Arrows denote traces with no ML based picks. For this shot gather, a total of 23 traces were affected by differing or no ML based picks, whereby the maximum mismatch does not exceed 10 ms however, as some seismic images and stratigraphic horizons are only available in the time domain, the velocity models are additionally converted into the time domain to allow a joint interpretation.
In order to assess the resolution and reliability of the obtained velocity models, we perform checkerboard tests for selected profiles. Thereby, we use two different checkerboard grids to assess the velocity models for different scales. Figure 4b,c shows the starting checkerboard model and the recovered checkerboard model for the same profile as illustrated in Fig. 4a with a checkerboard grid of 2100 × 80 m. Within the first 200 m depth, the checkerboard grid is well recovered with mostly sharp and distinguishable edges and only minor smearing effects. With increasing depth, the quality of the recovered grid decreases. Nonetheless, the general pattern remains visible in most parts of the recovered checkerboard model. Similar results are obtained using a checkerboard grid of 4000 × 120 m. We therefore assume the velocity models within the first approximately 200 m to be reliable and well resolved. However, with increasing depth the reliability decreases and the velocity models may only be able to represent the magnitude of the velocities and large-scale velocity variations.

Observations and analyses
Building upon the work of Frahm et al. (2020), we analyze the derived velocity models to investigate the Cenozoic and especially the Quaternary velocity structure of the Baltic sector of the NGB. We strive to identify velocity variations in the models that can be linked to bent-up structures and glacial tunnel valleys by superpositioning the final models with the corresponding stratigraphically interpreted seismic images. Additionally, we analyze velocity variations with depth within the Upper Cretaceous units to estimate burial depth anomalies for the Bay of Kiel and Bay of Mecklenburg. This aims to obtain estimated exhumation values for both parts of the study area during the Cenozoic.

Velocity analysis of Quaternary deposits
Almost all velocity models reveal laterally constraint positive velocity anomalies within the upper parts of the velocity models. Figures 5 and 6a,b show several examples of these velocity anomalies within the Bay of Kiel. All of these anomalies spatially correlate to observable bent-up structures in the seismic images (time domain) directly below the anomalies (Figs. 5 and 6a,b). Additionally, some of the  a). c Recovered checkerboard model with areas with no ray coverage displayed in white velocity anomalies also correlate to valley structures visible in the seismic images close to the seafloor (blue box in Figs. 5 and 6a). However, this is not the case for all anomalies. In some cases, there are no valleys recognizable in the seismic (red box in Figs. 5 and 6b) and in other cases, the velocity anomalies correlate to only parts of visible valleys as it is the case for the anomaly highlighted by the blue box in Fig. 5. In this case, the anomaly is constraint to the western part of the valley, while in the eastern part no distinctive anomaly can be recognized. None of the observed velocity anomalies is constraint at the base of a visible valley or the base Quaternary horizon but they gradually integrate into the velocities observed within the Eocene and Paleocene successions (Figs. 5 and 6a).
The amplitudes of the observed positive velocity anomalies exhibit considerate variations. Figure 6a,b shows two anomalies from the Bay of Kiel and corresponding lateral 1D velocity profiles through the velocity models at a depth of 150 ms TWT. While the anomaly on the left (Fig. 6a) reveals a local velocity increase of more than 200 m/s, the anomaly in the center corresponds only to a velocity increase of about 100 m/s (Fig. 6b). The dimension of the bent-up structures in the seismic images correlate with these differences with the bent-up structure on the left being much larger than the one in the center. Figure 6c on the right shows an example of a bent-up structure with a larger lateral extension. It is accompanied by a very faint positive velocity anomaly that is only recognizable in the lateral 1D velocity profile. A tunnel valley is only visible at the location of the largest velocity anomaly (Fig. 6a). Within the 19 investigated velocity models, we observed a total of 35 positive velocity anomalies with varying shapes and amplitudes (Fig. 8).

Burial depth analysis
To investigate exhumation in the Bay of Kiel and Bay of Mecklenburg during the Cenozoic, we analyze velocity variations with depth for the Upper Cretaceous unit. Due to the basin configuration, the strata can be encountered at various depths, which allows a depth dependent analysis of the velocities. For the analysis, we divide the Bay of Kiel into a southern and a northern part to include possible north-south variations across the northern basin margin of the NGB. For each location along the profiles (CMP positions), we utilize the mean layer depth (Z) and the mean layer velocity (V Z ). Similar to the approach of Marsden (1992) and Schnabel et al. (2021), we determine linear velocity-depth relations by using Eq. 1: Applying a linear regression yields the velocity-depth gradient k and the surface velocity V 0 . The results for both parts of the study area are summarized in Table 1.
In order to obtain velocity anomalies from the derived velocity-depth relations, they are compared to the velocitydepth baseline trend for the Upper Cretaceous Chalk Group of Japsen (2000). It is assumed that this trend corresponds to strata that has not been subject to any uplift or erosion during the Cenozoic. Figure 7 illustrates the Z−V Z tuples for each CMP, the derived velocity-depth relations and the reference baseline trend for the northern and southern part of the Bay of Kiel (left) and the Bay of Mecklenburg (right). All velocity-depth relations from this study considerably differ from the reference baseline trend. Measured at the mean depths covered by the Upper Cretaceous unit, the velocity anomalies within the Bay of Kiel represent an increased velocity of approximately 320 m/s in the south and 440 m/s in the north when compared to the reference baseline trend. In the Bay of Mecklenburg, the observed velocity is 360 m/s higher than predicted by the reference trend. Using the derived velocity anomalies as well as the calculated velocity-depth gradients k (Table 1), we obtain (1) V Z = V 0 + kZ burial depth anomalies of approximately − 385 ± 12 m in the southern Bay of Kiel, − 480 ± 6 m in the northern Bay of Kiel and − 475 ± 5 m in the Bay of Mecklenburg. Thereby, the error margin of these anomalies only represents the uncertainty imposed through the linear regression and is not able to include errors related to possible uncertainties within the velocity models or to assumptions made for the method itself.

Tunnel valley locations
In the Bay of Kiel, the locations of observed positive velocity anomalies correlate closely to previously published bent-up structures by Al Hseinat and Hübscher (2017;their Fig. 23) (Fig. 8). Frahm et al. (2020) linked these bent-up structures underneath visible tunnel valleys (green structures in Fig. 8) to high velocity infill within these valleys and showed that the bent-up structures represent velocity pull-ups and thus are seismic imaging artifacts. The authors argued that the comparably high valley infill velocities result in shorter travel-times for rays travelling through the valleys and consequently cause the upwards bent structures underneath the valleys. Given the direct spatial correlation of bent-up structures underneath visible valleys and the observed positive velocity anomalies (Fig. 8), this study confirms the conclusion of Frahm et al. (2020) by the application of the automated refraction-tomography. Additionally,  Regarding the bent up structures without visible tunnel valleys above (purple structures in Fig. 8), Frahm et al. (2020) proposed that they could be velocity pull-ups caused by small and seismically unresolved tunnel valleys. This assumption can be strongly supported with our tomography based velocity analysis. In about two thirds of the locations of bent-up structures without visible valleys above, we observed positive velocity anomalies, which indicate shallow high velocity deposits (Fig. 8). For the other bent-up structures, there are no or only very faint anomalies recognizable (see Fig. 6c and circle in Fig. 8). As the dimensions of these bent-ups are considerably smaller than for the bent-up structures with distinct velocity anomalies above (Fig. 6a,c), it is likely that the thickness of the shallow high-velocity deposits is insufficient for the tomography to resolve velocity anomalies. Nonetheless, the velocities derived through our automated refraction tomography enable a much broader and time-efficient analysis of the bent-up structures and their classification as velocity pull-ups. We regard the combination of bent-up structures and positive velocity anomalies as an adequate indicator for glacial tunnel valleys in the southwestern Baltic Sea. Velocities derived through our automated refraction tomography are especially valuable in regions with shallow water, as the method is not limited by physical constraints such as NMO-stretch or over-critical reflections, which commonly hampers the processing and velocity analysis of multi-channel seismic reflection data. Given the observed velocity anomalies above bent-up structures in the Bay of Mecklenburg, we assume that several glacial tunnel valleys exist in that area as well.

Cenozoic exhumation
The observed velocities for the Upper Cretaceous strata are high relative to their present burial depth and thus indicate that these strata are overcompacted. This suggests that the derived burial depth anomalies are measures of the magnitude of the Cenozoic exhumation, implying an exhumation of about 475 m in the Bay of Mecklenburg and between 385 and 480 m in the Bay of Kiel. Even though the analysis of Yellow marker denote locations where we observed positive velocity anomalies in the velocity models. The circle indicates a part of a bent up structure with only one observed velocity anomaly despite several crossing seismic profiles burial depths is well established (e.g . John 1975;Scherbaum 1982;Japsen 1998;Schnabel et al. 2021), the interpretation and conclusion is not straightforward. Assumptions made during the calculation of the burial depth anomalies include an assumed linear velocity increase with depth, thus not considering the influence of pore fluids, geothermal conditions or chemical processes on the velocities. Additionally, uncertainties of the velocity models in greater depths as indicated by the checkerboard tests may influence the derived burial depth anomalies as well. The error margins of the linear regression cannot represent these uncertainties. The actual uncertainties related to the obtained exhumation values may be larger and therefore we interpret the obtained burial depth anomalies as estimates of the general magnitude of the exhumation.
Other studies who investigated exhumation during the Cenozoic at the northern margin of the NGB and adjacent areas have come to comparable estimations (Fig. 9). This highlights the value of our refraction tomography based velocities and the applicability of the burial depth analysis. Schnabel et al. (2021) derived burial depth anomalies from migration velocity analysis and estimated 450 m of Cenozoic exhumation for the Bay of Mecklenburg. Rasmussen (2009) analyzed 2D and 3D seismic data and well information to estimate the uplift of the Ringkøbing-Fyn High during the Cenozoic at about 600 m. Similar results were derived by Japsen (1998), who used well information to create a map of burial depth anomalies for the North Sea Basin and its adjacent areas based on velocity-depth data for chalk. The author interpreted this as measure of uplift and erosion and assumed a magnitude of exhumation between 500 and 750 m for the northern part of the Bay of Kiel and the area just north of the Bay of Mecklenburg. The general outline of the burial depth map indicates that the exhumation in the Bay of Kiel tends towards 500 m rather than 750 m (see Fig. 9 in Japsen 1998). The exhumation during the Cenozoic at the Rødby-1 and Ørslev-1 wells located north of the Bay of Mecklenburg was determined with 634 m and 483 m, respectively (Japsen 1998).
In general, all presented studies indicate basin exhumation in the Baltic sector of the NGB and its adjacent areas during the Cenozoic. While the magnitude of exhumation from Schnabel et al. (2021) for the Bay of Mecklenburg and from the Ørslev-1 well correlate very well with the results from this study, the results from Rasmussen (2009), Japsen (1998 and the Rødby-1 well are slightly larger. As the results from Rasmussen (2009) andJapsen (1998) were not specifically derived for our study area, a variation due to the spatial distance is possible. The analysis within the Bay of Kiel reveals an increase of the derived apparent exhumation from about 385 m in the south to about 480 m in the north. Hence, this trend shows increased exhumation towards the basin margin, whereas the southward directed basin center is less affected (Scheck and Bayer 1999). Additionally, this trend is in agreement with higher exhumation values north of our study area from Rasmussen (2009) andJapsen (1998).
The applied method does not allow to quantitatively assign or split the obtained exhumation values to specific Fig. 9 Summary of the obtained apparent exhumation values from this study and from Japsen (1998), Rasmussen (2009), Schnabel et al. (2021 as well as from two Danish onshore wells (see Japsen 1998) events since the late Cretaceous. From the literature, those exhumation events are known from other areas in Central Europe, which could be considered to interpret the here estimated exhumation values within the regional tectonic framework. The late Eocene -Oligocene marks a phase of basin inversion, uplift and erosion in Central Europe (e.g. Kley 2018). However, contemporaneous E-W extension is observed in our study area  and it is therefore unlikely that the estimated Cenozoic large-scale exhumation occurred during this time. Other exhumation events are reported for the Paleocene and within the Neogene (von Eynatten et al. 2021;Rasmussen 2009;Japsen et al. 2007). Paleocene exhumation is characterized by largescale domal uplift centered around central Germany and related to dynamic topography (von Eynatten et al. 2021). The causes and the timing for the exhumation events during the Neogene are not clearly determined. Japsen et al. (2002) related the exhumation of southern Scandinavia to two phases, an early phase during the Miocene and a later and better constrained period of exhumation during the Pliocene. The authors attribute the exhumation to the rise of the southern Swedish Dome. Jones et al. (2002) suggested that the relative movement of the Iceland Plume towards Europe might have influenced Miocene uplift in northern Europe. Other authors related early Neogene exhumation to the reactivation of former structural elements such as the Sorgenfrei-Tornquist Zone (Fig. 1) (Rasmussen 2009;Eidvin et al. 2014). As a driving mechanism for this reactivation, Rasmussen (2009) discussed a relation to transmitted stress within the European foreland originating from the Alpine orogeny but also noted that an influence from the opening of the North Atlantic cannot be ruled out. The early Neogene exhumation coincides with the reactivation of salt structures in the Danish area to the northwest of our study area and in the Glückstadt Graben to its west (Rasmussen and Dybkjaer 2020;Hinsch 1986). While Ahlrichs et al. (2021) assumed that significant salt remobilization in our study area took place during the late Eocene, they did not rule out an additional phase of salt movement during the Neogene. Following exhumation and reactivation of salt structures in the early Neogene (early Miocene), there are indications for subsidence within several graben structures in the vicinity of our study area during the middle Miocene (e.g. graben of the Poznań-Oleśnica Fault Zone, Glückstadt Graben) (Widera et al. 2008;Maystrenko et al. 2005a). This subsidence may be related to local extension connected to the European Cenozoic Rift System. For the late Neogene (early Pliocene), Japsen et al. (2007) observed exhumation in the Norwegian-Danish Basin by analyzing fission track data. They suggest that exhumation of the hinterland continued during regional tilting in the late Pliocene. This is in agreement with observations from Graversen (2002), who noted exhumation of the southern Swedish Dome area in the late Neogene and into the Quaternary. Recent studies specified the onset of Neogene uplift within and north of the Sorgenfrei-Tornquist Zone to the early Miocene, while Neogene uplift south of the Sorgenfrei-Tornquist Zone began only in the Pliocene (Green et al. 2022). Thereby, the Sorgenfrei-Tornquist Zone is thought to have acted as hinge zone between the uplifted Fennoscadian Shield and subsiding Danish Basin (Green et al. 2022;Japsen et al. 2018). Given the complex tectonic evolution, we cannot constrain our observed exhumation to a specific event. However, our observed exhumation is in agreement with km-scale exhumation reported for the Paleocene to the south and Neogene to the north of our study area and represents the superposition of the impact of these events on the burial history of the region. For the Danish Basin and Ringkøbing Fyn High, Rasmussen (2009) proposed that Neogene exhumation exceeded the magnitude of the Paleocene exhumation by a factor of three. Following this interpretation, we assume that a significant amount of the observed exhumation in our study area took place during the Neogene, even though the impact of Paleocene uplift likely increases towards the south. End members for the underlying processes controlling Cenozoic exhumation in Central Europe might be either dynamic topography induced by mantle circulation or far-field stress induced by plate motions (Green et al. 2022). To clarify the underlying process of vertical motions in Central Europe during the Cenozoic and the possible link of local to regional observations of exhumation and burial, further work on the interaction of plate tectonics and mantle circulation is needed.

Conclusions
We presented a data driven scheme to pick refracted firstarrival travel-times by applying a ML algorithm and use these travel-times to calculate velocity models for 19 seismic profiles covering the Bay of Kiel and Bay of Mecklenburg. Based on the derived seismic velocity models, we investigated local velocity anomalies located above previously observed bent-up structures. In addition, we also used the derived velocities as well as a reference baseline trend for Upper Cretaceous chalk to calculate Cenozoic burial depth anomalies for both parts of the study area. Our main conclusions are: • Previously observed bent-up structures in seismic images are imaging artifacts, which are caused by high velocity infill within glacial tunnel valleys creating velocity pull-ups below them. The absence of visible valleys and positive velocity anomalies above some pull-ups is thought to be a result from small, shallow valleys, which could not be resolved due to the limited vertical resolution of the seismic data and the refraction travel-time tomography, respectively. • The combination of bent-up structures and positive velocity anomalies can serve as an indicator for glacial tunnel valleys in the south-western Baltic Sea. We therefore conclude that several tunnel valleys exist in the Bay of Mecklenburg as well. • Machine Learning aided seismic refraction tomography is an appropriate tool to overcome physical limitations of multi-channel seismic data in shallow waters, such as over-critical reflections and NMO-stretch, and enables a data-driven calculation of velocity models. • The study area experienced large-scale exhumation during the Cenozoic. The amount of apparent exhumation is estimated to about 475 m in the Bay of Mecklenburg and between 385 and 480 m in the Bay of Kiel. Within the Bay of Kiel, the exhumation increases from the south towards the basin margin in the north. • The obtained apparent exhumation estimates are in agreement with other studies on burial depth anomalies, which were based on well data or migration velocity analysis. This study shows the applicability of refraction tomography derived velocity analysis to infer overcompaction trends. Exhumation estimation in this study is solely based on seismic images and stratigraphic horizons in the time domain. For depth ranges that can be imaged with the tomography, this method represents an efficient alternative to timeconsuming velocity analysis and depth migration of seismic images. • Joint interpretation with other studies (Kley 2018;Ahlrichs et al. 2021;Rasmussen 2009) indicates that a significant amount of the uplift and erosion may have taken place during the Neogene. of the DFG-funded "StrucFlow" project (396852626). We would like to thank Jan Walda for access to and help with the Machine Learning algorithm. Developing the code was supported by the Federal Ministry for Economic Affairs and Energy of Germany, project number 03SX504A. Further, we thank Ari Tryggvason and Yaocen Pan for access and support with the PStomo_eq travel-time tomography. We gratefully acknowledge Schlumberger and IHS Markit for providing the 2016 VISTA seismic processing software and the Kingdom seismic interpretation software, respectively.
Author contributions AW and CH designed the study concept. AW is the corresponding author and carried out the travel-time picking, velocity model calculation and exhumation estimation. NA performed the seismic processing and the stratigraphic integration. All authors were part of the interpretation and discussion. AW created the figures and prepared the manuscript. All authors contributed to the editing and enhancement of the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. Data acquisition during RV MARIA S. MERIAN expedition MSM52 was financed by the German Research Foundation (DFG) together with the Federal Ministry of Education and Research (BMBF). This study was conducted in the course of the DFG-funded "Struc-Flow" project (396852626).

Data availability
The processed seismic data is available via the pan-European data portal Geo-Seas, https://www.geo-seas.eu/search/welcome.php?query=1678&query_code={558D29F1-65F5-4441-8CAF-8E15926546E9}. Additional raw data is available from the Federal Institute for Geosciences and Natural Resources (BGR) upon reasonable request. The Machine Learning code can be obtained from the corresponding author upon reasonable request for research purposes in Germany.

Conflict of interest
The authors declare that they have no competing interests.
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 http:// creat iveco mmons. org/ licen ses/ by/4. 0/.