Potential links between Baltic Sea submarine terraces and groundwater seeping

Submarine groundwater discharge (SGD) influences ocean chemistry, circulation, and the spreading of nutrients and pollutants; it also shapes sea floor morphology. In the Baltic Sea, SGD was linked to the development of terraces and semicircular depressions mapped in an area of the southern Stockholm archipelago, Sweden, in the 1990s. We mapped additional parts of the Stockholm archipelago, areas in Blekinge, southern Sweden, and southern Finland using high-resolution multibeam sonars and sub-bottom profilers to investigate if the sea floor morphological features discovered in the 1990s are widespread and to further address the hypothesis linking their formation to SGD. Sediment coring and sea floor photography conducted with a remotely operated vehicle (ROV) and divers add additional information to the geophysical mapping results. We find that terraces, with general bathymetric expressions of about 1 m and lateral extents of sometimes >100 m, are widespread in the surveyed areas of the Baltic Sea and are consistently formed in glacial clay. Semicircular depressions, however, are only found in a limited part of a surveyed area east of the island of Askö, southern Stockholm archipelago. While submarine terraces can be produced by several processes, we interpret our results to be in support of the basic hypothesis of terrace formation initially proposed in the 1990s; i.e. groundwater flows through siltier, more permeable layers in glacial clay to discharge at the sea floor, leading to the formation of a sharp terrace when the clay layers above seepage zones are undermined enough to collapse. By linking the terraces to a specific geologic setting, our study further refines the formation hypothesis and thereby forms the foundation for a future assessment of SGD in the Baltic Sea that may use marine geological mapping as a starting point. We propose that SGD through the submarine sea floor terraces is plausible and could be intermittent and linked to periods of higher groundwater levels, implying that to quantify the contribution of freshwater to the Baltic Sea through this potential mechanism, more complex hydrogeological studies are required.

1 Introduction

The influence of groundwater on sea floor morphology has been discussed for more than 80 years within the geoscientific community (Robb, 1990). Early examples include Stetson (1936), who suggested that the formation of submarine canyons in the southern flank of Georges Bank off the northeastern US coast was related to groundwater seeping aided by currents. Submarine canyons are today linked to a combination of geological processes including erosion by turbidity currents, slumping, and mass wasting (Harris and Whiteway, 2011; Shepard, 1981). However, groundwater as a shaping agent of sea floor morphology appears in several other marine geological settings (Robb, 1990). For example, mapped terraced walls and irregular courses of valleys in the continental slope off New Jersey have been interpreted to be formed by groundwater seeping during periods of lower sea level (Robb, 1984). Other examples of sea floor morphological expressions related to groundwater are depressions formed around submarine fresh and/or brackish springs in the Mediterranean, where extensive nearshore carbonate formations caused the development of submarine karstic aquifers (Rousakis et al., 2014). Depressions, which are similar in appearance to pockmarks formed by gas seeps, have also been found where freshwater escapes through the sea floor (Whiticar and Werner, 1981; Khandriche and Werner, 1995; Virtasalo et al., 2019). Groundwater discharge into the ocean is recognized as a widely occurring process and is commonly referred to as submarine groundwater discharge (SGD) (Moore, 2010). SGD is estimated to contribute about 6 %–7 % of the total hydrological discharge to the world oceans (Zektser, 2000), although different quantification methods have yielded varying results at specific localities (Prieto and Destouni, 2011). In the Mediterranean Sea, SGD has been shown to be a major source of nutrients (Rodellas et al., 2015) and the process has generally been raised as a potentially underestimated provider of chemical elements, including pollutants and nutrients, from land to coastal waters (Destouni et al., 2008).

In Sweden, terraces and semicircular depressions in the Baltic Sea floor were mapped in the 1990s along some islands in the southern Stockholm archipelago (Söderberg and Flodén, 1997) (Fig. 1). These features were interpreted to be formed by processes related to SGD (Söderberg and Flodén, 1997). The proposed mechanism was groundwater flowing through siltier permeable layers in glacial varved clay, and where the clay outcropped at the sea floor, the escaping fresh water led to the development of a terrace when the above layers of clay were undermined and collapsed. With respect to the semicircular depressions, Söderberg and Flodén (1997) found morphological similarities with sea floor features mapped near Kiel in Eckernförde Bay, southwestern Baltic Sea (Khandriche and Werner, 1995), where SGD is documented to occur from geochemical analyses of sea floor sediment and the water column (Schlüter et al., 2004).

Figure 1 Overview map showing the locations of the studied areas in this work with white dots. Corresponding figures for each area are shown with white text. The locations of previous studies showing SGD discussed here are shown with black dots. I and II: Söderberg and Flodén (1997); III: Schlüter et al. (2004); IV: Virtasalo et al. (2019). The bathymetry is from EMODnet 2018 (EMODnet Bathymetry Consortium, 2018)

