shoredate: Shoreline dating Stone Age sites on the Norwegian Skagerrak coast

Introduction

The concept of dating sites based on their present-day elevation with reference to past shoreline displacement has been an important tool for archaeologists in northern Scandinavia since the early 1900s (e.g. Brøgger 1905). This is based on the observation that Stone Age sites in the region tend to have been located close to the contemporaneous shoreline when they were in use. Furthermore, following the retreat of the Fennoscandian Ice Sheet, the isostatic uplift has been so severe that despite corresponding eustatic sea-level rise, the relative sea-level has been falling throughout prehistory in large parts of this region. Within any given area where this is the case, the general pattern is thus that older sites will be located at higher altitudes than younger sites.

This vignette describes the R package shoredate which provides tools for shoreline dating Stone Age sites located on the Norwegian Skagerrak coast using the approach presented in Roalkvam (2023). This is based on an empirical evaluation of the likely elevation of the sites above sea-level when they were in use, and the local trajectory of past relative sea-level change. Due to the geographical contingency of the method, and the dependency on geological reconstructions of shoreline displacement, the functionality of the package was developed for the purposes of dating coastal sites located in the region stretching from Horten county in the north east, to Arendal county in the south west. However, in the second vignette, ways in which the core functionality of the package could be applied in other areas of the world are presented.

Geographical and temporal coverage

Shoreline dating of Stone Age sites in northern Scandinavia is based on the premise that coastal sites within the region tend to have been located on or close to the shoreline when they were in use. By reconstructing the past trajectory of shoreline displacement, this thus allows us to ascribe an approximate date to when the sites were in use, based on their present-day elevation. The method is therefore dependent on local relative sea-level change and the likely elevation of the sites above sea-level when they were in use.

As the method is dependent on good geological reconstructions of shoreline displacement, the full functionality of shoredate is at present limited to being applicable in the region of south-eastern Norway between Horten in the north east to Arendal in the south west. This region has recently compiled shoreline displacement curves for Skoppum in Horten (Romundset 2021), Gunnarsrød in Porsgrunn (Sørensen et al. 2014a; Sørensen et al. in prep), Hanto in Tvedestrand (Romundset 2018; Romundset et al. 2018) and Bjørnebu in Arendal (Romundset 2018). This region also formed the study area for Roalkvam (2023), in which the specification of the method and the parameters used here were derived. The core spatial coverage of shoredate is indicated in the maps below.

plot of chunk limit

Furthermore, human occupation in Norway only occurred some time after the retreat of the Fennoscandian Ice Sheet, and the oldest known sites are dated to around 9300 BCE (e.g. Glørstad 2016). 9469 BCE is the latest start date among the geologically derived displacement curves in the region, and thus marks the oldest possible age to achieve with shoredate here, although no sites are yet known to be that old. If a site has an elevation that implies a date older than the latest start date of the displacement curves, the date is returned as NA and a warning is given.

In the Roalkvam (2023) study it was found that sites tend to be situated at more removed and variable distances from the shoreline from around 2500 BCE. By default, this therefore marks the upper limit of the method, where dates with a start date after 2500 BCE are returned as NA.

Installing and loading shoredate

The package can be installed from CRAN:

install.packages("shoredate")

The development version of the package can be installed from the GitHub repository. This requires the devtools package. Note that the development version can be unstable.

# Uncomment and run to install devtools
# install.packages("devtools")

# Installation requires devtools
devtools::install_github("isakro/shoredate", build_vignettes = TRUE)

When installed the package can be loaded.

library(shoredate)

Interpolating shoreline displacement to a site location

The function shoreline_date() forms the basis for shoredate, and is presented in more detail below. Unless a pre-compiled displacement curve is provided when calling shoreline_date(), the trajectory of shoreline displacement at the location of the sites to be dated is interpolated under the hood. This is based on the four different geological reconstructions of shoreline displacement along the Norwegian Skagerrak coat. Each of these displacement curves are associated with a shoreline isobase which indicates contours along which the relative sea-level change has followed the same trajectory. The variation in the rate of sea-level change is mainly determined by the degree of isostatic uplift, which is higher along a gradient than runs from the south west to the north east within the region. The shoreline displacement curve for any given location within the region is interpolated using the interpolate_curve() command, using inverse distance weighting, where the distance to the four isobases is weighted by the square of the inverse distances.

