A synthetic map of the north-west European Shelf sedimentary environment for applications in marine science

Seabed sediment mapping is important for a wide range of marine policy, planning and scientific issues, and there has been considerable national and international investment around the world in the collation and synthesis of sediment datasets. However, in Europe at least, much of this effort has been directed towards seabed classification and mapping of discrete habitats. Scientific users often have to resort to reverse engineering these classifications to recover continuous variables, such as mud content and median grain size, that are required for many ecological and biophysical studies. Here we present a new set of 0.125^∘ by 0.125^∘ resolution synthetic maps of continuous properties of the north-west European sedimentary environment, extending from the Bay of Biscay to the northern limits of the North Sea and the Faroe Islands. The maps are a blend of gridded survey data, statistically modelled values based on distributions of bed shear stress due to tidal currents and waves, and bathymetric properties. Recent work has shown that statistical models can predict sediment composition in British waters and the North Sea with high accuracy, and here we extend this to the entire shelf and to the mapping of other key seabed parameters. The maps include percentage compositions of mud, sand and gravel; porosity and permeability; median grain size of the whole sediment and of the sand and the gravel fractions; carbon and nitrogen content of sediments; percentage of seabed area covered by rock; mean and maximum depth-averaged tidal velocity and wave orbital velocity at the seabed; and mean monthly natural disturbance rates. A number of applications for these maps exist, including species distribution modelling and the more accurate representation of sea-floor biogeochemistry in ecosystem models. The data products are available from https://doi.org/10.15129/1e27b806-1eae-494d-83b5-a5f4792c46fc.

1 Introduction

Knowledge of the geographic variation in the sedimentary environment of the seabed is required for a wide variety of marine planning and science tasks. Benthic species have differing sediment requirements and seabed mapping can therefore help identify ecologically distinct habitats (Robinson et al., 2011). Sediment type and wave and tidal regime are important determinants of the rate of natural disturbance of the seabed (Aldridge et al., 2015; Bricheno et al., 2015). The composition of sediments also has a large influence on the consequences of anthropogenic disturbance on the seabed, particularly those due to trawling (Diesing et al., 2013). The evolution of deltas is strongly influenced by sediment composition (Edmonds and Slingerland, 2009; Falcini and Jerolmack, 2010). Mapping the sediment composition and physical environment of the seabed is therefore an integral part of understanding and managing benthic environments.

The north-west European Shelf is one of the world's sea regions most impacted by human activities (Halpern et al., 2014). These impacts are dominated by fishing, and it has been estimated that over 99 % of human impact on the seabed is from trawling (Foden et al., 2011). Existing maps of seabed sediments for this region have almost exclusively focused on the territorial waters of individual states (e.g. the British Geological Survey's DigSBS250 product) or subregions (e.g. the North Sea; Basford et al., 1993). Currently the EU Mesh project is mapping benthic habitat classes across the north-west European Shelf (Vasquez et al., 2015). However, no existing research has mapped the continuous properties of sediments across this region. Here we map key parameters related to the sediment composition and the physical environment of the seabed in an area extending from the Bay of Biscay to the northern limits of the North Sea.

This study was motivated by the need for openly available datasets of the sedimentary environment for parameterizing shelf sea ecosystem models (e.g. Baretta et al., 1995; Blackford, 1997; Heath, 2012; Ruardij and Van Raaphorst, 1995) and for habitat mapping. Hence, we set out to map mud, sand and gravel percentage compositions and a set of parameters which are of particular relevance for ecosystem modelling and habitat mapping.

A key challenge to mapping seabed sediments across the north-west European Shelf is that sediment data are unavailable across the entire region. In areas with high-quality spatial sediment data, it is relatively easy to provide credible maps of sediment composition using statistical interpolation techniques. However, an alternative method is needed where there is poor or no data coverage. Recently, Stephens and Diesing (2015) demonstrated that the mud, sand and gravel percentages of the seabed in British territorial waters and a large part of the North Sea can be predicted using random forest models (Liaw and Wiener, 2002) which have environmental conditions at the seabed as predictors. Further, other work has shown clear relationships between the sediment composition of the sea floor and the energetic regime at the sea floor (Porter-Smith et al., 2004; Ward et al., 2015; Heath et al., 2016).

We extend this method by predicting the sediment composition of the seabed across the entire north-west European Shelf. However, our approach to mapping differs from that taken by Stephens and Diesing (2015), who only mapped predictions of sediment composition. Since these predictions will be less reliable than interpolated values in regions with good data coverage, we interpolate sediment composition where data are available and predict it where it is not, thus creating a synthetic picture of the seabed over the north-west European Shelf. Further, we expand the approach of Stephens and Diesing (2015) and map a number of other key parameters including seabed rock, median grain sizes of the whole sediment, sand and gravel fractions, and porosity and permeability; the outputs of these models are combined with time series of tidal and wave orbital velocities and a model of natural disturbance to provide a map of natural disturbance rate on the shelf.

The motivation for the choice of seabed parameters is as follows. Mud, sand and gravel percentages and rock cover are key determinants of the suitability of a habitat for benthic species (Gray, 2002; Thrush et al., 2003), and they strongly influence the median grain size of sediments, which plays a key role in determining the natural disturbance of the seabed (Aldridge et al., 2015). Similarly, the median grain size of the sand and gravel fraction play key roles in the properties of sandy and gravelly sediments. The median grain size of the mud fraction is critical for cohesion in muddy regions, but there are insufficient data for this to be mapped credibly. A complete representation of seabed biochemistry in ecosystem models requires knowledge of porosity and permeability, which are related to whole-sediment median grain size (Ruardij and Van Raaphorst, 1995; Lohse et al., 1993). The carbon and nitrogen content of seabed sediments were mapped because of their importance to benthic communities, sediment resuspension and the potential importance of sediment carbon stores in national carbon inventories (Avelar et al., 2017). Quantitative information about the physical environment on the seabed is necessary for the production of benthic habitat maps (Vasquez et al., 2015) and as a means to compare rates of natural and physical disturbance (Diesing et al., 2013).

Table 1 Data products created at a spatial resolution of 0.125^∘ by 0.125^∘.