More recently, about 1 m high sea floor terraces, extending from a few metres to >100 m in length, and semicircular depressions 10–30 m wide and about 1 m deep were mapped using a high-resolution multibeam echo sounder and sub-bottom profiler east of the island of Askö in the southern Stockholm archipelago (Jakobsson et al., 2016). Both the terraces and depressions closely fit the descriptions by Söderberg and Flodén (1997), suggesting that these kind of sea floor features may be widespread along the Swedish coast. Here we show from extended high-resolution geophysical surveys in the Stockholm archipelago and in Blekinge, southern Sweden, as well as in southern Finland, that sea floor morphological features in the form of terraces and semicircular depressions are widespread in the Baltic Sea. Sediment cores and bottom inspection with ROVs and divers reveal that the terraces consistently form in nearly identical geological settings, permitting the prediction of their occurrence if information on the sea floor geology is available. In order to further test the formation hypothesis involving groundwater seeping, we performed geochemical analyses of a carbonate concretion found in the glacial clay unit of a submarine terrace. This study provides a framework for continued investigations involving in situ observations of potential groundwater seeping at selected terraces and semicircular depressions along Baltic Sea coasts. The widespread occurrence of these features suggests that SGD in the Baltic Sea could potentially be an important influence on the chemical composition of its waters with implications for circulation and the spread of nutrients and pollutants.

2 Material and methods

2.1 Survey areas

Results are presented from four different regions in the Baltic Sea: (1) east of the island of Askö in the southern Stockholm archipelago; (2) west of the island of Kastellholmen in Stockholm Harbour; (3) the southern Blekinge archipelago; and (4) southeast of Tvärminne Zoological Station in the southern Finland archipelago (Fig. 1). Only a small subset of the acquired multibeam bathymetry is shown from the surveyed area in the southern Blekinge archipelago, and it is without geographic coordinates due to military restrictions of revealing detailed depth information. Geophysical mapping results from all other surveyed areas in this study are granted permission to be shown in their full extent and resolution.

2.2 Geophysical mapping, processing, and analyses