Note that spatial geometries representing sites to which interpolate_curve() is to be applied must be set to the coordinate reference system (CRS) WGS 84 / UTM zone 32N (EPSG:32632). Furthermore, a warning is given if one attempts to interpolate the trajectory for shoreline displacement at a site located outside the spatial extent outlined above. However, as there might still be use-cases where it could be useful to extrapolate the method outside this limit, the procedure is still performed (see the second vignette Applying shoredate to other regions for suggestions on how to apply the package to entirely different regions).

# Create example point using the required coordinate system
target_point <- sf::st_sfc(sf::st_point(c(517250, 6527250)), crs = 32632)

# Interpolate the displacement curve
target_curve <- interpolate_curve(target_point)

To visualise the interpolated curve alongside the geological reconstructions, this can be passed to displacement_plot(). A simple map showing the location of one or more target sites relative to the isobases can also be displayed by using the command target_plot().

target_plot(target_point)

plot of chunk plot_location

# The opactity of the geological displacement curves is reduced with
# displacement_alpha
displacement_plot(target_curve, displacement_alpha = 0.4)

plot of chunk plot_intcurve

Shoreline dating

The following section outlines how to perform shoreline dating and ways in which to visualise the results. Various parameters and helper-functions underlying the dating procedure are also presented, along with ways to manipulate these if one wishes to explore the sensitivity of the dates or believe other parameters than the defaults are more sensible in a given setting.

Dating a single site

To perform a basic shoreline date, a site has to be provided as a spatial sf object if the shoreline displacement is to be interpolated to its location. In a addition to this, the elevation of the site above the present sea-level has to be provided. This can either be done by manually specifying the elevation, or by providing an elevation raster (see below), in which case the elevation of the site will be derived from the raster, provided the site is given as a spatial object. In Roalkvam (2023) it was found that the vertical distance from the sites to the contemporaneous shoreline within the study region could be reasonably approximated by a gamma distribution. This forms the foundation for the implementation of shoreline dating in the package.

# Date the example point, manually specifying it's elevation
target_date <- shoreline_date(target_point, elevation = 58.5)

The results can then be plotted using the shoredate_plot() command:

shoredate_plot(target_date)

plot of chunk plot_example

The blue gamma distribution on the y-axis indicates the likely elevation of the site above sea-level when it was in use, which is described by an empirically derived gamma distribution with the parameters shape (\(\alpha\)) = 0.286 and scale (\(\sigma\)) = 20.833. These parameters can be adjusted by specifying the values passed to model_parameters when calling shoreline_date(). The red envelope is the shoreline displacement curve as interpolated to the site location with interpolate_curve(), which is run under the hood if a curve is not provided to the target_curve argument of shoreline_date().

Transferring the probability from the gamma distribution to the calendar scale using the displacement curve gives the resulting shoreline date in grey, which is underlined with the 95% highest density region (HDR) in black. The coverage of the HDR can be adjusted with the hdr_prob parameter. The dating procedure involves stepping through the cumulative version of the gamma distribution, and at each step subtracting the probability at the previous step from the current one. The probability at each step is divided uniformly across the calendar years in the range between the lower and upper limit of the displacement curve. The default resolution on the calendar scale is 10 years, but this can be adjusted by specifying the cal_reso argument. The gamma distribution is stepped through using increments of 0.01m, which can be adjusted by using the argument elev_reso. By default, the resulting shoreline date is normalised to sum to unity.

The function shoreline_date() returns an object of class shoreline_date. By calling the object the dates are printed to console, displaying site name, site elevation and the HDR of the date. Furthermore, the first column of a data frame beyond the geom column holding the sf geometry will be used as site names. If no such column is present, the sites will simply be numbered as they are passed to shoreline_date().

