Abstract
We investigate the landscape development of the early Mesolithic hunter-gatherer sites of Duvensee (10000–6500 cal. BCE). Based on ground-penetrating radar (GPR) and geoarchaeological drillings, we present for the first time a three-dimensional (3D) reconstruction of the palaeoenvironment of 63 ha covering subarea of the former lake during the Mesolithic. The archaeological aims were (1) to detect the location of former islands possibly hosting hunter-gatherer settlements and (2) to reconstruct the ancient landscape development for understanding prehistoric land use. The research in Duvensee lasts almost 100 years, providing vivid illustrations of early Mesolithic life. Clusters of Mesolithic camps have been found located on small sand hills that formed islands in the prehistoric lake. For this environment, we present depth maps of the three most important sedimentary facies interfaces of the ancient Lake Duvensee. Interface1 represents the transition between coarse organic sediments (peat and coarse detritus gyttja) and fine-grained organic sediments (fine detritus gyttja, calcareous gyttja), Interface2 represents the transition to the underlying clayish-loamy sediments, and Interface3 marks the top of the basal sand deposits at the lake bottom. From Interface3, we identified the location and extent of five former islands with Mesolithic camps. Stratigraphic information from the corings enabled us to create a 3D model of the spatio-temporal development of the Duvensee bog. The locations of the islands and their estimated dive-up times agree with the spatio-temporal pattern of the previous archaeological finds. The model shows where hunter-gatherers could settle and move from one island to another following the shorelines of the overgrowing lake. The 3D stratigraphic model provides growth and shrinking rates of the island and lake areas in the Mesolithic, and volumes of organic and non-organic deposited lake sediments. Besides, it provides a basis for a sustainable groundwater management needed for heritage preservation.
Introduction
The evolution of a landscape is important for understanding human behavior within it. Aspects like geomorphology and anthropogenic influences have to be taken into account for a full comprehension of palaeolandscapes. The conventional techniques to figure out the stratigraphy of sites and their surroundings include coring and near-surface geophysical surveys (e.g. De Smedt et al., 2011 or Gourry et al., 2003). Drilling allows a detailed registration of the vertical layering but usually with a low lateral sampling density due to the time-consuming character of this method. This limitation in the interpretation of lateral variation can be improved using non-invasive geophysical methods like ground-penetrating radar (GPR) which allow a more continuous mapping of the lithological change.
In particular, GPR has been extensively used in archaeological studies for mapping and imaging subsurface objects. The application of this method is based on detecting the contrast in electromagnetic properties (dielectric permittivity) and on recognizing the different signatures between man-made and natural features (Campana and Piro, 2008; Conyers, 2012; Zhao et al., 2013). GPR has been employed for the investigation of sedimentary geometries, stratigraphic units (e.g. Davis and Annan, 1989; Jol and Smith, 1991; Neal, 2004; Stern, 2008), facies analysis (e.g. Carreon-Freyre et al., 2003; Grant et al., 1998; Heinz and Aigner, 2003; Mellet, 1995; Pipan et al., 2000; Ruffell et al., 2004), and peatland stratigraphy (e.g. Doolittle and Butnor, 2009; Holden et al., 2002; Lowry et al., 2009; Menotti and O’Sullivan, 2013; Warner et al., 1990), but it has also some limitations. Clay layers and a high level of water saturation produce a strong attenuation of GPR signals, which especially restricts the application of this technique on alluvial deposits (Ferring, 2001; Moorman, 1990; Tolksdorf et al., 2013).
In this paper, we present a geophysical investigation of the Duvensee palaeolake in south-eastern Schleswig-Holstein, Germany, which is one of the most relevant micro-regions for early Mesolithic hunter-gatherer archaeology on the Northern European Plain (Groß et al., 2018). The former lake occupies a depressed area within late Pleistocene sandy moraines that was formed by melting dead-ice after the retreat of the last Weichselian glaciation. Its irregular topography resulted in several islands scattered along the ancient water body, which were used by early Mesolithic hunter-gatherers groups to establish temporary camps.
The current study is part of an interdisciplinary project (which is part of the CRC 1266 ‘scales of transformation’, Kiel University) that focuses on reconstructing the environmental evolution of ancient Lake Duvensee against the background of its archaeological importance through prehistory. The main goal of the geophysical investigation is to find the location of former islands hosting early Mesolithic camps, enabling an evaluation of their relevance impact in the human occupation. It is very likely that in a transforming environment, the hunter-gatherers changed the locations of their camps from an island to another, following the lakeshore of the overgrowing lake. Thus, recreating the situation during that period by means of GPR and geoarchaeology is our goal. Following the background of intensive research in the Northern European Plain (Groß, 2017; Groß et al., 2018, 2019; Hartz et al., 2014; Lübke et al., 2011; Schmölcke, 2016; Schmölcke and Nikulina, 2015), this paper provides a first reconstruction of the paleotopography and stratigraphy of parts of the bog during the early Mesolithic using geophysics.
After giving a geological and archaeological background of the investigated area, we describe the geophysical survey carried out with GPR, its interpretation, and how it contributes to build a new picture of the lake development connected to human occupation. In the last section, we propose a three-dimensional (3D) model of the lake using the depth information of the different facies from GPR and stratigraphic information from corings. This allows to trace the extension of five former islands where hunter-gatherers’ settlement has been discovered. Based on the topographic information of the lake bottom in ancient time, the accessibility of the different islands could be inferred, explaining how people moved across the islands during the early Mesolithic.
Geological setting
The landscape surrounding the Duvensee basin was modeled by late Pleistocene glacial activity (source: Geologische Übersichtskarte 1:200.000, www.bgr.bund.de). The basin itself is the result of melting dead-ice blocks following the retreat of the Weichselian glacier front. Past stratigraphic investigations revealed a complex basin topography, characterized by numerous deeps and sand banks (Averdieck, 1986). Lacustrine sedimentation began approximately in the Older Dryas, while maximum areal extension was achieved during the early Holocene (Averdieck, 1986; Günther, 1986). The gradual infilling of the basin was accompanied by a general decrease in water level over time. By the 19th century CE, most of the basin was occupied by a bog. Two remnant open-water areas were drained in 1850 CE to facilitate peat extraction and agricultural activities (Averdieck, 1986). A simplified model of the gradual sediment deposition in the lake area is shown in Appendix S1 (available online). The various stages of sedimentation show characteristic profiles in connection with the construction of Mesolithic camps in the organic layer.
Archaeological background
In 1923, the first archaeological finds in the 4.3–km2 large Duvensee bog were made by the geologist Karl CJ Gripp who found early Mesolithic artifacts in drainage ditches during an excursion (Schwantes et al., 1925, 1939). Subsequently, the archaeologist G Schwantes, started the first investigation on record concerning the area and identified four different sites which he named Wohnplatz (abbreviated WP), meaning dwelling site. The result of this first investigation was that only two of the discovered sites (WP 2 and WP 5) provided informative material in terms of flint assemblage and a paddle made of pine wood (Bokelmann, 1971; Schwantes et al., 1925). In 1946, H Schwabedissen (1951) excavated another site (WP 1) and further intensive research was done by K Bokelmann whose work started in 1966–1967 (Bokelmann, 1971, 1991) and lasted until 2001.
Nowadays 23 different locations of early Mesolithic hunter-gatherer groups and Neolithic farmers are known in the Duvensee bog and 17 of them have been partly excavated. While the duration of each camp was rather limited in time, the full extent of the Mesolithic occupation at Duvensee spans a few millennia and has witnessed several landscape transformations. Radiocarbon dating sets the establishment of the earliest known Mesolithic camps within the Preboreal time period (WP 8-9; ca. 9000 cal. BC), whereas the youngest Mesolithic evidence is dated to the early Atlantic period (WP 19; ca. 6500 cal. BC). In addition, one site belongs to the northern early Neolithic Funnel Beaker culture (WP 17) and one to the late Neolithic Single Grave culture (WP 15; Figure 1b).