The high-resolution multibeam bathymetric data shown here were acquired using Stockholm University research vessels RV Electra (24.3×7.2 m) and R/V Skidbladner (6.4×2.4 m). A Kongsberg EM2040 0.4^∘ × 0.7^∘, 200–400 kHz, multibeam echo sounder is hull-mounted on R/V Electra, while R/V Skidbladner has a pole-mounted Kongsberg EM2040 0.7^∘ × 0.7^∘, 200–400 kHz multibeam. Both R/V Electra and R/V Skidbladner receive position, heading, and attitude data from Kongsberg Seatex Seapath 330+ navigation units with attached MRU5+ motion and reference sensors. The Seapath 330+ uses both GPS and GLONASS satellites and is capable of applying real-time kinematic (RTK) corrections, which for all surveys in Swedish waters were received from SWEPOS (https://swepos.lantmateriet.se/, last access: 20 December 2019) over the internet. In Finland, RTK corrections were provided over the internet by the Finnish equivalent, FinnRef (https://www.maanmittauslaitos.fi/en/node/1881, last access: 20 December 2019). The Seapath systems on both vessels generally indicate positions with horizontal accuracies better than ±5 cm and slightly coarser vertical accuracies. A Valeport MiniSVP sound velocity probe was used at discrete stations to acquire sound speed profiles in all surveyed areas for processing the multibeam bathymetry. In addition, both R/V Electra and R/V Skidbladner have AML Oceanographic^TM sound velocity probes mounted near the transducers for continuous reading of surface sound speed.

Sub-bottom profiles were acquired with R/V Electra using the hull-mounted Kongsberg Topas PS40 24-channel parametric sub-bottom profiler operating with primary and secondary frequencies of 35–45 and 1–10 kHz, respectively. The sub-bottom profile shown from the eastern part of the surveyed area east of the island of Askö was acquired using the system in chirp mode with a 4 ms long 4–10 kHz pulse. This system receives positions from the Seapath 330+. The second sub-bottom profile shown from the Askö area, located along the island, was collected in 2009 with a towed Edgetech SB-216S chirp sonar operated with a 1 ms long 3–9 kHz pulse. Positions were received from a Hemisphere VS100 GPS augmented with a satellite-based augmentation system (SBAS), yielding horizontal accuracies better than ±2 m. However, since the SB-216S was towed at about 3 m of water depth, georeferencing of the profile is considerably less accurate. The sub-bottom profile shown from Stockholm Harbour was acquired using a Kongsberg EA640 and a towed surface device; an Air15-17 transducer with a centre frequency of 15 kHz was mounted 20 cm below the water surface. A 13–17 kHz, 1.024 ms long chirp pulse was used. Positions were received from a Hemisphere VS100 GPS with the antenna mounted on the towing device directly above the transducer to avoid offsets. There were no sub-bottom profiles collected from the area in Blekinge, southern Sweden, surveyed with R/V Skidbladner in 2012. The terraces shown from the area east of Tvärminne in southern Finland were not captured properly with sub-bottom profiles during the two R/V Electra field campaigns 2017 and 2018 because they are located close to shoals, which prevented a profile from being captured perpendicular to the terrace extensions.

In 2013, we towed a Klein 3000 100∕500 kHz side-scan sonar over a stretch of sea floor with identified terraces east of the island of Askö. Positions were received from a Hemisphere VS100 GPS augmented with SBAS, i.e. a navigation setup similar to the SB-216S sub-bottom profiles described above.

Post-processing of the multibeam bathymetry was done using the QIMERA software by QPS version 1.7.2 (https://www.qps.nl/, last access: 20 December 2019). Bathymetric grids were produced in QIMERA with resolutions ranging from 1×1 to 0.25×0.25 m. These grids were subsequently imported to the open-source geographic information system QGIS, version 3.6.0 Noosa (QGIS Development Team, 2018), for further analyses and map making. The sub-bottom profiles were analysed using the open-source software OpendTect version 6.2.0 by dGB Earth Sciences (https://www.dgbes.com/, last access: 20 December 2019). However, the sub-bottom profile images shown in this article were produced using software tools provided by the Geological Survey of Canada, courtesy of Bob Courtney. Seeps, sometimes called “seep flares”, were identified in the multibeam midwater data logged in the area east of Askö using the FMMidwater suite of Fledermaus software. The “flare hunter” plug-in was used to batch-process the data from individual days, and identified flares were manually spot-checked and subsequently exported as a series of georeferenced points for the production of maps in QGIS.

2.3 Bottom inspection and photography

In the Askö and Tvärminne areas, terraces were photographed and filmed using two different remotely operated vehicles (ROVs): (1)a Saab Seaeye Falcon and (2) BlueROV2 by Blue Robotics. The photos shown from Stockholm Harbour and Blekinge were captured by divers. A 3-D model over an approximately 6 m long stretch of the terrace in Stockholm Harbour was assembled using the Agisoft Metashape software (https://www.agisoft.com/, last access: 20 December 2019).

2.4 Sediment coring, measurements, and analyses

Numerous sediment cores have been retrieved east of Askö from the areas with terraces during the field component of a Stockholm University course in marine geophysical mapping. This course used Askö yearly for fieldwork from 2009 to 2018. In Tvärminne, the mapped terraces were cored in 2017 and 2018. We present analyses of core Asko2018HT-2GC, retrieved at 58^∘50.76^′ N, 17^∘41.94^′ E in 16 m of water depth east of Askö, because it illustrates the characteristic sediment of the terraces we cored. The core was taken with Stockholm University piston–gravity coring system adapted for R/V Electra. This corer is capable of handling a core head with a maximum weight of 563 kg. The 6 m long core barrel is loaded with PVC liners having an outer to inner diameter of 110∕98.5 mm. Core Asko2018HT-2GC was retrieved with the corer rigged in gravity mode.

2.5 Sediment physical properties

Core Asko2018HT-2GC was subjected to high-resolution (1 cm) logging of sediment physical properties, including bulk density, magnetic susceptibility, and p-wave velocity, using a Geotek Multi-Sensor Core Logger (MSCL). The 4.3 m long core was cut into 1.5 m long sections that were logged before being split, described, and digitally imaged. Sediment grain size (0.2 µm–2 mm) was measured on 29 samples using a Malvern Mastersizer 3000 laser diffraction particle size analyser. Wet samples were immersed in a dispersing agent (<10 % sodium hexametaphosphate solution) and placed in an ultrasonic bath to facilitate particle disaggregation before analyses. Particle classes were defined using the European Standard EN ISO 14688-1:2018, according to which clay is the 0–2 µm fraction, silt 2–63 µm, and sand >63 µm.

2.6 Geochemical analyses

A carbonate concretion found in core Asko2018HT-2GC was subjected to geochemical analyses at the Stable Isotope Lab of the Department of Geological Sciences, Stockholm University, for determining the water source(s) during its formation. Two pieces of the concretion were analysed: one from the centre and one from the outer edge. The samples were milled to a powder and analysed for total carbon using a Thermo Delta V mass spectrometer and for phosphoric acid reactive carbon with a Gasbench II MAT 253 mass spectrometer. The elemental analysis provides δ¹³C of the total carbon, while the Gasbench II provides δ¹³C and δ¹⁸O of the carbonate. Paired t tests were applied to determine if the two methods provided the same δ¹³C results, while t tests were used to test for significant differences between the results from analysing samples from the inner and outer parts of the concretion. A confidence interval of 95 % was applied in all t tests. The δ¹³C and δ¹⁸O are reported with respect to the Vienna Pee Dee Belemnite (VPDB) standard (Coplen, 1996).

3 Results

3.1 Geophysical mapping

Terraces in the sea floor were mapped by Jakobsson et al. (2016) using a multibeam echo sounder in the vicinity of the island of Skåren, located east of Askö (Fig. 2). The terraces were identified in shaded relief images of the processed multibeam bathymetry. These terraces have steps of about 1 m high, extend from a few metres to >100 m in length, and are most abundant in water depths shallower than 15 m. In addition, a few semicircular depressions were mapped that resemble pockmarks with about 1 m of bathymetric expression, although with an opening towards the down-sloping sea floor (Jakobsson et al., 2016). Their semimajor axes are between 10 and 30 m. Here we have expanded the mapped area east of Askö considerably eastward and northward and derived descriptive statistics of the terrace depth distribution (Fig. 2). Most terraces are located in 12 m of water depth, while the mean and median depths are 16 and 15 m, respectively. The shallowest depth of a mapped terrace is 4 m, while the deepest is 28 m. However, it should be noted that our mapping did not extend much shallower than about 4 m of water depth. Many of the terraces extend laterally in rather systematic winding sinusoidal patterns, although there are some that take the form of a single hyperbola (Fig. 2c and d). The height of the terrace steps varies but is generally <1 m (Fig. 2e).

Figure 2 Multibeam bathymetry of the studied area east of the island of Askö, southern Stockholm archipelago. (a) Overview map with brown lines representing digitized terraces in the sea floor and blue dots showing locations where seeps were identified in the multibeam water column data. Black heavier lines show the locations of the sub-bottom profiles in Fig. 5a (SBP1) and b (SBP2). The sub-bottom profiles in Fig. 5 are portrayed between the two dots, with the profile running from left (black dot) to right (white dot). The location of core Asko2018HT-2GC is marked by a yellow star (abbreviation 02GC used on the map). (b) Depth statistics of the occurrence of terraces and seeps. (c) Detailed view of the terraces in the sea floor around the island of Skåren. (d) The location of the bathymetric profile across two terraces from X to X^′ shown in (e).

The multibeam bathymetry shows that the sea floor morphology of the terraces in the area east of Askö, Tvärminne, and the southern Blekinge archipelago closely resemble one another (Figs. 3 and 4). There are only a couple of mapped terraces in the latter two areas preventing a meaningful statistical comparison of their depth distribution. However, the terraces we mapped outside the area of Askö occur in the deeper depth range: (1) Tvärminne at 15–20 m and (2) the Blekinge southern archipelago at 23–25 m. The terrace in Blekinge is the longest we mapped; it is possible to trace continuously for >550 m (Fig. 4a). We do not have multibeam bathymetry of the terrace found in Stockholm Harbour, preventing a comparison of its morphology and spatial extent to the terraces in the other three areas.

Figure 3 Multibeam bathymetry of the studied area east of Tvärminne Zoological Station in the southern Finland archipelago. (a) Overview of the sea floor bathymetry in the vicinity of the only set of terraces mapped in this area, shown in (b). (c) Bathymetric profile from Y to Y^′ across two terraces. The profile location is shown in (c).

Figure 4 Multibeam bathymetry of a terrace mapped in the southern Blekinge archipelago. (a) Only a smaller part of the full survey is shown without coordinates because detailed bathymetric information from this area is under military restrictions. (b) Detail of the >500 m long terrace with the location of the bathymetric profile from $Z-Z^{\prime }$ shown in (c).

Multibeam water column information was logged and analysed for two-thirds of the surveyed area east of Askö. Seeps from the sea floor were found to be a common feature (Fig. 2a), and the question immediately arose if the seeps were related to either terraces or depressions in the sea floor. There is an abundance of seeps in the northern part where no terraces are identified. East and southeast of the island of Skåren, seeps begin to occur at about 20 m of water depth; i.e. with a few exceptions this is where the deepest terraces occur (Fig. 2b). The mean and median depths of the seeps are 21 and 20 m, respectively, while the shallowest is located at a water depth of 3 m and the deepest at 40 m. The multibeam backscatter shows that the terraces in the area east of Askö systematically appear in a relatively harder sea floor characterized by high backscatter, while the seeps generally occur in a softer seabed represented by lower backscatter (Fig. S1 in the Supplement). The side-scan data acquired in 2013 show high-resolution imagery of the terraces, with no apparent difference in signal intensity across them but with clear shadows present due to the bathymetric expressions (Fig. S2).

Figure 5 Sub-bottom profiles and bottom photographs portraying sea floor terraces. (a, b) Sub-bottom profiles SBP1 (a) and SBP2 (b) from the area east of Askö. The locations of terraces clearly identified in the multibeam bathymetry are marked T1–T9 (see Fig. 2a). The site of core Asko2018HT-2GC (02GC on profile SBP1) is shown with a black arrow. Note that it is located on top of terrace T1. (c) Sub-bottom profile across the terrace mapped west of the island of Kastellholmen in Stockholm Harbour. Bottom photographs are shown in (d, e) and a 3-D photomosaic of this terrace in (f).

There is a semi-regular grid of sub-bottom profiles covering the entire area east of Askö (Fig. S3). These profiles were acquired during the field component of a Stockholm University course in marine geophysical mapping held yearly in this area since 2009. From this database two sub-bottom profiles, extending across terraces identified in the multibeam bathymetry, are shown in Fig. 5a and b. The terrace west of the island of Kastellholmen in Stockholm Harbour was mapped by a sub-bottom profile perpendicular to its extent (Fig. 5c). Common for all terraces imaged by sub-bottom profiles is that the acoustic stratigraphy indicates well-stratified sediments that outcrop at the sea floor where the terrace is formed. This is particularly clear in the profile from Stockholm Harbour (Fig. 5c). In the area east of Askö, where the sub-bottom profile coverage is most comprehensive, an acoustically semitransparent surface unit with few internal noticeable reflections is commonly found in sections with water depths deeper than the general occurrence of terraces (Fig. 5b).

Figure 6 Bottom photographs of sea floor terraces in the studied areas. (a, b) East of the island of Askö, (c, d) east of Tvärminne Zoological Station in the southern Finland archipelago, (e, f) near Boön, Blekinge archipelago.

3.2 Bottom photography

The bottom photos confirm that distinct terraces are formed in the sea floor in all mapped areas (Figs. 5d–f and 6a–f). It is also possible to identify that the terraces developed in stratified sediments outcropping at the sea floor, which was particularly evident in the acoustic stratigraphy of the sub-bottom profile from Stockholm Harbour (compare Fig. 5c with Fig. 5e and f). Holes with diameters between about 1 and 2 cm are abundant in the nearly vertical terrace walls in Stockholm Harbour and Tvärminne (Figs. 5d, e and 6c, d). Some of these holes appear to be cavities from stones that were embedded in the sediments and eventually fell out during the erosional process forming the terrace. However, we cannot confirm if this is always the case because the holes sometimes appear to extend rather deep into the terrace walls. Hence, they may be zones of piping and erosion that developed in response to focused groundwater flow.

3.3 Sediment stratigraphy and physical properties

The 4.3 m long core Asko2018HT-2GC, retrieved from a terrace in a water depth of 16 m in the area east of Askö, consists of rhythmically alternating 0.5–2.5 cm thick silty–clay layers (Fig. 7). The upper 3.5 m is composed of inclined (∼10–16^∘) rhythmites with a gradual transition into the lower 0.8 m thick interval with nearly horizontal rhythmites. Micro-faults offsetting some rhythmites by ∼0.5 cm appear in the core from about 1.5 m of core depth and are present to the bottom of the core. The silt content within the 0.5–2.5 cm thick rhythmites varies as does the colour, which changes from greyish brown (Munsell: 5/2 10YR) to reddish brown (Munsell: 5/3 5YR). Over the interval at which grain size measurements were obtained (2.80–4.20 m b.s.f.) the sediments are composed almost exclusively of fine-grained material (>95 % is <20 µm) (Fig. 8). At about 3.5 m, p-wave velocity, bulk density, and magnetic susceptibility all transition towards slightly higher values, with peak values occurring in zones slightly enriched in silt (Fig. 7). Sediment porosity, calculated from the bulk density logs by assuming a grain density of 2.71 g cm⁻³, also decreases from ∼60 %–70 % to ∼55 %–60 % below 3.5 m. No disconformity is recognized across this transition. The porosity reduction does not appear to be driven by normal down-hole compaction but may be related to the slight coarsening of the sediments.

Figure 7 Lithology and sedimentology of sediment core Asko2018HT-2GC. Image includes a lithologic log (a), sediment physical properties measured on the MSCL (b), and results from grain size analyses showing variations in clay and silt content across the zone where the overall porosity reduction is seen (∼3.5 m) (c). Characteristic images from four different depths are shown (d), including the 2 cm diameter concretion recovered at 2.64 m and the smaller concretion at 1.77 m.

Table 1 Geochemical analyses of the concretion found in the glacial clay unit of core Asko2018HT-2GC. See Fig. 2a for core location and Fig. 7 for an image of the concretion. The upper section shows elemental analyses made using the Thermo Delta V mass spectrometer, and the lower analyses were made on phosphoric reactive carbon with a Gasbench II MAT 253 mass spectrometer (see “Material and methods” for further details).

Figure 8 Cumulative grain size distributions for the laminated glacial clays in Asko2018HT-2GC. Sediments are exclusively composed of fine fraction material, with 50 %–80 % of the sediments being <4 µm.

Two carbonate concretions were found at core depths of 1.77 m (diameter 1 cm) and 2.64 m (diameter 5.3 cm) (Fig. 7). The lowermost one has a disc-shaped appearance, displaying a series of concentric rings growing outwards from a central spherical concretion. It is a classic marleka, also known as a “fairy-stone” or “imatra stone”, which forms in Quaternary clay-rich sediments as calcium carbonate precipitates around a small pebble or organic matter (Neuendorf et al., 2011).

3.4 Geochemistry

The results from the geochemical analyses on the concretion are shown in Table 1 (δ¹³C and δ¹⁸O are reported in per mille vs. VPDB). The elemental analyses of the outer edge of the concretion give an average δ¹³C of −20.23 ‰ (SD = 0.11), while the Gasbench II yields an average of −19.87 ‰ (SD = 0.02). A paired t test (p=0.04) indicates a significant difference between the two methods for the outer edge samples considering a 95 % confidence interval (α=0.05). From the central sample, averages of −19.35 ‰ (SD = 22) and −19.32 ‰ (SD = 0.03) were provided by the elemental analyses and Gasbench II methods, respectively. A paired t test (p=0.88) suggests that there is no significant difference between the measured centre samples using the two different approaches. If we compare the average of all measurements made on the samples from the edge of the concretion ($M=-\mathrm{20.05}$ SD = 0.21) with the average of all measurements made on the samples from the centre ($M=-\mathrm{19.33}$ SD = 0.14), a t test suggests a significant difference (p=0.00) in δ¹³C, although it is small. The average δ¹⁸O from the central sample is −10.54 ‰ (SD = 0.19) and the edge −8.02 ‰ (SD = 0.14), with a t test suggesting a significant difference (p=0.00) between the central and outer edge of the concretion with respect to δ¹⁸O. The total carbon has an average of 6.98 % (SD = 0.36) at the edge and 10.46 % (SD = 0.32) in the centre. A t test (p=0.00) demonstrates that these means are significantly different from each other.

Figure 9 Schematic illustration showing how siltier layers in glacial clay could act as a conduit for groundwater, eventually escaping at the sea floor as SGD and leading to the formation of a bathymetric terrace. The critical parameters dh and dL in Darcy's law (Eq. 2) are illustrated in the sketch as they dictate the hydraulic head that would drive the flow.

4 Discussion and conclusions

The first discoveries of terraces formed in the sea floor along the Swedish coast were made in the 1990s with a conventional 30 kHz echo sounder and a 100∕500 kHz side-scan sonar (Fig. 1; Söderberg and Flodén, 1997). Our detailed surveys east of the island of Askö showed that for a terrace to be identified with a single-beam echo sounder, the profiling direction cannot deviate much from being perpendicular to the orientation of the terrace. The main reason for this is that the bathymetric expressions of the terraces rarely exceed 1 m. They are easier to map with side-scan sonar, although as no bathymetric information is provided by a conventional side-scan sonar, they may also be misinterpreted as other surface patterns in the sea floor (Fig. S2). It is only with the latest generation of shallow-water high-resolution multibeam echo sounders that the terraces are irrefutably recognized as prominent bathymetric features in the sea floor (Figs. 2–4). This may explain why it took around 2 decades to realize how common the type of terraces mapped by Söderberg and Flodén (1997) are along the Swedish coast, as well as elsewhere along the coasts of the Baltic Sea. Recognizing that the features are widespread, the hypothesized formation mechanism involving SGD makes them important from an environmental point of view. The terraces could comprise focal points at which groundwater enters the Baltic Sea to influence its brackish waters by providing not only freshwater, but also potentially pollutants and nutrients. If this is the case, quantification and chemical analyses of the SGD from the terraces will provide new information for assessments of the nutrient budget of the Baltic Sea.

Our study generally supports the hypothesis put forward by Söderberg and Flodén (1997) on the terrace formation mechanism. They argued that groundwater flows through siltier layers in glacial clay to eventually escape at the sea floor. Erosion from the flowing water undermines the overlying layers, causing them to collapse and form a sharp terrace in the sea floor (Fig. 9). The underwater photos in Fig. 6 include examples in which we believe that this process can be readily envisioned. It is possible to see how cavities are formed in the varved clay at the bottom of some of the terraces as well as blocks of the overlying clay that have collapsed from being undermined. However, we cannot exclude the possibility that mechanisms other than SGD could produce terraces in the sea floor similar to those we mapped in this study, and therefore alternative formation mechanisms are discussed below.

We have also found from our geophysical mapping and coring results that the terraces are systematically formed in glacial clay throughout the studied areas. Glacial clay in the Baltic Sea sediment stratigraphy is commonly found draped on top of till or glaciofluvial material; in some cases it rests directly on bedrock (Andrén et al., 2011, 2015) (Fig. 9). This type of clay was for the most part deposited during the last deglaciation in front of the retreating Scandinavian Ice Sheet. Glacial clays left from previous older glaciations are extremely rare and found only at a few locations (e.g. Björck et al., 1990). Swedish geologist Gerard De Geer discovered that glacial clay is comprised of rhythmites, wherein the layers composed of varying proportions of clay and silt are annual depositions of erosional material from the retreating ice sheet (De Geer, 1912). He introduced the term “varve” for one annual layer of glacial clay and noted that its thickness and silt content varied depending on the proximity to the retreating ice margin; i.e. thicker and siltier varves were deposited close to the ice margin. He further proposed that there would be a higher degree of silt content in the part of the varve representing the meltwater-rich summer period. Thus, grain size variations in the glacial clays are found on a number of scales: from millimetre- to centimetre-scale variations across rhythmites, to longer-scale (decimetre) variations related to climatically driven variations in subglacial discharge, to even longer-scale (> metre) variations related to the proximity of the ice margin. The use of glacial varved clay as a record documenting the ice retreat was also adopted early on the Finnish side of the Baltic basin (Sauramo, 1926). From the knowledge gained from these studies, it follows that there is a higher chance of finding more silt-rich layers in the older sections of glacial clay deposited close to the ice margin. This would be in the lower sections of the glacial clay units. In the case of Asko2018HT-2GC, this larger-scale variation appears to be captured in the coarsening of the grain size below 3.5 m of depth, which has length scales that exceed the duration of individual rhythmites (Fig. 7). Although fine-grained glacial clays are not commonly considered to be highly transmissible sediments, variations in the silt content of clay-rich sediments have a dominant effect on their permeability (Schneider et al., 2011). Experimentally, it has been shown that at a given porosity for silt–clay mixtures, increasing the clay content (<2 µm) from 36 % to 57 % decreased the permeability by an order of magnitude (Schneider et al., 2011).

It should be noted that the glacial clay sequences are time-transgressive throughout the Baltic basin, with older clay in the south and younger in the north because the Baltic basin was deglaciated from south to north (Hughes et al., 2016; Stroeven et al., 2016). The first small freshwater body in which glacial clay could be deposited during the last deglaciation formed in front of the ice margin around eastern Denmark and the northern coasts of Germany and Poland at about 16–15 ka (Houmark-Nielsen and Henrik Kjær, 2003). This water body grew as the ice sheet retreated northward to become the Baltic Ice Lake, which appears to have lasted until the end of the Younger Dryas cold period at about 11.7 ka, when it catastrophically drained westward north of Mount Billingen in south–central Sweden (Björck and Digerfeldt, 1986; Andrén et al., 2002; Swärd et al., 2015). Mapping of preserved palaeo-shorelines in the 1920s showed that the drainage occurred when the ice sheet margin reached north of the damming high terrain in the west (Lunqvist, 1921). A brackish water phase called the Yoldia Sea, likely constrained to the central Swedish side of the Baltic (Schoning, 2001), followed the Baltic Ice Lake (Björck, 1995); however, it would take several hundred years for the Baltic to become brackish after the drainage and deposition of varved glacial clay continued close to the retreating ice margin (Andrén et al., 2011). The Baltic Sea basin was completely ice free at about 10 ka (Hughes et al., 2016). As the ice retreated (Stroeven et al., 2016), conditions may have developed for terrace formation at different places around the Baltic Sea depending on local sea level in relation to glacial clay deposits. The mapped terraces in the different regions may therefore be of different ages; some may be inactive, while others are active today.

The geochemical analyses of the concretion from core Asko2018HT-2GC do not irrefutably determine whether or not groundwater flow occurred through the glacial clays, but they provide valuable insights into the formation environment. The δ¹³C carbonate isotope values of the concretion between about −19 and −20 ‰ are low compared to Baltic Sea δ¹³C DIC (dissolved inorganic carbon), which is usually between 0 and 1 ‰ (Filipsson et al., 2017). Carbonates formed from Baltic Sea water should thus have similar values. Therefore, the most plausible explanation for the observed isotope values is that the carbon source is respired organic matter in the sediments. Organic matter in the Baltic Sea is usually in the range −23 to −28 ‰ depending on if it is coming from terrestrial sources or from primary production in the Baltic Sea (Alling et al., 2008; Deutsch et al., 2012). The small difference in isotope values between the concretion and organic matter in the Baltic Sea could be due to mixing with DIC from Baltic Sea water during formation.

The δ¹⁸O values of samples taken from the edge and the centre of the concretion are significantly different from one another by about 2.5 ‰. Here we investigate possible causes for the isotopic differences between these samples using the Shackleton (1974) equation for relating temperature to δ¹⁸O_c and δ¹⁸O_w:

$$\begin{array}{}\text{(1)}& t=\mathrm{16.9}-\mathrm{4.0}\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{w}}\right),\end{array}$$

where t is temperature, O_c is from carbon dioxide extracted from carbonate, and O_w is from carbon dioxide equilibrated with water. From Eq. (1) it follows that the δ¹⁸O_w could be different in the samples from the edge compared to the samples from the middle if the formation temperatures were different. However, the difference between the middle and the centre is too large to be explained only by a change in temperature during the formation. It more likely indicates a change in the composition of the water source, perhaps in combination with a temperature change.

Precipitation in the Stockholm area has a yearly mean δ¹⁸O of about −10 ‰ (the Global Network of Isotopes in Precipitation (GNIP) database, IAEA/WMO). The Baltic Sea currently has a δ¹⁸O value between about −8 and −9 ‰ at the salinity where the core is taken, east of the island of Askö (Deutsch et al., 2012). However, the Baltic Sea stages during which glacial clay was deposited most likely had much lower values considering the influence from the meltwater of the Scandinavian Ice Sheet. For example, isotopic δ¹⁸O values of the Greenland Ice Sheet are of the order of −35 to −40 ‰ (Andersen et al., 2004). This is in line with isotopic measurements made on pore water from a sediment core retrieved in Lake Vättern, which show a prominent down-core progression to δ¹⁸O values of less than −35 to −40 ‰ in the sections of the core representing the Baltic Ice Lake (Fig. S4). This 74 m long sediment record captured the Baltic Sea stages from the Baltic Ice Lake to the time Lake Vättern was isolated at 9530±50 cal yr (Swärd et al., 2018). The generally low δ¹⁸O values from the concretion (−7.93 to −10.70 ‰) are much less negative than modern Baltic Sea water and even further removed from the inferred δ¹⁸O composition of the Baltic Ice Lake ($<-\mathrm{35}$ to −40 ‰). This suggests that groundwater took part in the formation of the concretion because if it had been water from a Baltic Sea stage heavily influenced by meltwater, lower δ¹⁸O values should be expected.

Isostasy plays an important role in estimating the potential hydrogeological connection in the Baltic region between land and sea through glacial clay. The Swedish terrain was isostatically depressed by the Scandinavian Ice Sheet, which was several kilometres thick (Lambeck et al., 2010), and the highest coastline is therefore found in several areas far inland of the present coast (Björck, 1995). Isostatic rebound eventually caused glacial clay sequences to be lifted above the level of the Baltic Sea, thereby creating a hydraulic head between the landward end of permeable glacial clay layers and their seaward continuation (Fig. 9). To estimate SGD into the Baltic Sea through these permeable layers one could apply Darcy's law according to

$$\begin{array}{}\text{(2)}& Q=-KA{\displaystyle \frac{dh}{dL}},\end{array}$$

where Q is the water flow (m³ s⁻¹) in a saturated porous material, K is the hydraulic conductivity (m s⁻¹), A is cross section area (m²), and dh and dL represent the difference in height (pressure drop) with respect to the distance over which the water flows. Future estimates of the SGD rates into the Baltic Sea from permeable layers in glacial clays will require knowledge of where terraces occur on the sea floor, the hydraulic head, i.e. dh and dL, and the hydraulic conductivity K of the permeable siltier layers, which is governed by their composition. Such estimates will require high-resolution mapping efforts to identify potential locations of SGD, cross-sectional areas, and emission depth.

The SGD through the terraces will most likely be intermittent considering the relatively large variation of the seasonal groundwater table around the Baltic Sea, as shown by geological survey monitoring on land. Looking at data from monitoring stations provided by the Swedish Geological Survey in the vicinity of the Stockholm archipelago, we note that the seasonal variation in the groundwater table of some locations exceeds 3 m, and the highest groundwater table is generally found from the late fall to late spring (Fig. S5). Groundwater discharge to the Baltic Sea through the sea floor is known to occur through geological formations other than the terrace formations discussed here, in particular where glaciofluvial aquifers on land connect with the sea floor (Peltonen, 2002). However, the prevalence of the potentially SGD-related terraces in the mapped regions suggests that SGD to the Baltic Sea is likely underestimated and remains an unconstrained source of pollutants and nutrients, as previously argued by Destouni et al. (2008). The ecological influences of the terraces is therefore unknown. It is also clear that there is a general lack of data and understanding of SGD emission mechanisms for proper assessments of emission rates and/or volumes (Taniguchi et al., 2002).

There are other mechanisms that potentially could have played a role in the formation of the sea floor terraces mapped in this study. For example, sliding and slumping of glacial varved clays has been suggested to occur due to the liquefaction of layers during palaeoseismic events (Hutri and Kotilainen, 2007; Virtasalo et al., 2007). This could leave behind terraces at the sea floor formed in glacial clay. However, we do not observe any morphological evidence of sliding, and most of the terraces we mapped occur in areas where the sea floor slopes at <1^∘ and the terraces have nearly flat bases, as evident in the bathymetric profiles in Figs. 2e, 3c, and 4c. We also note that the terraces we mapped are widespread across the Baltic and systematically appear in glacial varved clay. It seems unlikely that slides would occur over such spatially large areas in several regions. Finally, the processes responsible for the formation of some of the terraces seems to be ongoing judging from the bottom photographs showing that small blocks of clay are presently falling down to form sharp terraces (Fig. 6). While we cannot exclude the possibility that other processes formed the terraces mapped in this study, we interpret our results to be in support of the formation mechanism proposed by Söderberg and Flodén (1997). Our study provides a geological and morphological framework for further research involving longer-term monitoring of potential SGD from the terraces.