target_date

## ===============
## Site:  1
## Elevation:  58.5 
## 
## 95% HDR:
## 8600 BCE-6900 BCE
## 5730 BCE-5070 BCE

While a range of graphical parameters can be adjusted to change the appearance of the output plot, setting greyscale = TRUE offers a short-cut to plot a date in greyscale:

shoredate_plot(target_date,  greyscale = TRUE)

plot of chunk greyscale

A more sparse version can also be plotted by specifying what elements are to be excluded from the plot:

shoredate_plot(target_date, elevation_distribution = FALSE,
               displacement_curve = FALSE, highest_density_region = FALSE)

plot of chunk sparse

Finding the terminus post quem date

It is also possible to not account for the likely distance between the site and the shoreline. This is done by setting the model parameter to “none” when calling shoreline_date(). This thus effectively provides a termnius post quem date, the earliest possible date for when the site could have been in use, under the assumption that when the sea receded from the site location marks the earliest possible point in time for the occupation of the site.

# Date the example point, manually specifying it's elevation and setting the
# HDR coverage to 1.
tpq_date <- shoreline_date(target_point,
                           elevation = 58.5,
                           model = "none",
                           hdr_prob = 1)

# Plot and adjust the position of the HDR label on the y-axis
shoredate_plot(tpq_date, hdr_label_yadj = 0.6)

plot of chunk tpq_date

Dating multiple sites

As exemplified in the code chunks to follow, it is also possible to date multiple sites at once.

# Create multiple example points to be dated
target_points <- sf::st_sfc(sf::st_point(c(576052, 6567955)),
                              sf::st_point(c(554212, 6538835)),
                              sf::st_point(c(516599, 6512142)))

# Specify the correct CRS and make the points a sf data frame to be able to
# add the column below
target_points <- sf::st_as_sf(target_points, crs = 32632)

# Adding example site names
target_points$names <- c("Site 1", "Site 2", "Site 3")

# Create a plot showing the location of the points within the spatial limit
target_plot(target_points)

plot of chunk multipoint

When dating multiple sites it can be useful to track the progress by setting verbose = TRUE to print the progress to console (this output is not reproduced here).

# Dating the target points, specifying the elevations in a
# vector of numeric values
target_dates <- shoreline_date(target_points,
                               elevation = c(68, 98, 73),
                               verbose = TRUE)

# Print results
target_dates

## ===============
## Site:  Site 1
## Elevation:  68 
## 
## 95% HDR:
## 7630 BCE-6110 BCE
## 6000 BCE-5430 BCE
## ===============
## Site:  Site 2
## Elevation:  98 
## 
## 95% HDR:
## 9020 BCE-8470 BCE
## ===============
## Site:  Site 3
## Elevation:  73 
## 
## 95% HDR:
## 9170 BCE-8340 BCE

The default behaviour when passing multiple dates to shoredate_plot() is to plot a series of individual plots. However, setting multiplot = TRUE collapses the dates on a single, more sparse plot. The sites are ordered from from earliest to latest possible start date for the occupation of the sites.

shoredate_plot(target_dates, multiplot = TRUE)

plot of chunk multiplot

Using an elevation raster to find site elevations

The below code demonstrates how to find the elevation of one or more sites by passing a digital elevation model (DEM) to shoreline_date(). The example here uses data from Amazon Web Service Terrain Tiles, retrieved with the package elevatr. However, the best elevation data available for the study area can be retrieved freely from the Norwegian Mapping Authority at https://hoydedata.no/LaserInnsyn2/. To use this with shoredate, manually download a raster covering the location of the sites to be dated and read this in to R with terra::rast() before passing it to shoreline_date().

# Reproject target_point for retrieving raster data with elevatr
target_wgs84 <- sf::st_transform(target_point, crs = 4326)

# Retrieve raster data
elev_raster <- elevatr::get_elev_raster(target_wgs84, z = 14,  src = "aws")