Area of investigation including the archaeological and geophysical research. (a) Location of the Duvensee area considering also the maximum extent of the lake in the early Holocene (Groß et al., 2018). (b) Dating, location of the sites and the supposed positions of the islands (yellow dashed lines) in the western part of the ancient Lake Duvensee. (c) Ground-Penetrating Radar Survey (red lines), CMP measurements (green stars), and corings done in the area of interest (color-coded dots). (d) Locations of the corings done in 2016 and 2017 along the reference profile.
Coherent overviews of the research of ancient Lake Duvensee have been given by Bokelmann (2012) and Groß et al. (2018). All of the sites investigated so far are located on small sand banks that formed multiple islands (called islets by K Bokelmann) in the prehistoric lake. The archaeological layers are not directly located on the mineral soils, but within the overlaying organic sediments, because these islands started to be overgrown with peat already during the early Holocene. Island 1 (Figure 2b), for instance, is located closest to the western shore of the former lake. The known Boreal sites WP 1-6 as well as WP 11 and 12 are located at 150 m to the east on island 2, a sandy island, which was exposed by a lowering of the lake water level during the Boreal. A thin gyttja layer shows that it must have still been under water in the Preboreal (Bokelmann, 1971; Groß et al., 2018). Two other islands – island 3 and island 4 – are also known to be occupied during the Boreal and the early Atlantic periods. From the distribution of surface finds, at least two more islands (5 and 6) are expected (Figure 1a and b): island 5 in the south and island 6 in the north (the area including island 6 is not part of this study).