Table 2 Summary of data sources used in sediment analysis. Datasets 1 (http://www.bgs.ac.uk/geoindex/wms.htm) and 8 (http://www.vliz.be/vmdcdata/nsbs/) were open access. Datasets 3, 4, 5, 9, 10 and 12 were available from the transnational database of North Sea sediment data (Valerius et al., 2015), which is a collation of data compiled by the EMODnet-Geology (http://www.emodnet.eu/geology), TOLES (http://www.belspo.be/belspo/brain-be/projects/TILES_en.PDF) and AufMod (http://www.kfki.de/de/projekte/aufmod) projects. Datasets 2, 6 and 7 were available from institutional contacts. Dataset 11 was downloaded from https://jetstream.gsi.ie/iwdds/index.html.

Figure 1 Region where the sedimentary environment was mapped. We defined the north-west European Shelf as the region between 17^∘ W and 9^∘ E and 44 and 63^∘ N, where bathymetry was less than 500 m. The solid black line demarcates the region where tidal velocities were taken from the Scottish Shelf Model, as described in Sect. 2.2.4.

2 Methods

2.1 Overview

Our goal was to produce synthesized maps of the sedimentary environment of the north-west European Shelf, which we define to be areas shallower than 500 m within the longitude and latitude range 17^∘ W to 9^∘ E and 44 to 63^∘ N (Fig. 1). There are minimal sediment data for deeper areas within this region and almost all of the observations are dominated by mud (George and Hill, 2008), so it is reasonable to assume that these regions are comprised largely of mud and will have negligible natural disturbance rates. Data products were created with a spatial resolution of 0.125^∘ longitude by 0.125^∘ latitude and are listed in Table 1.

Figure 2 Locations with field estimates of each seabed sediment parameter. Data sources are listed in Table 2.

Seabed sample coverage of this shelf region is highly heterogeneous with large expanses of the domain lacking accessible data. Hence, our strategy was to fill these voids in the sample coverage with statistically modelled values. The steps involved in mapping the sedimentary environment were therefore as follows.

1. Sediment data from a number of sources (Table 2) were compiled to create a composite dataset of mud, sand and gravel percentages, rock cover, carbon and nitrogen content of sediments, and median grain sizes.

2. In areas where we have data, we spatially interpolate the relevant statistic onto the study grid.

3. Using observations, we developed random forest (RF) models to predict sediment composition using wave and tidal velocities, bathymetric properties of the seabed and distance from the coast.

4. We then used RF-predicted values to infill regions of the mapping domain where the observed data density was insufficient for direct gridding.

5. Sediment porosity and permeability at each map grid point were derived from the whole-sediment median grain size using empirically based relationships assembled from literature data.

6. The natural disturbance rates of sediments at each gridded location were then calculated from wave and current bed shear stress and grain size estimates using sediment dislocation theory.

2.2 Data sources

2.2.1 Raw sediment data and processing

We compiled data on the sediment composition of the seabed from a large number of sources. Our analysis uses the following data: mud, sand and gravel percentages, rock cover and the median grain size of the whole sediment, sand fraction and gravel fraction. The data sources are summarized in Table 2 and the geographic locations where sediment data were available are shown in Fig. 2.

The British Geological Survey (BGS, 2013) provides mud, sand and gravel percentages and the median grain size of the sand fraction for locations in most of the United Kingdom's territorial waters. Data were downloaded from the BGS website using the offshore GeoIndex tool (http://www.bgs.ac.uk/geoindex/wms.htm). The raw BGS dataset included 26 259 records of sediment composition. However, to provide a consistent measure of mud, sand and gravel content we only used grab samples. This reduced the total number of records of sediment percentages and sand D₅₀ records to 20 857 and 13 289 respectively.

An extensive dataset of surface mud, sand and gravel percentages was compiled for the transnational database of North Sea sediments (Valerius et al., 2015). This provides 36 997 records of sediment composition, with data coming from historical records of the Federal Maritime and Hydrographic Agency (Germany), the Geological Survey of Denmark and Greenland, the Geological Survey of the Netherlands, Rikswaterstatt (the Netherlands), and the Royal Belgian Institute of Natural Sciences.

Records of whole-sediment median grain size were available from the North Sea Benthos Survey (NSBS) (Basford and Eleftheriou, 1988; Basford et al., 1993). NSBS data were downloaded from the http://www.vliz.be/vmdcdata/nsbs/ website. In total there were 219 records of the whole-sediment median grain size. These data are available as separate webpages for each location and we used the rvest package in R to convert the html code into columned csv format.

The Centre for Environment, Fisheries and Aquaculture Science (Cefas) provided sediment data which included the mud, sand and gravel percentages and the distribution of sediments by grain size. Data provided by Cefas covered a large part of English and Welsh waters. In total, Cefas provided 3814 records of mud, sand and gravel percentages and sediment distribution. However, to provide a consistent estimate of sediment type we restricted our analysis to sediments analysed using laser methodology and from the top 10 cm of the seabed. This resulted in a total of 1879 sediment records being used. Cefas did not provide estimates of median grain size. We therefore calculated the median grain size as follows. For the sediment record at each location, a cumulative curve of sediment weight percentage was calculated. We then calculated the median point of this curve and classified this as the median grain size. This was carried out for the entire sediment and also for the gravel fraction.

Two datasets were provided by Marine Scotland. The first included mud, sand and gravel percentages and whole-sediment median grain sizes for a large part of the North Sea. In total, this dataset had 1214 sediment records. The second dataset included estimates of the median grain size of the combined mud and sand fraction. These grain size data were not directly usable, so we filed out samples in which the percentage of gravel was small enough that the median grain size of the mud–sand fraction was close to that of the whole sediment. To do this we analysed Cefas data and established that when the whole-sediment D₅₀ is calculated with and without the gravel fraction for sediments with less than 10 % gravel, there is negligible difference between the estimates of D₅₀. We therefore used the BGS dataset to identify regions where the gravel fraction was below 10 %. This was carried out by first calculating the number of BGS observations in each 0.5^∘ N by 1^∘ W cell. We excluded cells with fewer than 10 observations. We then further excluded all cells in which more than 10 % the observations had 10 % of higher gravel content. In these regions we accepted the Marine Scotland data as a reasonable estimate of the whole-sediment D₅₀.

The Infomar project (http://www.infomar.ie) is mapping the seabed in Ireland's territorial waters. It has compiled a historical dataset of grab samples which show the surface mud, sand and gravel percentages in many locations in Irish waters. Data were downloaded in shape file format from the https://jetstream.gsi.ie/iwdds/map.jsp website. In total, there were 1392 records of surface sediment composition.

2.2.2 Rock data

Our aim was to classify locations as non-rock, rock at surface (i.e. approximately the top 10 cm of sediment) and rock in the approximately top 50 cm of sediment and to map the percentage of surface area in each rock classification. Historically, areas have only been mapped in a discrete fashion (e.g. the British Geological Survey's Digirock map; Gafeira et al., 2010), with relatively broad areas placed in one rock category or another. Further, there are no published large-scale datasets explicitly identifying whether locations have rock at or near the sea floor. We therefore created a composite dataset using historical survey logs for British territorial waters and borehole records for the territorial waters of Denmark, Germany and the Netherlands.

The British Geological Survey provides a database of downloadable historical logs of sediment sampling surveys (available from http://www.bgs.ac.uk/geoindex/wms.htm) with good spatial coverage for British territorial waters. These logs come in the form of scanned PDFs, and they provide written summaries of each sampling event. The analysis was restricted to corer records, which typically provided sufficient information to determine if there was strong evidence of rock at or near the seabed. Grab sample survey data were initially analysed; however, the use of grab sample records will underestimate rock levels as the grab can return sediment despite there being rock at or close to the surface. We therefore ignored grab samples.

Before analysing the PDFs we created the following categories for the records: (1) evidence shows there is no rock at the location; (2) written logs are consistent with rock at the surface or rock covered by a thin skin of sediment (approximately 10 cm); (3) written logs show that there is probably a significant layer of sediment covering rock; (4) ambiguous or an unreadable record. A Python script was written that will move through each PDF and allow an analyst to classify it. This process was randomized to ensure there was no spatial bias in classification error. In total there were 20 709 initial PDFs. Of these 149 could not be classified as rock or non-rock and were discarded. There were 18 871 records with no evidence of rock, 747 with evidence of rock at or near the surface and 942 records showing rock in the top 50 cm.

Borehole records provide reliable records of the rock composition of the seabed and the layers below it. German borehole data are available from the Geopotenzial Deutsche Nordsee project. The http://www.gpdn.de/ website provides visual records of borehole logs at a large number of locations in German territorial waters in the North Sea. A total of 862 records were visually inspected and we found no evidence of rock at or near the surface in any record. The Geological Survey of the Netherlands provides extensive borehole data. These were downloaded as individual text files from the https://www.dinoloket.nl/en/subsurface-data website. Each text file provided a record of the sediment type in each layer of the borehole in a consistent format. We first identified whether there was rock in the top 50 cm of any of the core records and found none. We therefore found no evidence of surface rock in Dutch waters. The Geological Survey of Denmark and Greenland provides borehole data for Danish waters in the North Sea. Data were available from the http://www.geus.dk website. Each borehole record is available as a separate webpage, and we therefore used the R package rvest to save the relevant html code and convert the depth profile of sediment type to csv format. We were then able to identify the sediment type in the top 10 and 50 cm at each location.

2.2.3 Carbon and nitrogen content

Diesing et al. (2017) showed that the carbon content of sediment could be credibly predicted based on a series of environmental predictors. Here we take a similar approach to the predictive mapping of the carbon and nitrogen content of sediments. Particulate organic carbon (POC) and total nitrogen (TN) content were downloaded from the Cefas Data Hub (https://doi.org/10.14466/CefasDataHub.32) and taken from Serpetti et al. (2012). It is clear that carbon and nitrogen levels in sediment are strongly determined by mud content (Serpetti et al., 2012), and each record of carbon and nitrogen content is associated with a field estimate of mud content. We therefore used mud as a predictor. However, the mud measurements in the Cefas dataset alternate between using laser and sieve methodology and therefore do not provide consistent and comparable estimates of mud content. The Cefas dataset contained 182 sediment samples for which the mud content was estimated using laser and sieve methodology, which showed a strong statistical relationship between each measure. We therefore converted each sieve estimate of mud content to a laser equivalent using a statistical relationship modelled using the lm function in R (laser mud = 3.157 × (sieve mud)^0.7225, p value: $<\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{16}}$, r²=0.93) (Fig. 3).

2.2.4 The physical environment

Depth-averaged tidal velocities were calculated as follows. For most of the study region tidal velocities were taken from the output of the Scottish Shelf Model, which is an implementation of the unstructured, finite-volume 3-D hydrodynamic model FVCOM. The spatial domain of this model covers approximately 80 % of our study domain (Fig. 1). A full description of the model is provided by De Dominicis et al. (2017), and here we use the same model run described therein. A 1-year climatology (for the years 1990–2014) of atmospheric forcings was used to run the model.

For the rest of the model domain we derived tidal velocities as follows. The Oregon State University Tidal Prediction Software (OTPS) is a well-known open-source barotropic tidal model based on the Oregon State University tidal inversion of TOPEX/POSEIDON altimeter data and tide gauge data (Egbert et al., 2010). This model was used to derive the relevant tidal components. The model can be obtained from http://volkov.oce.orst.edu/tides/otps.html. We obtained a regional tidal solution using the Oregon State University Tidal Inversion Software (OTIS) with a spatial resolution of 1/30^∘. The model satisfies the depth-integrated two-dimensional shallow water equations describing momentum balance as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathit{f}\times \mathit{u}+\mathit{u}\cdot \mathrm{\nabla }\mathit{u}+F+A_{\mathrm{H}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{u}=-g\mathrm{\nabla }(\mathit{\eta }-{\mathit{\eta }}_{\mathrm{EQ}}),\end{array}$$

and volume conservation

$$\begin{array}{}\text{(2)}& -{\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}=\mathrm{\nabla }\cdot (H+\mathit{\eta })u,\end{array}$$

where η is sea surface elevation, u is the horizontal velocity vector, f is the Coriolis parameter, F is the fractional damping, A_H is an eddy coefficient, which is assumed to be constant, H is bathymetry and η_EQ is the equilibrium tide allowing for the body tide, tidal self-attraction and loading.

Wave conditions were acquired from the ERA-Interim reanalysis (Dee et al., 2011). Significant wave height, mean wave period and mean wave direction were downloaded from the ECMWF website at http://www.ecmwf.int/en/research/climate-reanalysis/. The ERA-Interim reanalysis has a spatial resolution of approximately 79 km and a temporal resolution of 6 h. Orbital velocities at the seabed were calculated using the equations of Soulsby (2006), and the relevant equations are given in this paper's Appendix. Bathymetry for the wave and tidal model runs was attained from the high-resolution (30 arcsec) General Bathymetric Chart of the Oceans (GEBCO). With the exception of the Scottish Shelf Model output we used 2012 as the year for analysis of wave and tidal conditions.

To calculate the bed shear stress we used the equations of Soulsby and Clarke (2005) under combined wave and currents conditions on smooth and rough beds. This set of equations is reproduced in the Appendix to this paper. Root mean square shear stresses for waves plus currents were used. The calculation of bed shear stress requires the bathymetry, depth-averaged current speed, current direction, significant wave height, wave period and wave direction.

For the statistical modelling of sediment composition we used EMODnet bathymetry data. These have a spatial resolution of 1/8 arcmin by 1/8 arcmin and were downloaded from the EMODnet website (http://www.emodnet.eu/bathymetry/). Data processing and calculations were carried out in R using the packages dplyr (Wickham et al., 2017) and Rcpp (Eddelbuettel et al., 2011).

2.3 Spatial gridding and predictive modelling

The synthetic maps of mud, sand and gravel percentages, and rock cover were created as follows. First we identified regions where a statistical interpolation of the relevant parameter would give a reasonable estimate across that region. In other regions we used statistical models to predict the parameter. We assume that the environmental drivers of sediment composition are consistent across space.

Sampling coverage of sediment composition covered almost all of the North Sea, the United Kingdom's territorial waters and parts of Ireland's territorial waters (Fig. 2). Observations almost universally come from sampling programmes that aimed to provide consistent spatial coverage of a specific region (e.g. the North Sea Benthos Survey; Basford et al., 1993), and parameters can be interpolated in those regions. These regions were selected by creating an alphahull around each unique set of coordinates using the R package alphahull (Pateiro-l and Rodr, 2010). An alphahull is a convex envelope around the data points which will exclude areas outside the sampled regions and exclude large holes in the data coverage. Data were first interpolated onto a 1/16^∘ by 1/16^∘ grid and then means were calculated for each 0.125^∘ by 0.125^∘ cell. Parameters were spatially interpolated using bilinear spline interpolation using the interpp function from the R package akima (Akima and Gebhardt, 2016).

For areas outside the alphahulls we used random forest (Liaw and Wiener, 2002) models to predict each parameter. This class of model has been used to predict seabed sediments (Diesing et al., 2014; Huang et al., 2012; Li et al., 2011) and carbon and nitrogen content (Diesing et al., 2017). Random forest was developed by Breiman (Breiman, 2001). It is an ensemble-based modelling approach that makes no assumptions about the form of the relationships between predictor and response variables, does not require extensive parameterization, performs internal cross-validation and avoids over-fitting. Random forest takes an ensemble-based approach to regression. This is carried out by first growing a number of regression trees (Loh, 2011). Each tree is composed of a bootstrapped sample from, and of the same size as, the fitting data. Bootstrapped samples are drawn with replacement. Each split in the tree-building process only uses a subset of the predictor variables. Splitting the trees in this way reduces the dominance of individual variables and thus decorrelates the trees, making the trees less variable and more reliable (James et al., 2013). The average across all trees is then used for predictions. This ensemble averaging makes random forest robust to over-fitting (Breiman, 2001).

The observed mud, sand and gravel percentages summed to 100. However, there is no guarantee that separately predicted mud, sand and gravel percentages will sum to 100. We therefore predicted the mud, sand and gravel percentages separately and then a multiplier was applied to each prediction so that the predictions were adjusted to total 100. Random forests were created in R using the ranger package (Wright and Ziegler, 2017), which is a computationally efficient implementation of random forest for high-dimensional data. The number of trees was set to 2000, with mtry set to 3.

A similar process was carried out for median grain sizes and carbon and nitrogen content. Grain size data were available for large parts of the United Kingdom's territorial waters and some parts of the North Sea, while carbon and nitrogen content were exclusively available in parts of the United Kingdom's territorial waters (Fig. 2). First we used the alphahull approach to identify regions where we can interpolate the parameter. We then used statistical models (discussed in Sect. 2.3.1) to predict each parameter. In each case the sediment percentage maps discussed above were used as predictors in the mapping exercise. Maps and figures were produced using the R package ggplot2 (Wickham, 2016) and ternary diagrams were produced using the R package ggtern (Hamilton, 2017).

Table 3 Predictors used for statistical models for predicting sediment parameters. When mud, sand and gravel percentages and whole-sediment median grain sizes were used as predictors, raw field data were used in the creation of the statistical models, whereas the synthetic maps created in this study were used for model predictions.

2.3.1 Environmental predictors for random forest models and model validation

The environmental predictors used for the random forest models that predicted mud, sand and gravel percentage, rock cover, and carbon and nitrogen content are listed in Table 3. Predictors were chosen based on a review of evidence on the environmental influences on the seabed and the requirement that data were available at the necessary spatial resolution.

Tidal and wave energy levels at the seabed should strongly influence mud, sand and gravel percentages. Large grain sizes require more energy to dislodge from the seabed, and therefore high bed shear stress is associated with increases in average grain size and reductions in mud content (Ward et al., 2015; Heath et al., 2016). There is scarce evidence to determined if seabed composition is influenced by year-round bed shear stress or individual high-energy events. We therefore used mean and maximum annual tidal and wave orbital velocities as predictors in the models of sediment composition and carbon and nitrogen content. The supply of sediment from river discharges and coastal erosion influences seabed sediment composition and carbon and nitrogen. We therefore included distance from the coast as a model predictor. The distance from the coast was calculated as follows. Shape files of coasts were attained from the Global Self-consistent, Hierarchical, High-resolution Geography (GSHHG) Database (Wessel and Smith, 1996). Distance of each data point from the coast was then calculated using the R package geosphere (Hijmans et al., 2012).

Smoothness of the seabed will influence seabed disturbance and sediment accumulation and is likely an indicator of the existence of rocky outcrops. We therefore included measures of seabed roughness as predictors in each random forest. A number of methods exist to quantify the roughness of the seabed (Wilson et al., 2007). However, many of them are not independent of the slope of the sea floor and are arguably not purely measures of roughness. For example, the standard deviation of bathymetry would classify a steeply sloping but smooth part of a continental shelf as being very rough. We therefore used the standard deviation of slope and the standard deviation of the residual topography as predictors in the random forests. Residual topography is the difference between the bathymetry at a specific point and the mean bathymetry within a specified spatial window. The residual topography was calculated using a 25-cell moving window. First the mean bathymetry was calculated within each window. The standard deviation of residual topography (σ) was then calculated using the formula of Cavalli et al. (2008): $\mathit{\sigma }=\mathrm{1}/\mathrm{25}\sqrt{{\sum }_{i=\mathrm{1}}^{\mathrm{25}}(x_i-x_m{)}^{\mathrm{2}}}$, where x_i is the bathymetry in a specific cell in the moving window and the respective moving window mean bathymetry. Slope was calculated using the slope function from the R package SDMTools (VanDerWal et al., 2014). We then calculated the standard deviation of slope in a similar 25-cell moving window.

The above predictors were used for the mud, sand and gravel percentage and rock cover models. For the models of carbon and nitrogen content we also included chlorophyll, salinity and seabed temperature. Carbon and nitrogen content are influenced by biological activity and should thus be influenced by primary production levels and temperature at the seabed. The MetO-NWS-REAN-PHYS-bed-daily reanalysis was used for seabed temperature. These data were downloaded from the Copernicus Marine Environmental Monitoring Service website (http://marine.copernicus.eu/services-portfolio/access-to-products/). Daily seabed temperatures from 1995–2014 were interpolated onto each location and an annual climatology was calculated for each model grid point. Climatological (1997–2015) annual mean chlorophyll (mg m⁻³) data were derived from the level 4 North Atlantic chlorophyll concentration from satellite observations reprocessed data product, which is available from the Copernicus website. Proximity to river outflows likely influences levels of carbon and nitrogen, and salinity levels act as a proxy for this. We therefore calculated an annual climatological mean (1985–2014) of salinity from the MetO-NWS-REAN-PHYS-monthly-SAL reanalysis product available from the Copernicus website.

Our methodology involves predicting the sedimentary environment in geographically distinct regions. We therefore tested the ability of random forest models to do this credibly by using a cross-validation technique involving spatially disaggregated training and test datasets. Spatial disaggregation has been shown to be a reasonable method to avoid the excessive overconfidence that can possibly result from other training and testing methodologies of spatial models (Bahn and McGill, 2013; Roberts et al., 2017). The cross-validation method was as follows. We chose to use the spatial blocking method from Roberts et al. (2017). This places data into consistently sized and spatially separate blocks or bins. We chose to bin data at a resolution of 1^∘ longitude by 1^∘ latitude. We then used 100 iterations in which each bin was randomly assigned to training and test datasets. In each iteration the random forest was trained using the training dataset and this model was then used to predict the relevant parameter using the test data. We therefore evaluated the predictive ability of the model by calculating the mean value of each statistic in the test data for each 0.125^∘ by 0.125^∘ cell. The number of observations, and thus the observation reliability, in each cell varies significantly. We therefore calculate the weighted r² between predicted and observed values in each cell, with the number of observations used as the weighting value. Weighted correlation coefficients were calculated using the function corr from the R package boot (Canty, 2002). For the full predictive models over the entire European shelf we retrained the random forests using all available data.

Figure 3 The relationship between mud estimated from laser and sieve methodology for the same samples. For estimates of carbon and nitrogen content with only sieve-based estimates of mud content, we estimated what the mud percentage would be when calculated using laser methodology. The dashed red line shows this relationship (laser mud = 3.157 × (sieve mud)^0.7225, p value: $<\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{16}}$, r²=0.93).

Figure 4 (a) Assembled data on sediment porosity and median grain size (filled circles) and the fitted relationship (solid line). (b) Annual average permeability m⁻² of sediments from seven sites off the north-east coast of Scotland; data from Serpetti et al. (2016).

2.3.2 Median grain sizes

Sufficient median grain size data were available to provide a spatial interpolation of whole-sediment D₅₀ in most of the North Sea and large parts of the English Channel and Irish Sea. We therefore interpolated whole-sediment D₅₀ in these regions. This was carried out in the same way as for the distribution of sediment percentages using bilinear spline interpolation and interpolating solely within the alphahull which surrounds the relevant data points. Outside the alphahulls we predict the relevant D₅₀ using the mud, sand and gravel percentages in the synthetic maps created in this study.

In contrast to the mud, sand and gravel percentages, we chose not to predict median grain sizes using environmental variables. Predicting both the sediment percentages and median grain sizes separately is likely to result in contradictory predictions. For example, a model might predict a much higher median grain size than is possible given the predicted sediment percentages. We therefore chose to create a statistical model which predicts the median grain size using mud, sand and gravel percentages.

The median grain size of the gravel fraction has previously been shown to relate strongly to the mud to sand ratio (Aldridge et al., 2015). We therefore chose to model the whole-sediment D₅₀, the sand D₅₀ and the gravel D₅₀ in relation to mud, sand and gravel percentages. In all cases we used general additive models (GAMs) (Wood, 2006), which marginally outperformed random forests in terms of predictive ability.

The median grain size of the whole sediment varied by 4 orders of magnitude. Consequently, a GAM which uses the D₅₀ unaltered was incapable of credibly predicting the D₅₀ for the small-grained muddy sediments. We therefore used the following log transformation for the general additive model of the total sediment median grain size; log₁₀(D₅₀)∼te(mud, sand, gravel), where the interactions between mud, sand and gravel percentages are accounted for using a tensor product smooth (te). As with the sediment percentages, data were split into training and testing data. We randomly selected 70 % of the data and used it as the training data, and then used the remaining 30 % as the test data. Likewise, the final predictive model was created using all of the data.

For the sand and gravel fractions we used a GAM of the form D₅₀∼te(mud, sand, gravel), with a log-link function to ensure predictions were never negative. Finally, a small number of predictions for the sand and gravel D₅₀ were outside the grain size boundaries for gravel or sand respectively. In these cases we forced the modelled D₅₀ to be the largest or smallest possible grain size where appropriate. General additive models were created using the R package mgcv (Wood, 2001).

Wiesner et al. (1990)Lohse et al. (1993)Ruardij and Van Raaphorst (1995)Serpetti et al. (2012)

Table 4 Published literature with porosity estimates. These data were used to statistically model porosity in terms of whole-sediment median grain size.

2.3.3 Porosity and permeability

The porosity and permeability of sediments are quantitatively related to grain size distribution, with coarser-grained sediments having lower porosity and higher permeability. We evaluated the relationship between porosity and whole-sediment median grain size by compiling published data (Table 4). Porosity is conventionally expressed as the percentage volume of sediment occupied by void spaces of water. However, some data (Wiesner et al., 1990) expressed water percentage by weight. In this case we converted the water content data (by weight) to porosity assuming a solid material density of 2.65 g cm⁻³ and a fluid density of 1.025 g cm⁻³. There was a sigmoidal relationship between log-transformed porosity and log₁₀ grain size (mm). We therefore fitted a logistic relationship between them using Nelder–Mead optimization in the optim package in R (Fig. 4). This equation is shown below and the parameters are given in Table 5.

$$\begin{array}{}\text{(3)}& {\text{log}}_{\mathrm{10}}\text{porosity}=p_{\mathrm{1}}+p_{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{\frac{-({\text{log}}_{\mathrm{10}}\text{D}\mathrm{50}-p_{\mathrm{3}})}{p_{\mathrm{4}}}}}}\end{array}$$

To our knowledge the best dataset available on the relationship between whole-sediment permeability and median grain size is that of Serpetti et al. (2016). This dataset covered muddy sand, sand and mixed sediments sampled at approximately monthly intervals over 1 year at seven sites off the east coast of Scotland. Permeability and median grain size were measured on cores from the upper 5 cm and upper 10 cm of the seabed at each site. Most sediments are sampled at a depth of 10 cm and we therefore chose to only map permeability at this depth. The differences in annual average permeability (m⁻²) can be explained using a power function of median grain size (D₅₀, mm) (r²=0.999 for 10 cm cores). The equation was as follows:

$$\text{Permeability}={\mathrm{10}}^{-\mathrm{9.213}}{\text{D}}_{\mathrm{50}}^{\mathrm{4.615}}\text{(10\hspace{0.17em}cm cores)}.$$

Porosity and permeability were mapped across the study region using the above equations and the synthetic map median grain size.

Table 5 Fitted values and standard errors of the four parameters required for the function relating sediment porosity to median grain size.

We then used the porosity estimates and the maps of POC and TN to derive additional maps of the density of carbon (kg C m⁻²) and nitrogen (kg N m²) stored in the surface sediment layer across the shelf. This was derived from the carbon and nitrogen percentages of sediment and porosity values using the following equation.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{Carbon density }\left(\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }=\text{POC}\times \text{sediment depth (m)}\\ {\displaystyle }& {\displaystyle }\times \text{Dry sediment density }\left(\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\right)\times (\mathrm{1}-\text{porosity})\\ {\displaystyle }& {\displaystyle }=\text{TN}\times \mathrm{0.1}\times \mathrm{2650}\times (\mathrm{1}-\text{porosity})\\ {\displaystyle }& {\displaystyle }\text{Nitrogen density }\left(\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }=\text{TN}\times \mathrm{0.1}\times \mathrm{2650}\times (\mathrm{1}-\text{porosity})\end{array}$$

2.3.4 Natural disturbance

We modelled the extent to which the surface layers of the sediment were disturbed by waves and tides during the year. Disturbance was defined as an event which results in physical movement of the surface sediments due to the effects of bed shear stress. We then estimated the average percentage of area disturbed per month in each 0.125^∘ by 0.125^∘ cell over our model region. We assumed that sediments are mobilized when the bed shear stress exceeds a critical Shields threshold and that this threshold is given by the equation provided by Wilcock et al. (2009).

Disturbance could be heterogeneous in space and time within each of our 0.125^∘ by 0.125^∘ cells due to variations in grain size and shear stress. We accounted for this heterogeneity as follows. The bed shear stress on the seabed is determined by the wave and tide conditions and the whole-sediment D₅₀. However, for the mud, sand and gravel fraction the critical threshold is determined by the D₅₀ of the relevant fraction.

We therefore estimate natural disturbance using the following procedure for each day of the year.

1. Calculate the bed shear stress at each 15 min time interval using equations shown in this paper's Appendix and the whole-sediment D₅₀.

2. Determine the critical threshold at each time step for mud, sand and gravel using the respective D₅₀ and Eq. (A50).

3. Percentage of area disturbed = (${\text{Mud}}_{\mathit{\%}}\times {\text{Mud}}_{\mathrm{dist}})+({\text{Sand}}_{\mathit{\%}}\times {\text{Sand}}_{\mathrm{dist}})+({\text{Gravel}}_{\mathit{\%}}\times {\text{Gravel}}_{\mathrm{dist}})$, where Mud_dist, Sand_dist, Gravel_dist denote whether the Shields stress exceeded the critical threshold for mud, sand and gravel respectively.

We follow Aldridge et al. (2015) and use a 1-day time window to classify disturbance events. Monthly disturbance rates are then calculated by aggregating the areas disturbed in each day of the month. It is important to note that the modelled disturbance rate ignores the existence of rock at the surface. We are therefore only modelling the disturbance rate in regions with sediment cover.

Figure 5 The derivation of the synthetic map of sediment percentages. The interpolated map uses bilinear spline interpolation using sediment data over the region. The random forest map predicts the sediment percentages using a random forest model which relates the percentage to the bed shear stress and the distance to the coast. The synthesized map merges the two by using spatial interpolations where we have data and the random forest predictions where we do not.

3 Results

3.1 Sediment percentages and median grain sizes

Figure 5 shows the derivation of the synthetic map of sediment percentages. The interpolated map shows that mud (regions with greater than 50 % mud) is largely concentrated in the deep Norwegian Trench, an area in the north-western North Sea, part of the western Irish Sea and in patches on the Scottish west coast. Sandy sediments (greater than 50 % sand) dominate in the North Sea, except for those areas with high mud and a small region on the south-eastern English coast with high gravel levels. High gravel levels are seen exclusively in shallow coastal regions, with most of the English Channel having more than 50 % gravel.

The predictions of the random forest models reproduce the large-scale geographic patterns of sediment composition. The R² values of the predictions of mean sediment percentage in each 0.125^∘ by 0.125^∘ grid cell on the test data were 0.444, 0.412 and 0.476 for mud, sand and gravel percentages respectively. The models pick up most of the key geographic features revealed by the spatially interpolated map. The high levels of mud in the Norwegian Trench, west of the Isle of Man and the region of the northern North Sea are reproduced. Regions of the western North Sea with relatively high mud levels are also well represented. Similarly, the model predicts the existence of relatively high levels of mud south of Ireland.

The GAM of whole-sediment D₅₀ created using the training dataset performed well against the test data. R² was 0.85 on the D₅₀ values and 0.95 on the log10(D₅₀) values. This model had an R² of 0.98. Figure 6 shows the modelled relationship between percentages of mud, sand and gravel and the median grain size of the whole sediment. The GAM relating the sand D₅₀ to the mud, sand and gravel percentages, which was trained on the training dataset, had an R² 0.42 when compared with the test data. The R² for the GAM relating the gravel D₅₀ to the mud, sand and gravel percentages was 0.38.

Figure 7 shows the derivation of the synthetic maps of median grain sizes. Whole-sediment median grain size can be interpolated for most of the North Sea, English Channel, and the Irish and Celtic seas. It varies by approximately 3 orders of magnitude, with median grain sizes above 10 mm in the gravelly regions in the English Channel and other coastal regions and median grain sizes close to 0.01 mm in muddy regions such as that in the north-western North Sea. The median grain size of the sand fraction can be interpolated for most British territorial waters and is highest in regions which are predominantly gravelly. The median grain size of the gravel fraction can only be interpolated for parts of southern British territorial waters, and it is highest in regions of high gravel content.

Figure 8 shows the derived maps of porosity and permeability. Porosity is similar across most regions, with the exception of the muddy areas in the Norwegian Trench, north-western North Sea and the Irish Sea. Permeability varies by 18 orders of magnitude. It is highest in the gravelly regions in the English Channel and some coastal regions, and it is lowest in muddy regions.

Figure 6 (a) Predictions of a GAM that relates whole-sediment D₅₀ to the mud, sand and gravel percentages. (b) Relationship between total sediment median grain size and percentage of mud, sand and gravel. The relationship was derived using a general additive model which relates the D₅₀ to the mud, sand and gravel percentage.

Figure 7 Summary of the derivation of the synthetic median grain size maps. Where we have sufficient median grain size data we spatially interpolated a map of D₅₀. In other locations we used the synthetic map of mud, sand and gravel percentages and a GAM which relates the D₅₀ to the mud, sand and gravel percentages to predict the D₅₀.

Figure 8 Maps of porosity and permeability. The relationship between porosity and permeability and median grain size was estimated using published field data. We then predicted porosity and permeability using the synthetic map of median grain size.

Figure 9 Proportion of area in each rock classification. Areas were classified by whether there was rock at the surface or a surface sediment layer plus rock in the top 50 cm. Historical survey logs and borehole records were first interpolated to provide a map of rock cover where we have sufficient data. Random forests were used to predict rock cover elsewhere using wave and tidal velocities, bathymetry, measures of bathymetry variation and distance from the coast as predictors.

Figure 10 Derivation of the synthetic maps of particulate organic carbon (POC) and nitrogen (TN). Data were interpolated based on field observations in areas with good spatial coverage. In other regions, parameters were predicted using a random forest which had mud content and physical environmental variables as predictors.

Figure 11 Modelled monthly disturbance rate. A disturbance event was defined as a time when the bed shear stress exceeded the threshold required to move either the mud, sand or gravel portion of the sediment. The monthly disturbance rate was defined as the mean fraction of the total mud, sand and gravel area disturbed per day.

The synthetic maps of rock cover are shown in Fig. 9. Observed data indicate that the eastern North Sea is almost entirely free of surface rock. There are large concentrations of surface rock in the English Channel, south-west of England, the Bristol Channel and west of the Hebrides Islands on the west coast of Scotland. The predictions of the random forest model of rock provide credible large-scale reproductions of the geographic patterns of rock cover. Predictions of surface rock and rock in the top 50 cm have r² of 0.104 and 0.1991 when compared with mean values in each 0.125^∘ by 0.125^∘ grid cell. The random forest predictions in Fig. 9 reproduce the key rock areas shown by the spatially interpolated map. Regions where we rely on predictions are largely rock free, with the notable exceptions of the high-energy English Channel, north-west of France and west of the Faroe Islands.

The mapped carbon and nitrogen content of sediment are shown in Fig. 10. The random forest predictions show close agreement with observations. Across 100 iterations in which training and test data were spatially disaggregated, 70 % of data in the training data, there was a mean r² of 0.59 and 0.70 between predicted and observed POC and PON respectively. Carbon and nitrogen content are largely determined by mud content. Therefore the regions of high carbon and nitrogen content reflect those of large mud content.

3.2 Natural disturbance

Figure 11 shows modelled natural disturbance in each month. The deep Norwegian Trench is notable for lacking any disturbance year round. Disturbance is highest in the southern North Sea where sandy regions on the French, Belgian and Dutch coasts see disturbance events almost on a daily basis. There is a notable seasonal pattern in disturbance rates, with summer months seeing lower disturbance rates, which reflects the lower wind and wave regime in this time period.

4 Data availability

The data products listed in Table 1 can been be downloaded in csv, netcdf and ESRI grid format from https://doi.org/10.15129/1e27b806-1eae-494d-83b5-a5f4792c46fc (Wilson et al., 2017).

5 Discussion

The underlying goal of this study was to synthesize large-scale information about the physical environment of the seabed, both in terms of the characteristics of sediment and the wave and tidal regimes which cause disturbance. Using field estimates of the sediment composition of the seabed, we were able to map with high confidence the sediment composition of the North Sea and British territorial waters, and we were able to make credible statistical predictions of the sediment composition in other regions. The compiled datasets of sediment composition and disturbance regime are, as far as we know, the most extensive that exist over such a large spatial scale. A number of applications exist for these datasets, including habitat mapping and quantification of anthropogenic disturbance on the seabed.

Habitat mapping requires knowledge of the composition of seabed sediments (Galparsoro et al., 2012), and the maps we produced can be seen as complementary to previous work (e.g. the EU Mesh project; Vasquez et al., 2015). Existing habitat maps typically use discontinuous categories, and the continuous nature of the maps we have produced may be advantageous for some researchers.

Limitations and assumptions

A simplifying assumption of our study was that sedimentary environments are in a state of equilibrium or near equilibrium throughout the European Shelf. However, this is unlikely to be true everywhere. Ward et al. (2015) have argued that the coarser sediments found south-east of Ireland were inherited from prior stress regimes. Furthermore, the Irish Sea has linear tidal sand ridges, which are likely relics from an earlier more energetic stress regime (Uehara et al., 2006; Scourse et al., 2009). Reconstructions of historical tidal conditions on the European Shelf (e.g. Uehara et al., 2006; Neill et al., 2010, 2009) could potentially be included as model predictors in future modelling studies.

Our maps of rock area are broadly comparable with the existing hard substrate map for British territorial waters produced by the British Geological Survey (Gafeira et al., 2010). Both maps largely draw on historical British Geological Survey logs from sea-floor surveys; however, the philosophy and motivation of our study differed from that of the British Geological Survey. The British Geological Survey was motivated by mapping rocky reef areas for marine conservation planning purposes. Regions were classified as rock or non-rock, which inevitably leads to an overestimation of rock cover if analysts assume that all mapped rock regions are made up exclusively of rock. This is illustrated in the region west of the Hebrides Islands on the west coast of Scotland, where the British Geological Survey historical records show that the seabed is a complex mixture of rock-free seabed and rocky outcrops. However, the British Geological Survey substrate map classifies almost this entire area as rock. This classification was justifiable given the aim of identifying broad regions that may have rocky reefs. However, in applications such as species distribution modelling this approach is problematic. The classification of mixed habitats as rock could result in a priori ruling out a large amount of biological activity, such as fish spawning (Ellis et al., 2012), that is known to take place in these areas. In this case the continuous mapping approach taken by our study is likely more informative.

The confidence in our rock data products is significantly lower than that for mud, sand and gravel percentages. However, this was an expected result and was consistent with existing work (Diesing et al., 2015; Stephens and Diesing, 2015; Downie et al., 2016). The survey data we rely on were explicitly designed to estimate mud, sand and gravel percentages. In contrast, the rock data were based on interpretations of historical survey logs, which creates an additional level of uncertainty. Furthermore, our raw data revealed that rock cover shows large levels of heterogeneity. The low numbers of samples (1 or 2) available in most 0.125^∘ by 0.125^∘ grid cells means that our available estimates of rock cover are highly uncertain, which inevitably leads to a model with lower levels of predictability. Predictive modelling is also complex due to the array of conditions that appear to result in a rocky seabed. The English Channel and Bristol Channel are rocky due to the strong tidal energy regime, whereas the region west of the Hebrides Islands on the Scottish west coast is relatively rocky due to the existence of rocky outcrops. It is also possible that underlying geology plays a key role in determining rock levels. A previous study that took a similar predictive modelling approach in British waters used information about rock formations as predictors (Diesing et al., 2015; Downie et al., 2016); however, we were unable to find any comparable datasets that covered the entire north-west European Shelf.

We excluded the influence of rivers from predictive models because of a lack of large-scale data. However, it is likely that this is a key influence near large estuaries. This can be seen in the high-energy Bristol Channel, where there is both a high level of rock and a relatively high level of mud due to the contradictory influences of strong tidal currents and the sediment deposits from the river Severn (McLaren et al., 1993). The influence of river outflows is implicitly captured by the inclusion of distance from the coast as a predictor. For example, there is a large increase in the carbon content of sediments close to coasts, which is likely influenced by sediment deposits from rivers. We therefore cannot rule out the possibility that certain parameters were over- or under-predicted in coastal regions due to the influence of estuaries. Similarly, we did not include the potential effects of the horizontal transport of sediment by currents (Tiessen et al., 2017) or the cross-shore transport of wave-induced resuspended sediment due to the effects of gravity (Wright and Friedrichs, 2006; Falcini et al., 2012).

Previously, Aldridge et al. (2015) mapped the natural disturbance rates of the seabed in English territorial waters and a large part of the North Sea. Despite using different methodology and assumptions, our modelled disturbance rates were broadly similar for sandy and muddy regions. However, they deviated drastically for gravelly sediments, in particular in the English Channel. Our model typically predicted disturbance events to occur at least 10 times more often in gravelly sediments compared with Aldridge et al. (2015). This difference likely results from the assumption for median grain size. A key difference in assumptions between our work and Aldridge et al. (2015) is that we used the whole-sediment D₅₀ as the basis for the bed roughness term in the shear stress and disturbance rate calculations, whereas Aldridge et al. (2015) used the D₅₀ of the gravel fraction only. Where the seabed sediments are composed of mixtures of mud, sand and gravel fractions this leads to large differences in estimated disturbance rates. It is not clear which approach is more correct. The critical Shields stress calculations are parameterized from empirical studies of sorted sediments. The extent to which these calculations apply to poorly sorted sediments is uncertain. In fact, there is a lack of a theoretical and empirical basis for estimating the suspension and transport dynamics of sediments comprising mixtures of mud, sand and gravel. Sensitivity analysis (not shown) indicated that only using the gravel D₅₀ to determine disturbance in our model resulted in comparable disturbance levels to those in Aldridge et al. (2015). Further research is therefore necessary to reduce the level of uncertainty in our knowledge of the disturbance of mixed coarse sediments.

Appendix A: Bed shear stress and natural disturbance

The bed shear stress and sediment dislocation rates were calculated by combining the equations of Soulsby and Clarke (2005), Soulsby (2006) and Wilcock et al. (2009).

Bed shear stress was calculated using the equations of Soulsby and Clarke (2005), who developed equations which calculate combined bed shear stress under waves and currents. Finally, we use the equations of Wilcock et al. (2009) to estimate the critical threshold required for the bed shear stress to cause dislocation of sediment from the sea floor. Model inputs and outputs are given in Tables A1 and A2 respectively.

A1 Calculation of wave orbital velocities

We calculate the wave orbital velocities using the equations of Soulsby (2006). All waves are assumed to be irregular (spectral).

The the zero-crossing period T_z (s), the peak period T_p (s) and the mean wave crossing period T_m (s) are calculated depending on whether there is a JONSWAP spectrum:

$$\begin{array}{}\text{(A1)}& T_{\mathrm{z}}=\left\{\begin{array}{ll}{P}_x/\mathrm{1.28}& \text{ if }\text{JONSWAP}=\mathrm{1},\\ {\text{P}}_x& \text{ if }\text{JONSWAP}=\mathrm{0}.\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(A2)}& T_{\mathrm{p}}=\left\{\begin{array}{ll}{P}_x& \text{ if }\text{JONSWAP}=\mathrm{1},\\ \mathrm{1.28}\times P_x& \text{ if }\text{JONSWAP}=\mathrm{0}.\end{array}\right.\end{array}$$

We calculate the natural scaling period, T_n, as follows:

$$\begin{array}{}\text{(A3)}& T_{\mathrm{n}}={\left({\displaystyle \frac{h}{g}}\right)}^{\mathrm{0.5}}.\end{array}$$

Soulsby and Smallnan (1986) formulated equations which approximate the wave orbital velocity at the seabed, U_w ($\mathrm{1.48}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹), as follows.

$$\begin{array}{}\text{(A4)}& U_{\mathrm{w}}=\mathrm{0.25}\left({\displaystyle \frac{H_{\mathrm{s}}}{T_{\mathrm{n}}{\left(\mathrm{1}+A_{\mathrm{pw}}{t}_{\mathrm{pw}}^{\mathrm{2}}\right)}^{\mathrm{3}}}}\right),\end{array}$$

where

$$\begin{array}{}\text{(A5)}& t_{\mathrm{pw}}={\displaystyle \frac{T_{\mathrm{n}}}{T_{\mathrm{z}}}}\end{array}$$

and

$$\begin{array}{}\text{(A6)}& A_{\mathrm{pw}}={\left[\mathrm{6500}+{\left(\mathrm{0.56}+\mathrm{15.54}\times t_{\mathrm{pw}}\right)}^{\mathrm{6}}\right]}^{\mathrm{1}/\mathrm{6}}.\end{array}$$

Table A1 Calculation inputs.

Table A2 Calculation outputs.

A2 Intermediate calculations

We must then calculate a number of intermediate terms for the shear stress calculation.

We relate the kinematic viscosity υ (m² s⁻¹) to sea water density (1026.96 kg m⁻³) and kinematic velocity μ (kg (ms)⁻¹).

$$\begin{array}{}\text{(A7)}& \mathit{\upsilon }={\displaystyle \frac{\mathit{\mu }}{\mathit{\rho }}}\end{array}$$

If we define ϕ_c as the current direction and ϕ_w as the wave direction, then we can calculate ϕ_c as the angle between the current and wave direction as follows:

$$\begin{array}{}\text{(A8)}& \mathit{\varphi }=\left\{\begin{array}{ll}({\mathit{\varphi }}_{\mathrm{c}}-{\mathit{\varphi }}_{\mathrm{w}})\cdot \mathit{\pi }/\mathrm{180}& \text{ if }{\mathit{\varphi }}_{\mathrm{c}}>{\mathit{\varphi }}_{\mathrm{w}},\\ ({\mathit{\varphi }}_{\mathrm{w}}-{\mathit{\varphi }}_{\mathrm{c}})\cdot \mathit{\pi }/\mathrm{180}& \text{ if }{\mathit{\varphi }}_{\mathrm{w}}>{\mathit{\varphi }}_{\mathrm{c}},\\ \mathrm{0}& \text{ if }{\mathit{\varphi }}_{\mathrm{w}}={\mathit{\varphi }}_{\mathrm{c}}.\end{array}\right.\end{array}$$

The bed roughness length, z₀, is calculated as follows.

$$\begin{array}{}\text{(A9)}& z_{\mathrm{0}}={\displaystyle \frac{{\text{D}}_{\mathrm{50}}}{\mathrm{12}}}\end{array}$$

The Reynolds number for currents is calculated as

$$\begin{array}{}\text{(A10)}& {\mathit{\text{Re}}}_{\mathrm{c}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{U}h}{\mathit{\upsilon }}}.\end{array}$$

The current drag coefficient for smooth turbulent flow is calculated as follows:

$$\begin{array}{}\text{(A11)}& {\text{CD}}_{\mathrm{s}}=\mathrm{0.0001615}\times \exp \left(\mathrm{6}\times {\mathit{\text{Re}}}_{\mathrm{c}}^{-\mathrm{0.08}}\right).\end{array}$$

The current drag coefficient for rough flow is calculated using the equation

$$\begin{array}{}\text{(A12)}& C_{\mathrm{Dr}}={\left({\displaystyle \frac{\mathrm{0.4}}{\log (h/z_{\text{0}})-\mathrm{1}}}\right)}^{\mathrm{2}}.\end{array}$$

The wave semi-orbital excursion is calculated using the equation

$$\begin{array}{}\text{(A13)}& A={\displaystyle \frac{U_{\mathrm{w}}{T}_{\mathrm{p}}}{\mathrm{2}\mathit{\pi }}}.\end{array}$$

The Reynolds number for waves is calculated as follows:

$$\begin{array}{}\text{(A14)}& {\mathit{\text{Re}}}_{\mathrm{w}}={\displaystyle \frac{U_{\mathrm{w}}A}{\mathit{\upsilon }}}.\end{array}$$

The wave friction factor for smooth flow is calculated using the equation

$$\begin{array}{}\text{(A15)}& f_{\mathrm{ws}}=\mathrm{0.0521}\times {\mathit{\text{Re}}}_{\mathrm{w}}^{-\mathrm{0.187}}.\end{array}$$

The wave friction factor for rough flow is calculated using the equation

$$\begin{array}{}\text{(A16)}& f_{\mathrm{wr}}=\mathrm{1.39}{\left({\displaystyle \frac{A}{z_{\text{0}}}}\right)}^{-\mathrm{0.52}}.\end{array}$$

The bed shear stress depends on whether there are currents or waves only or whether there is a combination of waves and currents.

Case 1: $\stackrel{\mathrm{‾}}{U}>\mathrm{0}$ and U_w=0; current only and no waves

The current-only shear stress is calculated as follows. When ${\mathit{\text{Re}}}_{\mathrm{c}}<=\mathrm{2000}$, we calculate the current bed shear stress under laminar flow equations.

$$\begin{array}{}\text{(A17)}& {\mathit{\tau }}_{\mathrm{m}}={\mathit{\tau }}_{\max }={\displaystyle \frac{\mathrm{3}\mathit{\rho }\mathit{\upsilon }\stackrel{\mathrm{‾}}{U}}{h}}\end{array}$$

When Re_c>2000, we calculate the current bed shear stress under turbulent equations.

$$\begin{array}{}\text{(A18)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{mr}}=\mathit{\rho }{C}_{\mathrm{Dr}}{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\text{(A19)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{ms}}=\mathit{\rho }{\text{CD}}_{\mathrm{s}}{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\text{(A20)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{m}}={\mathit{\tau }}_{\max }=max\left({\text{tc}}_{\mathrm{r}},{\text{tc}}_{\mathrm{s}}\right)\end{array}$$

Case 4: $\stackrel{\mathrm{‾}}{U}=\mathrm{0}$ and U_w>0; waves only, not currents

We calculate wave-only stress as follows.

4.1 Laminar flow: ${\mathit{\text{Re}}}_{\mathrm{w}}<={\text{Rec}}_{\mathrm{wr}}$

$$\begin{array}{}\text{(A21)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{m}}=\mathrm{0}\text{(A22)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{m}}={\mathit{\tau }}_{\max }=\mathit{\rho }{\mathit{\text{Re}}}_{\mathrm{w}}^{-\mathrm{0.5}}{U}_{\mathrm{w}}^{\mathrm{2}}\end{array}$$

4.2 Turbulent flow: ${\mathit{\text{Re}}}_{\mathrm{w}}>{\mathit{\text{Re}}}_{\mathrm{w},\mathrm{cr}}$

$$\begin{array}{}\text{(A23)}& {\mathit{\tau }}_{\mathrm{wr}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{f}_{\mathrm{wr}}{U}_{\mathrm{w}}^{\mathrm{2}},\end{array}$$

with f_wr calculated using Eq. (A15).

$$\begin{array}{}\text{(A24)}& {\mathit{\tau }}_{\mathrm{ws}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{f}_{\mathrm{ws}}{U}_{\mathrm{w}}^{\mathrm{2}},\end{array}$$

with f_ws calculated using Eq. (A16).

$$\begin{array}{}\text{(A25)}& {\mathit{\tau }}_{\mathrm{m}}=\mathrm{0}\end{array}$$

$$\begin{array}{}\text{(A26)}& {\mathit{\tau }}_{\max }=max\left({\text{tw}}_{\mathrm{r}},{\text{tw}}_{\mathrm{s}}\right)\end{array}$$

Case 5: $\stackrel{\mathrm{‾}}{U}>\mathrm{0}$ and U_w>0; combined wave and current flow

First we calculate the critical current Reynolds for transition from laminar to turbulent flow.

$$\begin{array}{}\text{(A27)}& {\mathit{\text{Re}}}_{\mathrm{c},\mathrm{cr}}=\mathrm{2000}+{\left(\mathrm{5.92}\times {\mathrm{10}}^{\mathrm{5}}\times {\mathit{\text{Re}}}_{\mathrm{w}}\right)}^{\mathrm{0.35}}\end{array}$$

The critical wave Reynolds for transition from laminar to turbulent flow:

$$\begin{array}{}\text{(A28)}& {\mathit{\text{Re}}}_{\mathrm{w},\mathrm{cr}}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{5}}.\end{array}$$

5.1 ${\mathit{\text{Re}}}_{\mathrm{c}}\le {\mathit{\text{Re}}}_{\mathrm{c},\mathrm{cr}}$ and ${\mathit{\text{Re}}}_{\mathrm{w}}\le {\mathit{\text{Re}}}_{\mathrm{w},\mathrm{cr}}$

$$\begin{array}{}\text{(A29)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{m}}={\mathit{\tau }}_{\mathrm{c}}\text{(A30)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\max }={\left(({\mathit{\tau }}_{\mathrm{m}}+{\mathit{\tau }}_{\mathrm{w}}|\cos \left(\mathit{\varphi }\right)|{)}^{\mathrm{2}}+({\mathit{\tau }}_{\mathrm{w}}\left|\sin \right(\mathit{\varphi }\left)\right|{)}^{\mathrm{2}}\right)}^{\mathrm{0.5}}\end{array}$$

5.2 Turbulent flow: Re_c>Re_w or ${\mathit{\text{Re}}}_{\mathrm{w}}>{\mathit{\text{Re}}}_{\mathrm{w},\mathrm{cr}}$

We must recalculate τ_c and τ_w since it is possible that one or the other of these could appear as laminar flow when estimated independently. The following are for the current-only component of stress under combined wave and current turbulent conditions.

$$\begin{array}{}\text{(A31)}& {\displaystyle }& {\displaystyle }{\text{tc}}_{\mathrm{r}}=\mathit{\rho }{C}_{\mathrm{Dr}}{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\text{(A32)}& {\displaystyle }& {\displaystyle }{\text{tc}}_{\mathrm{s}}=\mathit{\rho }{\text{CD}}_{\mathrm{s}}{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\text{(A33)}& {\displaystyle }& {\displaystyle }\text{tc}=max({\text{tc}}_{\mathrm{r}},{\text{tc}}_{\mathrm{s}})\end{array}$$

The following are for the wave-only component of stress under combined wave and current turbulent conditions.

$$\begin{array}{}\text{(A34)}& {\displaystyle }& {\displaystyle }{\text{tw}}_{\mathrm{r}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{f}_{\mathrm{wr}}{U}_{\mathrm{w}}^{\mathrm{2}}\text{(A35)}& {\displaystyle }& {\displaystyle }{\text{tw}}_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{f}_{\mathrm{ws}}{U}_{\mathrm{w}}^{\mathrm{2}}\text{(A36)}& {\displaystyle }& {\displaystyle }\text{tw}=max({\text{tw}}_{\mathrm{r}},{\text{tw}}_{\mathrm{s}})\end{array}$$

5.2.1 Rough turbulent flow: ${\mathit{\tau }}_{\max ,\mathrm{r}}>{\mathit{\tau }}_{\max ,\mathrm{s}}$

We calculate tm_r and tmax_r as follows.

$$\begin{array}{ll}\text{(A37)}& {\displaystyle }& {\displaystyle }{a}_{\mathrm{r}}=\mathrm{0.24}\text{(A38)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{1}}=max\left(a_{\mathrm{r}}{\left({\displaystyle \frac{f_{\mathrm{wr}}}{\mathrm{2}}}\right)}^{\mathrm{0.5}}\left({\displaystyle \frac{A}{z_{\mathrm{0}}}}\right),\mathrm{12}\right)\text{(A39)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{2}}={\displaystyle \frac{h}{T_{\mathrm{1}}{z}_{\mathrm{0}}}}\text{(A40)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{3}}={\left(C_{\mathrm{Dr}}^{\mathrm{2}}+{\left({\displaystyle \frac{f_{\mathrm{wr}}}{\mathrm{2}}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{U_{\mathrm{w}}}{\stackrel{\mathrm{‾}}{U}}}\right)}^{\mathrm{4}}\right)}^{\mathrm{0.25}}\text{(A41)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{1}}=T_{\mathrm{3}}{\displaystyle \frac{\log \left(T_{\mathrm{2}}\right)-\mathrm{1}}{\mathrm{2}\log (T_{\mathrm{1}}}}\text{(A42)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{2}}=\mathrm{0.4}{\displaystyle \frac{T_{\mathrm{3}}}{\log \left(T_{\mathrm{1}}\right)}}\text{(A43)}& {\displaystyle }& {\displaystyle }{\text{CD}}_{\mathrm{m}}={\left({\left(A_{\mathrm{1}}^{\mathrm{2}}+A_{\mathrm{1}}\right)}^{\mathrm{0.5}}-A_{\mathrm{1}}\right)}^{\mathrm{2}}{\displaystyle }& {\displaystyle }{\text{CD}}_{\max }=\left({\left({\text{CD}}_{\mathrm{m}}+T_{\mathrm{3}}{\displaystyle \frac{U_{\mathrm{w}}}{\stackrel{\mathrm{‾}}{U}}}{\left({\displaystyle \frac{f_{\mathrm{wr}}}{\mathrm{2}}}\right)}^{\mathrm{0.5}}\left|\cos \right(\mathit{\varphi }\left)\right|\right)}^{\mathrm{2}}\right.\\ \text{(A44)}& {\displaystyle }& {\displaystyle }{\left.+{\left(T_{\mathrm{3}}{\displaystyle \frac{U_{\mathrm{w}}}{\stackrel{\mathrm{‾}}{U}}}{\left({\displaystyle \frac{f_{\mathrm{wr}}}{\mathrm{2}}}\right)}^{\mathrm{0.5}}\left|\sin \right(\mathit{\varphi }\left)\right|\right)}^{\mathrm{2}}\right)}^{\mathrm{0.5}}\text{(A45)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{r}}=\mathit{\rho }{\text{CD}}_{\mathrm{m}}{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\text{(A46)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\max }=\mathit{\rho }{\text{CD}}_{\max }{\stackrel{\mathrm{‾}}{U}}^{\mathrm{2}}\end{array}$$

Shields number calculation

We calculate the bed shear velocity as

$$\begin{array}{}\text{(A47)}& u_{\mathrm{star}}={\left({\displaystyle \frac{{\mathit{\tau }}_{\max }}{\mathit{\rho }}}\right)}^{\mathrm{0.5}}.\end{array}$$

The Reynolds particle number is calculated as follows:

$$\begin{array}{}\text{(A48)}& \mathit{\text{Re}}=u_{\mathrm{star}}\times {\displaystyle \frac{{\text{D}}_{\mathrm{50}}}{\mathit{\upsilon }}}.\end{array}$$

The Shields stress is calculated as

$$\begin{array}{}\text{(A49)}& {\text{Tc}}_{\mathrm{S}}={\displaystyle \frac{{\mathit{\tau }}_{\max }}{({\text{se}}_{\mathrm{d}}-\mathit{\rho })\times g\times {\text{D}}_{\mathrm{50}}}}.\end{array}$$

The critical Shields stress is calculated as follows from the empirical relationship shown in Wilcock et al. (2009).

$$\begin{array}{}\text{(A50)}& {\text{TcS}}_{\mathrm{cr}}=(\mathrm{0.105}\times {\mathit{\text{Re}}}^{-\mathrm{0.3}}))+\mathrm{0.045}\times \exp (-\mathrm{35}\times {\mathit{\text{Re}}}^{-\mathrm{0.59}}\left)\right)\end{array}$$

Here there is a U-shaped relationship between grain size and the critical Shields stress because of the high stresses required to dislodge bigger grains and the cohesive nature of mud.

RMS shear stress for waves and currents

The root mean square shear stress is calculated as follows:

$$\begin{array}{}\text{(A51)}& {\mathit{\tau }}_{\mathrm{rms}}={\left({\mathit{\tau }}_{\mathrm{m}}^{\mathrm{2}}+{\mathit{\tau }}_{\mathrm{w}}^{\mathrm{2}}\right)}^{\mathrm{0.5}}.\end{array}$$

Information about the Supplement

High-resolution versions of the paper's figures have been made available as a supplementary file so that figures can be used in presentations and reports.