## Warning in matrix(unlist(x, use.names = FALSE), nrow = length(x), byrow = TRUE, : data length [3] is not a sub-multiple or multiple of the number of
## rows [2]

# Make the retrieved raster a SpatRaster and re-project
# it to the correct coordinate system
elev_raster <- terra::project(terra::rast(elev_raster), "epsg:32632")

# Plot the raster and target point for inspection
terra::plot(elev_raster)
plot(target_point, add = TRUE)

plot of chunk retrieve_raster

When a SpatRaster has been loaded into the R session this can be passed to shoreline_date() by specifying it in the elevation argument. This is then used to find the site elevations without having to provide these manually, as was done above. Note that if a site is represented by an object with a spatial extent (i.e. not a point feature) the default is to use the mean of the elevation on the raster. This can be changed with the elev_fun argument which accepts any function that can be passed to the fun argument in terra::extract() (see ?terra::extract).

raster_date <- shoreline_date(target_point, elevation = elev_raster)

# Print to console
raster_date

## ===============
## Site:  1
## Elevation:  58.47393 
## 
## 95% HDR:
## 8600 BCE-6890 BCE
## 5730 BCE-5070 BCE

Summed probability distribution of shoreline dates

Using sum_shoredates() it is also possible to sum the probability of multiple shoreline dates to evaluate the temporal intensity of the dates, as is frequently done with calibrated radiocarbon dates (e.g. Crema and Bevan 2021). Do note that as this procedure involves collapsing the frequency of dates and their uncertainty, the interpretation of the resulting sum is not straightforward, as it’s shape is determined not only by the number of sites but also the trajectory of shoreline displacement (see e.g. Bronk Ramsey 2017; Crema 2022 for perspectives and ways to approach this issue in the context of radiocarbon dates).

# Sum the three dates from above
summed_dates <- sum_shoredates(target_dates)

A simple line plot displaying the resulting summed probability can then be produced with shoredate_sumplot().

shoredate_sumplot(summed_dates)

plot of chunk sum_plot

Direction of shoreline gradient

The direction of the shoreline gradient is specified by the isobases that run perpendicular to this gradient. When interpolating the displacement curve used for dating a site, the direction of these isobases defaults to 327\(^\circ\), following Romundset et al. (2018:180, fig.1). However, some authors operate with different isobase directions, with Sørensen et al. (2014:fig.2.2.3) using 338\(^\circ\). Furthermore, the direction of the uplift gradient might have varied throughout the Holocene (Sørensen et al. 2014:42–44).

Specifying different isobase directions

While the range of the typically employed isobase directions generally result in minor differences in the resulting shoreline dates, it is possible to specify other and multiple directions for these, and perform the shoreline dating using each individual direction. This can be useful either for evaluating the sensitivity of a date or if one believes there is reason suspect a different direction than the default in a particular case.

# Add a name to the example site
target_point <- sf::st_as_sf(target_point, crs = 32632)
target_point$name <- "Example site"

# Using the same target point and elevation as above, but specifying two
# different directions for the isobases when dating
iso_date <- shoreline_date(site = target_point, elevation = elev_raster,
                           isobase_direction = c(325, 338))

iso_date

## ===============
## Site:  Example site
## Elevation:  58.47393 
## 
## Isobase direction:  325 
## 
## 95% HDR:
## 8610 BCE-6900 BCE
## 5710 BCE-5070 BCE
## 
## Isobase direction:  338 
## 
## 95% HDR:
## 8560 BCE-6900 BCE
## 5810 BCE-5070 BCE

In the call to plot it can then be specified that the direction of the isobases is printed with each date. Note that it is not possible use multiplot when having used multiple isobase directions when dating the sites. The following thus also shows the default behaviour when plotting multiple dates at once and not setting multiplot = TRUE.

shoredate_plot(iso_date,
               site_name = TRUE,
               isobase_direction = TRUE,
               hdr_label_yadj = 0.35)

plot of chunk plot_example_isobases

Sum the probability of dates across isobase directions