Velocity calculation. Example of two different CMP measurements and NMO velocity determination.
At a few sites, specific fireplaces were excavated, defining the main anthropogenic feature of a camp, which can be connected with hazelnut roasting and other subsistence activities (Figure 7 in Bokelmann, 2012). These hearths, however, were not constructed following a uniform design scheme, probably this may be connected with different purposes of the features or chronological differences (Bokelmann, 2012). In several cases, the base of the roasting hearths in Duvensee consisted of a loamy/clayey layer which was interpreted as remnants of specially constructed clay plates (Lage, 2011). But in some cases, a compact layer of loam directly on a layer of non-heated sand was detected (Bokelmann, 2012), giving the impression that the loam layer was not heated on gleaming charcoal. Generally, it has to be assumed that several hearths were used for hazelnut processing but different construction techniques are present. Nevertheless, this makes a detection of anthropogenic features of the camps very difficult in terms of geophysics because the single remains are below the resolution of any technique. We thus focused on the detection and investigation of the island features as a whole.
Methodology
GPR is a non-invasive geophysical technique that can be used to detect electrical discontinuities in the shallow subsurface. The theory and methodology of GPR are comprehensively explained, for example, in Davis and Annan (1989). The electrical properties (dielectric permittivity and electric conductivity) of the different sediments associated with changes in water and clay content determine the electromagnetic wave propagation velocity and the amount of energy which is absorbed and back-reflected at the interfaces between different lithologies (Davis and Annan, 1989; Doolittle et al., 2007). In a medium with uniform moisture content, the conversion of reflection traveltime to interface bottom can be done using one average velocity without introducing errors in depth determination. This condition is rarely encountered in the field and using one average value may cause errors (Boll et al., 1996). Identifying the correct velocity distribution of electromagnetic waves in the subsurface is essential for accuracy and precision of the depth of investigation in a way to locate an archaeological object, feature or layer, and crucial for any imaging method applied to the data. Knowledge of subsurface GPR velocity is also crucial for guiding intensive probing and excavations that may follow a GPR survey or, in our case, to ground truth the geophysical results with the stratigraphy coming from the corings.
Numerous sedimentological studies have used GPR to reconstruct past depositional environments, but often with open questions regarding the resolution of the investigation (Neal, 2004). The vertical resolution increases with the used GPR center frequency, enabling the detection of thinner layers. From theory, we know that the best structural resolution that can be reached is about one quarter of the dominant wavelength (Neal, 2004; Sheriff, 1977). Due to energy absorption, both signal strength and center frequency of radar signals decrease during wave propagation (e.g. Wunderlich and Rabbel, 2013), leading not only to a limitation of depth penetration but also to a decrease in structural resolution with depth. For sedimentological interpretation, this implies that the lower stratification can be poorly imaged compared to the upper stratification. If the subsurface incloses layers that are thin compared to the dominant wavelength, the accuracy of the stratigraphic interpretation of GPR sections can be greatly improved by incorporating information from drillings. In this context, drilling enables the recognition of major interfaces and facies patterns in radargram and serve thus for ground truthing.
GPR profiling and interpretation procedure
A GSSI GPR unit with a 200 MHz center frequency antenna was used to investigate the Duvensee area (for acquisition settings and processing sequence, see Appendix S2, available online). In monostatic mode, this antenna enables the acquisition of profiles with a constant, close to zero, offset between transmitter and receiver. A rectangular grid of GPR profiles was established to investigate the palaeolandscape with 30-m spacing in the northern part of the area. In the southern sector, the profile spacing is ca. 100 m in NS direction and 30 m in WE direction (Figure 1c). For deriving a 3D stratigraphic subsurface model, it was necessary to recognize and pick in every radargram the reflections associated with the interfaces between the main sedimentary units identified in the drill cores (yellow dots in Figure 1c): basal sands, fine minerogenic sediments, fine and coarse organic deposits. The description of each sediment unit will be presented in Section 5.1. Picking was performed using the Kingdom IHS Software for the entire survey area. This Software was chosen because it allows to display GPR data together with drillings. During the picking procedure, a two-dimensional (2D) and 3D view of the dataset is possible comparing anytime the GPR reflectors with the stratigraphic information for each picked profile.
Kingdom enabled then the comparison between the GPR profiles with the stratigraphic information from the corings in order to understand which interfaces and facies units could be identified uniquely. The picked reflectors were finally interpolated and a map showing the depth of each interface was created. The interpolation was carried out using the Flex Gridding algorithm (more information is available in the IGS Kingdom online manual, https://kingdom.ihs.com) with a cell size of 25 m. For converting reflection time to depth, we used velocity-depth functions determined from common-midpoint (CMP) measurements that had been carried out separately. In this way, a specific wave velocity was attributed to each layer. At the locations of the corings, we compared the waveforms of the GPR traces to the drilled stratigraphic sequences in order to evaluate the resolution of the GPR section and to identify the wave patterns indicating the main stratigraphic transitions and facies (e.g. Figure 5c). The resulting 3D model of the investigated area was finally created and visualized using the Surfer software. The grid for the 3D model was set up by applying the Nearest Neighbor Method (https://support.goldensoftware.com).
GPR velocity model
For the determination of the velocity-depth functions, the local CMP method was applied as described, for example, in Annan (2004) or Jol (2009).
The CMP soundings were performed with two GSSI GPR antennae with 200-MHz center frequency at 16 CMP sample locations. The CMP gathers were measured up to 2-m transmitter-receiver distance with a spacing of 0.1 m per record (Figure 1, green stars). The two-way traveltimes of the major reflected arrivals were picked and the normal-moveout (NMO) velocity determined from the curvature of the respective reflection hyperbolae by least-squares fitting (Figure 2). Using Dix’s formula (Yilmaz, 1988), we then converted the NMO velocity values into interval velocity values for each layer.
As the stratification of the subsurface is not strictly horizontal, the effect of dip on the NMO velocity needs to be estimated (Jacob, 2016). It turned out that the average interface dip is 8°. The corresponding error in NMO velocity is about 1%, and is thus negligible.
Acquisition of stratigraphic data by drilling
The extraction of stratigraphic information from surface geophysical measurements is a task that should be preceded by calibration through actual lithological data. For this reason, the sedimentary infilling of Duvensee was probed with a series of cores along a 50-m long transect (Figure 1d). This transect has been taken as reference profile to verify the stratigraphy close to the shoreline. A total of 12 cores was extracted with an Usinger piston corer (Mingram et al., 2007), proceeding at 1-m increments in depth. During the years of research at the Duvensee bog area, other corings had been done previously by Bokelmann and Averdieck (internal unpublished report), and their results were incorporated in our study, too.
Results
Stratigraphy of the basin from drilling
The simplified lithology of all cores extracted along the 50-m transect (Figure 1d) is shown in Figure 3. The most distinctive stratigraphic units are peat, detritus gyttja, calcareous gyttja, organic mud, clay/silty sediments and sand. The distinction between these sediment types is primarily based on tactile/visual field observations and on loss in ignition data (LOI, Heiri et al., 2001).

Stratigraphy of the reference profile. Simplified lithology of the cores presented in Figure 1d. The different sediments are indicated with different colors.
Peat deposits present a distinct dark brown/black color. Distinguishable plant remains are generally limited to modern roots or sparse decomposed fragments. The structure of this sedimentary unit depends on its degree of decomposition and pedogenization, ranging from loose blocks of compact and rather homogeneous plant matter (decomposed peat) to sub-centimetric granules (pedogenized peat). Detritus gyttja is finer and more elastic than peat. It retains a predominant organic component (LOI values between 60% and 90%; carbonates are limited to ca. 1–5%) and is generally lighter in color than peat (brown to dark brown). It can be separated in fine, coarse and elastic detritus gyttja. The term fine detritus gyttja refers to a sediment composed by a rather homogeneous fine organic matrix, with elastic and compact texture. An increase in coarser plant remains, still encased in a predominantly fine and rather elastic matrix, determines the transition to coarse detritus gyttja. The term elastic detritus gyttja is used to describe sections of fine detritus gyttja with more accentuated elastic properties and visibly lacking recognizable plant macroremains. Calcareous gyttja is composed by fine and plastic sediments with a distinctly high carbonate component (LOI values between 20% and 35%). LOI values for organic matter are generally between 15% and 20%, although bands of irregular thickness with higher organic content are present (LOI values up to 45–55%). Organic mud is a sediment characterized by a high fine organic content (brown, dark brown color), visible but sparse plant tissue remains and higher clay/silt content compared to fine detritus gyttja. This sediment has been recorded in cores that lack a calcareous gyttja unit and represents a transitional element between clayish/silty sediments and fine detritus gyttja. Accordingly, its silt/clay content appears to decrease gradually from bottom to top, together with an increasing fine organic component. Clay/silty clay sediments have a predominant fine minerogenic component (organic matter between 1% and 10%; carbonates between 6% and 12%).
Although coring operations proceeded at 1-m increments in depth, the uppermost segments of the cores are sometimes shorter than 1 m (Figure 3). This discrepancy results from sediment compression produced by the extrusion of the sediment out of the Usinger corer barrel. Notably, these compression issues appear to affect only organic layers in the first meter below the surface. Sedimentary units with a distinct minerogenic or calcareous component are only marginally affected. Understandably, sediment compression affects the measurements of stratigraphic transitions. For this reason, we focus primarily on the stratigraphy below the superficial detritus gyttja/peat layer. The relation between GPR readings and gyttja/peat transitions are discussed as well, although with a higher degree of caution.
GPR survey
After the processing of the GPR data (see Appendix S2, available online), the picking of the main horizons was performed. This procedure was carried out in feedback with the comparison of the radargrams with the provided stratigraphy from the corings. The location of some corings was not always exactly on a GPR profile; however, this information was still helpful to understand the reflection in the GPR. Figure 4 shows two examples of different profiles crossing two islands. Clearly, a reflection correlated to this feature is indicated with a round shape that becomes deeper on the sides. From the corings at the investigated area, we know that the sediment below this reflection is the basal sand. The detection of this transition was easy because the associated reflection was clearly visible in the radargrams. The shallowest interface identifiable (called Interface1 from now on) separates the top of the stratigraphy peat and coarse detritus gyttja from the underlying fine organic deposits (calcareous gyttja, fine detritus gyttja), whereas the second interface (Interface2) marks the limit between these fine organic sediments and clay. Below these two transitions we find Interface3 defining the change between clay and the basal sands. The dashed lines in Figure 4 indicate, respectively, Interface1 (yellow line), Interface2 (green line), and the red line is Interface3, marking the clay/basal sands limit. To gain more information about the reflections visible in the survey, we compared several radargrams with the stratigraphy from the drillings available in the surrounding area. The dashed lines have been set considering the stratification along the profile, although sometimes the coring was few meters distant from the profile. Considering the entire survey, we recognize five locations in which we found reflections clearly indicating the former islands.

GPR Profiles crossing two different islands. The red line indicates the reflection coming from the transition between clay and sand (Interface3), green and yellow lines indicate Interface2 and Interface1. The only location with the best stratigraphic information was in correspondence of the reference profile because of the large number of coring in that area (Figure 1d). The estimation of the depths in the surrounding area, where few corings were available, could be then affected by an error due to the interpolation between the picked horizons of the entire survey. Also where the location of one drilling was not at the same position of a GPR profile, the estimation of the main transition could be approximated.
GPR velocity model
Using the 16 CMP measurements together with the stratigraphic information from the corings (Figure 2, Table 1), we could define an average interval velocity and its standard deviation for each layer: v1 = 0.045 ± 0.005 m/ns for coarse sediments (peat and coarse detritus gyttja), v2 = 0.039 ± 0.003 m/ns for the fine organic sediments (fine detritus gyttja and calcareous detritus gyttja), v3 = 0.037 ± 0.003 m/ns for the clay layer. From the evaluation of the CMP measurement, we obtained 14 values for v1, 9 for v2, and 7 for v3 and the relative standard deviation was calculated. The velocity associated to the basal sand is v4 = 0.038 ± 0.004 m/ns coming from three values of the CMP.
Results from CMP measurements.
CMP: common-midpoint.
TWT indicates the two-way reflection time and εr is the relative permittivity value.
These interval velocities were then used to convert the travel time to interface depth.
Knowing the travel time from the picked horizons in the radargrams and the three calculated interval velocities, the depth of each interface has been computed. The radargram in Figure 5a shows the fit of this conversion to the stratigraphic column. Maximum amplitude values are visible from the top of the section until about 0.5 m depth; below this limit, the amplitude decreases. In the first 0.5 m of the radargram, a series of reflections are visible coming from different reflectors, making the interpretation of the peat on top difficult. To better understand each reflections a more detailed stratigraphy is needed. Between 0.5 and 1.10 m the amplitude decreases, indicating lower contrast and attenuation of the radar waves. The comparison with the stratigraphy shows that this part of the radargram is correlated to the fine organic sediments (fine detritus gyttja and clay). A new increase of the amplitude appears at about 1.10 m depth. It is caused by the strong dielectric parameter contrast connected with the stratigraphic change from fine organic sediments to the sand forming the core of the island.

Comparison between GPR results and the stratigraphy. (a) 1D trace from the reference profile in the same position of the Coring I.2bis already converted using three different velocities from the CMP measurements, (b) single GPR trace and λ determination, (c) velocity distribution associated with the reference profile, (d) reference profile together with the stratigraphy and the estimation of the error for the depth calculation. The transitions between the sediment layers are indicated with dashed lines and the sediment with different colors. The yellow line (Interface1) indicates the transition between coarse organic sediments (peat, coarse detritus gyttja) and more fine organic sediments (fine detritus gyttja and calcareous gyttja). The green line (Interface2) indicates the transition between fine organic sediments (fine detritus gyttja and calcareous gyttja) and clay sediments. The third reflection (Interface3) represents the transition between clay and basal sands.
In summary, we can correlate the observed amplitude changes with the transitions between the following major units: the fine organic sediments (fine detritus gyttja and calcareous gyttja) below the coarse sediments (peat and coarse detritus gyttja), followed by the basal sand. For estimating the vertical resolution of the GPR survey, the combination of stratigraphic column and depth-converted radargram can be used too (Figure 5b). Given a radiated center frequency of 200 MHz and an average wave velocity of ~0.04 m/ns, the dominant wavelength is ~0.2 m, which shifts toward ~0.3 m at deeper levels due to absorption. Using the quarter-wavelength criterion, the resolution limit can be estimated 5- to 8-cm layer thickness, which agrees with the visual appearance of the GPR section.
Taking the reference profile (Figure 1d) as a starting point, we visualize the velocity distribution in 2D (Figure 5c) by extrapolating the velocity value along the respective stratigraphic units. The stratigraphy together with the radargram and the associated standard deviation are shown in Figure 5d. The effect of the errors in the depth determination increases with depth. Accordingly, the standard deviation of the depths of the interfaces increases from 6–8 cm in the shallower parts to 25–30 cm in the deepest parts of the section.
As evident from Table 1 and Figure 5c, the velocity decreases from the upper to the lower layers by about 10%. This is caused by an interplay of the content of fine-grained clayish sediments and water saturation. The peat on top of the stratigraphy is indeed characterized by a high water content (~70% by weight), but the velocities are also attenuated by clay, the volume portion of which increases below the gyttja deposition.
Regarding the intensity of the reflections visible in the reference profile, we can recognize that the most visible reflection belongs to the clay/basal sand transition. From the theory it is known that when a propagating electromagnetic wave encounters a significant discontinuity with respect to the dielectric permittivity, some energy is reflected. The amplitude of the reflection depends on the reflection coefficient, which is related to the contrast in dielectric permittivity. Typically, the values of dielectric permittivity in wet environment are around 26 for wet sand and 46 for wet clay (Conyers, 2013). This latter value agrees qualitatively with the value of 59 (3.9 cm/ns), which we found for the clayish bottom layer. The resulting reflection coefficient at the bottom interface is 0.2 (59 from this study vs 26 from Conyers, 2013).
The permittivity values of the gyttja and peat layers resulting from the velocity values of Table 1 are 54 (4.1 cm/ns) and 49 (4.3 cm/ns), respectively. The corresponding reflection coefficients for the peat–gyttja and gyttja–clay interfaces are about 0.025 both. In comparison with the observed reflection amplitudes, this calculation obviously underestimates the reflection coefficients (cf. Figure 5), which may effectively get enhanced by intercalated thin layers and related tuning effects. Yet, it explains why the transition between peat and clayish-loamy sediments looks less pronounced as compared to the clay–sand transition. This conclusion holds as well if permittivity values from the literature are used (e.g. 37 and 46 for loam and peat, following Conyers, 2013). Between the different organic sediments like calcareous and fine detritus gyttja, the variation of the dielectric permittivity is not big enough to create an evident reflection which can clearly distinguish them. Also, more gradual than abrupt transitions can occur between fine detritus gyttja and organic mud and also between silty sand and sand. Changes in the amount and type of fluid occupying pores, changes in grain shape and grain type give also significant reflections
To verify the calculated velocity from the CMP measurements, we used the known depths of Interface2 and Interface3 from the cores and the traveltime of each interface from the GPR. The results are (for core I.2bis) v2 = 0.042 m/ns for interface2 and v3 = 0.039 m/ns for Interface3. The values are similar to the obtained from the CMP (v2 = 0.039 ± 0.003 m/ns and v3 = 0.037 ± 0.003 m/ns) and they are in the range of uncertainty.
In summary, considering these conditions, GPR is able to detect the three major transitions presented above due to a change in the characteristics of the sediments. Interface1 delimits the change between coarse and fine organic sediment, Interface2 indicates the boundary between fine organic sediment and clay and Interface3 is associated with the clay–sand transition. The step from fine detritus gyttja to calcareous gyttja is not immediately recognizable in the GPR profiles despite the difference in composition.
Development of the Duvensee landscape and settlement places as derived from the GPR subsurface model
Using the Kingdom Software, it was possible to visualize all the GPR profiles and compare them with the stratigraphy from the corings. Locally clear rounded reflections coming from the transition clay/basal sand (Interface3) were visible, defining easily the location of former islands in the ancient landscape. From the comparison of the GPR profiles and the stratigraphy, it turned out that the visible reflectors correspond to the transitions between coarse and fine sediments described in the section, ‘GPR survey’. Interface2 belongs indeed to the transition between the fine gyttja sediments and the Lateglacial clay deposited at the bottom of the ancient lake. This transition was easy to follow through the profiles despite the smaller amplitude of the radar back-reflections. The interface was weaker in the radargrams than Interface3, but possible to follow. Interface1 represents the transition from the fine gyttja and the coarse gyttja and peat sediments which filled the lake. Recognizing this transition was difficult because the radargrams somewhere did not present a clear reflection associated with this layer. The stratigraphic information was helpful to detect this reflector though the survey and the visualization at the same time of the picked horizons and stratigraphy allowed to find Interface1 in the GPR profiles also where it was not clear. Because of an error due to the interpolation, the depth estimation of this layer should be considered with caution. Nevertheless, it does not impact the general result of this study, whose primary focus is on the detection of ancient islands. The time to depth conversion of the three main interfaces was carried out using the velocity values presented in the section, ‘GPR velocity model’ with the aim to create 2D maps showing the bottom of each layer. The interpolation between the profiles was finally done allowing the creation of 2D contour maps shown in Figure 6a–c. If we focus on Interface3, we recognize at least five areas characterized by brown colors, meaning that the reflections associated with the basal sands are at shallower depths with respect to the blue areas, in which they are deeper. This means that the brown features are supposed to represent the islands.

Maps of the transition between the main sediments visible in the GPR Survey: (a) Interface1 indicates the transition between coarse sediments and fine organic sediments, (b) Interface2 indicates the transition between fine organic sediment and clay, whereas (c) Interface3 shows the transition between clay and the basal sand associated with the bottom of the lake. Furthermore, in boxes (d) and (e), the match between the brown areas (associated with the islands) and the Mesolithic camps (the yellow squares) is presented.
Now we can compare the map of Interface3 with the positions of the dwelling sites already known. Referring to Figure 6c, we notice a very good correspondence between the likely locations of the islands as obtained from the GPR survey and the archaeological evidence. The Mesolithic camps are indeed concentrated in the brown areas (islands 1 to 5) in accordance with the archaeological field studies performed so far (Groß et al., 2018).
Based on the new data, the extension of each island can be determined: island 1, for instance, has an extension of 35 m in the east-west direction and 65 m in the north-south direction, hosting the oldest excavated sites (WP 8 and WP 9). Leaving island and moving approximately 100 m to the east, we find island 2, which was ~175-m long and ~85-m wide. Here, we find the early Boreal sites (WP 1, 2, 5, 6, 11, and 12). Moving 100 m to the south, almost in the middle of the investigated area, island 3 shows up with an extension of ~42 m in EW direction and ~74 m in NS direction carrying WP 13 and WP 19. Further south, 200 m from island 3, two more islands are found: island 4, about 85-m long and ~80-m wide, and island 5, about 82-m long and ~150-m wide. These islands are hosting early Neolithic settlements.
Figure 7 shows the three major interfaces and radar velocities of the corresponding layers in a 3D visualization using Surfer. Each velocity is shown with a different color and the depth of the equivalent interface is color coded. Maximum depths are about 3 m, showing the penetration depth of the GPR in this case study. The island topography shows up clearly in the exemplary vertical sections of Figure 7b.

3D model of the interfaces and velocity distribution. (a) 3D Visualization of the main interfaces and layers derived from the 3D GPR survey, (b) vertical sections from the 3D model, profiles 1 and 2 cross the location of two supposed islands in Figure 1b.
Taking into account each layer, the volume occupied by sediment unit could be estimated.
The coarse sediments (peat and coarse detritus gyttja) on top occupy a volume of ~30.8 × 106 m3, the fine organic sediments (fine detritus gyttja and calcareous gyttja) ~23.8 × 106 m3, the clay ~37.6 × 106 m3 and the basal sand (down to 3-m depth) ~68.6 × 106 m3. These results are important for a first estimation of the water reduction at the ancient Lake Duvensee because of the sedimentation. The fish supply of the hunter-gatherers could have been indeed influenced by this factor. We used these results to calculate the area of the islands between the first late Preboreal occupation and the Neolithic for estimating the growth rate of each island, which corresponds to an average value of ~0.50 m2/yr. The hypothesis at this point is that people moved probably southward to the deeper part of the lake where water was still available. This can be the reason why the oldest Mesolithic camps are concentrated in the northern part of the surveyed area earlier covered by sediments.
For investigating the evolution of the lake environment as living space of moving hunter-gatherers and the occupation of the area with time, we have to consider the changing water table of the lake. Using the new 3D model of the lake area, we can hypothetically simulate the regression of the water between the late Preboreal and the Subboreal.
Assuming a continuous decrease of the water level from the late Preboreal to Sub-Boreal, we obtain a spatio-temporal model of the landscape evolution showing the growth of the habitable island areas and the shrinking of the lake areas over time. From the archaeological research, we know that the early Mesolithic camps were located within the overlaying organic sediments and not on the mineral soils forming the islands; therefore, the model shown in Figure 8 includes Interface3, Interface2 and a changing water level. Its maximum depth is about 3 m, corresponding to the penetration depth of GPR. The stratigraphy visible in the corings shows that the transition clay/basal sand is in the eastern part of the lake below 3 m depth. Therefore, the model has to be considered with some care at the boundaries. Figure 8 shows clearly that the oldest Mesolithic camps (9000 cal. BC) are concentrated in the northern part of the ancient lake (WP 8 and WP 9). The 3D reconstruction suggests that island 1 was indeed a lone emerged area at that time. This may be the reason why it became the place of the first colonization during the lowering of the water level in the Preboreal. During the early Boreal, island 2 and island 4 became the location for new hunter-gatherer camps and further in time also islands 3 and 5 were chosen as a place to settle, even for Neolithic groups (WP 15 and WP 17).

3D reconstruction of the Duvensee area with a hypothetical regression of the water level and the occupation of the islands by the hunter-gatherers (color coded squares). A linear water decrease starting at 9000 BC and ending at 3000 BC.
Discussion
General assessment
In this paper, we presented the results of a spatio-temporal landscape transformation based on a GPR survey and coring analysis. For the first time the results allow a reconstruction of the development of parts of the Lake Duvensee landscape during its hunter-gatherers occupation. The derived four-dimensional (4D) model (3D spatial plus time) shows the growth of the habitable area on the islands over time, the development of pathways between the islands, the shrinking of the lake area, and it provides thereby estimates of organic and non-organic sediment volumes deposited in the lake depression. This may be used as an input or boundary condition for dynamic landscape modeling or to infer on the water volume hosting the fishing resources in this area. As a by-product of this archaeological investigation, the present model can be used as a tool for a sustainable groundwater management necessary in the frame of a general heritage preservation strategy.
Apart from methodical details, we see a general confirmation of the model in the agreement between the spatio-temporal pattern of the places found, areal growth of the islands and the time when they emerged. Clearly, details of the presented spatial model of the Duvensee bog are debatable because especially the organic parts of the sediment fill show complex forms of layering that could not be fully resolved in the vertical direction and were difficult to follow over large distances. Therefore, the interpretation had to concentrate on three major interfaces that could be traced through the whole investigation area. As to the reliability of the spatial part of the model, it has to be highlighted that these major interfaces and the included facies have multiply been checked and depth calibrated by drilling and stratigraphic analysis. We see this combination of GPR surveying and drilling as the essential part of the chosen methodology.
Core compaction
Because of their central role for ground truthing, the information obtained from drilling demands a special discussion. This is because the extrusion of the soil from the Usinger core barrel can lead to a strong compaction of the core, which complicates the depth matching of the stratigraphic column and the interfaces depths detected by GPR (giving an estimation error up to 35 cm). In our case, the compaction affected only the first layer made up of peat and coarse detritus gyttja, while the minerogenic and calcareous sediment were marginally or not affected. Therefore, the compaction leads to some uncertainty regarding the depth of Interface1. The depth of the clay-to-basal sand transition (Interface3) is therefore not affected.
Although the depth of the sandy base was the major target of this study, we find it also important to estimate the compaction ratio to clarify the depth of Interface1. For this purpose, we compared the GPR depth section along the reference profile and the core catena.
Considering an average depth of the peat coming from the stratigraphy, the corresponding GPR reflector has been set and compaction has been estimated. The compaction for core I.5, for instance, is 21%, for core I.8, is 53% and for core I.4bis, it is 38% and it does not present a homogeneous behavior; therefore, each core has to be taken into consideration for this calculation.
Timing aspects
As for the timing of sedimentation, it is important to distinguish between our reconstruction and the more complex sedimentary dynamics occurring within water basins, especially in the organic sections. Our model is capable of differentiating between broad sedimentary units, yet it treats each one of them as occurring synchronously across the whole study area. In reality, the deposition of different units co-occurred in different sedimentary environments (Averdieck, 1986; Bokelmann, 2012) and only pollen analysis and radiometric dating are capable of determining the time of deposition. Despite the lack of this chronological component, we find it plausible to assume that more elevated areas emerged above the water level earlier than the more shallow ones (e.g. after a gradual lowering of the water table, or due to a gradual accumulation of organic matter), thus allowing us to model the evolution of the Duvensee landscape through time with a simple water-level model. We see the plausibility confirmed by the already mentioned agreement of the spatio-temporal patterns of the archaeological finds and the appearance and growth of island areas.
Subsidence is also a factor that should be ideally considered when modeling ancient landscapes. This phenomenon is influenced by soil/sediment porosity (Allis, 2000), and as such it might have affected the thickness of the most compressible sedimentary units within our case study (e.g. peat and gyttja deposits). Given the highly irregular topography of Duvensee and the variable thickness of its sedimentary layers, it is reasonable to assume that ground subsidence affected different areas of the basin with variable degrees of severity over time. Modeling this complex behavior is a task beyond the capabilities of the current model, and therefore – while we acknowledge its importance – it will not be addressed further in the discussion of our results.
Conclusions
To study the landscape development at the early Mesolithic sites of ancient Lake Duvensee, we derived a 3D model of the major stratigraphic units covering an area of about 63 ha.
The 3D model is based on combining an areal GPR survey with geoarchaeological information from exemplary corings. In this approach, GPR was used to detect major units in terms of dielectric permittivity contrasts and to trace the respective layer interfaces all over the area, whereas the drill cores served for identifying the units lithologically, for calibrating the layer depths determined by GPR CMP sounding and for assessing the spatial resolution of the GPR depth sections. We estimated that the depth determined via CMP measurements were accurate within ± 10% and that the effective stratigraphic resolution (corresponding to half a wavelength) of the applied 200-MHz GPR unit was of the order of 0.10 m near the surface and 0.15 m at the maximum sounding depth of about 3 m.
We observed that the major reflections in the GPR records were generated from interfaces between layers that differ significantly in grain size. Therefore, the 3D model consists of four layers separated by the following three main interfaces (from top to bottom):
Interface1 represents the transition between the coarse organic sediments (peat and coarse detritus gyttja) at the surface and the underlying fine organic sediments (e.g. fine detritus gyttja, calcareous gyttja);
Interface2 represents the transition between fine organic sediments and underlying clayish-loamy sediments in the bottom of the previous lake;
Interface3 marks the transition between clayish-loamy layer and the basal sand deposits.
Based on the 3D stratigraphic model of the Duvensee bog and a hypothetical linear time model of the height of the water level, we come to the following conclusions:
There were five former islands that occurred successively in the late Preboreal and early Boreal, which grew and eventually connected in the Sub-Boreal.
The maximum sizes of the islands before they merged was between ~180,000 and ~3500 m2. The growth rates were of the order of ~0.25 to ~ 0.78 m2/yr during this period. Vice versa, the area of the surrounding lake shrunk at the same rate, thus possibly reducing the fish supply of the hunter-gatherers.
The locations of the islands and their estimated times of emergence from the water agree with the spatio-temporal pattern of the previous archaeological finds, which therefore support our model.
The growth pattern of the islands and the land bridges, which formed between them before the whole area landed up, highlights the growing extension of the hunter-gatherer occupation. With the regression of the water and the overgrowing peat from the shore to the center of the lake, the described islands became suitable for the settlement of new camps, as confirmed by both our topographic model and the occurrence of archaeological findings.
The 3D model enabled determining the volumes of the peat, gyttja, and clayish-loamy inorganic layers, thus providing important constraints for the numerical modeling of sedimentation in the bog area, which may be matter for future investigation.
The results show that GPR is quite a cost-effective and useful way for understanding past wetland areas. As the overgrowing processes smeared topographical information, modern remote sensing techniques present a useful way for understanding past landscapes and consequently human behavior and settlement strategies within. Even though people during the early Mesolithic were not keeping up long-term settlements, they needed an understanding of their surroundings and its potential. Therefore it is mandatory to understand how the existence of a settlement was intertwined with the surrounding environment and which potentials it provided.
Besides archaeological considerations, the developed 3D stratigraphical model can be used as an input for modeling software applied for sustainable groundwater management (Beall et al., 2011), which is a necessary component of a comprehensive heritage preservation strategy.
Supplemental Material
Supplementary_materials – Supplemental material for Reconstructing the palaeoenvironment at the early Mesolithic site of Lake Duvensee: Ground-penetrating radar and geoarchaeology for 3D facies mapping
Supplemental material, Supplementary_materials for Reconstructing the palaeoenvironment at the early Mesolithic site of Lake Duvensee: Ground-penetrating radar and geoarchaeology for 3D facies mapping by Erica Corradini, Dennis Wilken, Marco Zanon, Daniel Groß, Harald Lübke, Diana Panning, Walter Dörfler, Katharina Rusch, Rebekka Mecking, Ercan Erkul, Natalie Pickartz, Ingo Feeser and Wolfgang Rabbel in The Holocene
Footnotes
Acknowledgements
We would like to thank the longtime excavator and researcher of the Duvensee bog sites, Klaus Bokelmann, for detailed information. Furthermore, we thank the two reviewers for their valuable comments on the text.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: We are grateful to the German Research Foundation (DFG) grant 2901391021 (SFB 1266) for funding the presented research.
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