As an alternative to keeping the two shoreline dates of the site separate, the parameter sum_isobase_directions can be set to TRUE when calling shoreline_date() to sum the probability of the dates. Again, by default this is normalised to sum to unity.

sum_iso_date <- shoreline_date(site = target_point,
                               elevation = elev_raster,
                               isobase_direction = c(325, 338),
                               sum_isobase_directions = TRUE)

sum_iso_date

## ===============
## Site:  Example site
## Elevation:  58.47393 
## 
## Sum of isobase directions:  325 338 
## 
## 95% HDR:
## 8610 BCE-6880 BCE
## 5770 BCE-5060 BCE

This can then be passed to shoredate_plot(). Note that the interpolated displacement curves and the gamma distribution for the elevation of the site above sea-level is not included in the plot when summing the probability of dates across isobase directions.

shoredate_plot(sum_iso_date, site_name = TRUE, isobase_direction = TRUE)

plot of chunk plot_sum_isobases

References

Bronk Ramsey, C. 2017. Methods for Summarizing Radiocarbon Datasets. Radiocarbon 59(6):1809–1833. DOI: 10.1017/RDC.2017.108

Brøgger, W.C. 1905 Strandliniens beliggenhed under stenalderen i det sydøstlige Norge. Geological Survey of Norway, Kristiania. URI: https://hdl.handle.net/11250/2675184

Crema, E. and Bevan, A. 2021. Inference from Large Sets of Radiocarbon Dates: Software and Methods. Radiocarbon 61(1):23–39. DOI: 10.1017/RDC.2020.95

Crema, E. 2022. Statistical Inference of Prehistoric Demography from Frequency Distributions of Radiocarbon Dates: A Review and a Guide for the Perplexed. Journal of Archaeological Method and Theory 29:1387-1418. DOI: 10.1007/s10816-022-09559-5

Glørstad, H. 2016. Deglaciation, sea-level change and the Holocene colonization of Norway. Geological Society, London, Special Publications 411(1):9–25. DOI: 10.1144/SP411.7

Roalkvam, I. 2023 A simulation-based assessment of the relation between Stone Age sites and relative sea-level change along the Norwegian Skagerrak coast. Quaternary Science Reviews 299:107880. DOI: 10.1016/j.quascirev.2022.107880

Romundset, A. 2018. Postglacial shoreline displacement in the Tvedestrand–Arendal area. In The Stone Age Coastal Settlement in Aust-Agder, Southeast Norway, edited by Reitan, G. and Sundström, L.. Cappelen Damm Akademisk, Oslo, pp. 463–478. DOI: 10.23865/noasp.50

Romundset, A. 2021. Resultater fra NGUs undersøkelse av etteristidas strandforskyvning nord i Vestfold. Geological Survey of Norway, Trondheim.

Romundset, A., Lakeman, T.R. and Høgaas, F. 2018. Quantifying variable rates of postglacial relative sea level fall from a cluster of 24 isolation basins in southern Norway. Quaternary Science Reviews 197:175e192. DOI: 10.1016/ j.quascirev.2018.07.041

Sørensen, R, Henningsmoen, K.E. Høeg, H.I. and Gälman V. 2023. Holocen vegetasjonshistorie og landhevning i søndre Vestfold og sørøstre Telemark. In The Stone Age in Telemark. Archaeological Results and Scientific Analysis from Vestfoldbaneprosjektet and E18 Rugtvedt–Dørdal, edited by Persson, P. and Solheim, S., in press.

Sørensen, R, Henningsmoen, K.E. Høeg, H.I. and Gälman V. 2014. Holocene landhevningsstudier i søndre Vestfold og sørøstre Telemark – Revidert kurve. In Vestfoldbaneprosjektet. Arkeologiske undersøkelser i forbindelse med ny jernbane mellom Larvik og Porsgrunn. Bind 1, edited by Melvold, S. and Persson, P. Portal forlag, Kristiansand, pp. 36–47. DOI: 10.23865/noasp.61

Isak Roalkvam