Using Open-Data Portals, Remote Sensing and Computational Modelling to Investigate Historic Wreck Sites and Their Environments: 45 Years on from Muckelroy

ABSTRACT

We explored site formation processes at 549 shipwrecks in the North and Baltic Seas using multibeam echosounder (MBES) data and environmental/anthropogenic variables. Phase 1 classified the structural integrity of each wreck. Phase 2 used regression modelling to explore correlations between wreck class and environmental variables at individual sites. A best-fit model comprised depth, salinity, fishing intensity, and wreck year. In phase 3 multi-criteria analysis predicted preservation potential across the full model domain. Highest preservation potential is modelled for the southern Baltic Sea and Norwegian Trench, with low potential in shallower areas of the North Sea and off England, Netherlands and Germany.

RESUMEN

Exploramos los procesos de formación de sitio de 549 pecios de los mares del Norte y Báltico, empleando datos de ecosonda multihaz y variables ambientales/antropogénicas. En la Fase 1 se clasificó la integridad estructural de cada pecio. En la fase 2 se utilizaron modelos de regresión para explorar las correlaciones entre tipo de pecio y variables ambientales en sitios individuales. El modelo que mejor se ajusta incluye profundidad, salinidad, intensidad de pesca y el año de naufragio. En la fase 3, un análisis multicriterio predijo el potencial de preservación a lo largo de todo el dominio del modelo. Se calculó el mayor potencial de preservación para el sur del Mar Báltico y la Fosa Noruega, y un bajo potencial en áreas más superficiales del Mar Norte y frente a Inglaterra, Países Bajos y Alemania.

摘要

我们利用多波束測深探測儀（MBES）数据和环境/人为变量探查了北海和波罗的海549艘沉船的遗址形成过程。阶段1是对每艘沉船的结构完整性进行分类。阶段2使用回归模型探索单个遗址中沉船类别与环境变量之间的相互关联。最佳拟合模型包括深度、盐度、捕捞强度和沉船年代。在阶段3以多标准分析预测全模型领域的保存潜力。波罗的海南部和挪威海沟的保存潜力最高，北海较浅区域以及英格兰、荷兰和德国沿海的保存潜力则较低。

摘要

我們利用多波束測深探測儀（MBES）數據和環境/人為變量探查了北海和波羅的海549艘沉船的遺址形成過程。階段1是對每艘沉船的結構完整性進行分類。階段2使用回歸模型探索單個遺址中沉船類別與環境變量之間的相互關聯。最佳擬合模型包括深度、鹽度、捕撈強度和沉船年代。在階段3以多標準分析預測全模型領域的保存潛力。波羅的海南部和挪威海溝的保存潛力最高，北海較淺區域以及英格蘭、荷蘭和德國沿海的保存潛力則較低。

المُستخلص

لقد تم القيام باستكشاف عمليات تشكيل وتكوين المواقع لعدد ٥٤٩ من حطام السفن في بحر الشمال وبحر البلطيق وهذا باستخدام بيانات MBES والمتغيرات البيئية/البشرية. تم تصنيف المرحلة الأولى بالسلامة الهيكلية لكل حطام. بينما تم استخدام نمذجة الانحسار والتراجع في المرحلة الثانية وذلك لاستكشاف الارتباطات بين فئة الحطام والمتغيرات البيئية في المواقع الفردية. ومن الجدير بالذكر أن النموذج الأفضل الذي تم اتباعه في هذه الدراسة كان يتألف من العمق والملوحة وكثافة عمليات الصيد وعام التحطم. وفي المرحلة الثالثة، تنبأ التحليل متعدد المعايير بإمكانية توفير الحفظ والحماية عبر مجال النموذج الكامل. وتم تصميم أعلى إمكانات الحفظ في جنوب بحر البلطيق والخندق النرويجي، مع إمكانات منخفضة في المناطق الضحلة في بحر الشمال وقبالة إنجلترا وهولندا وألمانيا.‏

KEYWORDS:

• Shipwreck
• site formation processes
• site environment
• preservation potential

PALABRAS CLAVE:

• pecio
• proceso de formación de sitio
• ambiente del sitio
• potencial de preservación

关键词:

• 沉船
• 遗址形成过程
• 遗址环境
• 保护潜力

關鍵詞:

• 沉船
• 遺址形成過程
• 遺址環境
• 保護潛力

الكلمات الدلالیة:

• حطام السفينة
• عمليات تشكيل الموقع
• بيئة الموقع
• إمكانية الحفظ

Introduction

Shipwrecks provide empirical evidence of ancient trade (Foley et al., 2012), migration (Mentz, 2015), past naval warfare (Church, 2014), fishing (Grenier et al., 2007), technological change (Adams, 2001), maritime trade and transportation networks (Secci et al., 2021) and serve as microcosms of the wider society that produced them (Gibbins & Adams, 2001). In some instances, they are tangible monuments to historic events ranging from individual disasters like the sinking of RMS Titanic (Ballard & Archbold, 1987), while others serve as a memory of global conflicts such as the First (McCartney, 2016) and Second (Macleod, 2006) World Wars.

Shipwreck sites can serve as habitats for marine life, locally enhancing biodiversity and biomass (Balazy et al., 2019). On the other hand, they can pose environmental pollution threats, especially if they contain hazardous substances such as oil, chemicals or munitions (Ndungu et al., 2017; Szafranska et al., 2021). Moreover, marine spatial planning issues often arise as wrecks interfere with ever-increasing offshore sub-sea engineering efforts, such as the installation of cables, pipelines or wind farms (Papageorgiou, 2018). Shipwrecks can also be seen as geohazards and obstacles to marine engineering, as they can pose risks to navigation, may contain unexploded ordnance, and need to be avoided when planning developments. Thus, beyond their archaeological, historical and cultural significance their presence on the seabed incurs several problems and conflicts of interest (Firth, 2018). Common to all these issues is how cultural and environmental forces have, and will, affect the integrity and preservation of these sites in the past, present and future.

These forces can be assessed through the study of site formation processes, which are defined as ‘the factors that create the historic and archaeological records’ (Schiffer, 1987). Shipwreck site formation research encompasses the processes by which wreck sites are created and preserved (Keith, 2016) and provide valuable insights into the factors that influence the formation and preservation of shipwrecks and their associated artefacts and remains. Pre-depositional formation processes, which include the nature of the wrecking process itself, are important in understanding how and in what state a ship is deposited on the seabed. Post-depositional site formation processes are those that affect the long-term preservation of a site. These are commonly divided into natural processes, caused by the ambient physical, chemical and biological processes occurring in the marine environment (N-transforms) and cultural processes, related to human impact (C-transforms) (Oxley & Keith, 2016).

One of the pioneers of shipwreck site formation research was Keith Muckelroy, who attempted to systematically categorise site formation processes and those events which change a site over time (Muckelroy, 1976). He recognised that the site’s environment is ‘the one determinant external to the basic wreck system’ (Muckelroy, 1976, p. 281) and the processes, or ‘extracting filters’, that lead to the loss of material including: the wrecking incident, salvage operations and the disintegration of perishables. Subsequently, after a ship is wrecked, ‘scrambling devices’ affect the archaeological material that is displaced (or scrambled) on the seabed. Before Muckelroy’s work, most nautical archaeological projects concentrated on wreck sites in the Mediterranean and conclusions about types of seabed favouring shipwreck preservation were drawn from individual qualitative case studies (Dumas, 1962; Frost, 1972) rather than systematic assessment of multiple sites.

Muckelroy (1977), developed his work further by establishing a five-tier classification system for wreck site preservation based on 20 wooden wrecks off the British Isles. This incorporated various environmental variables and used statistical assessment of factors influencing the modern-day condition of wreck sites. Significant correlation was observed between the seabed sediment type and the preservation state of shipwrecks. Shipwrecks in Class 1 were the best preserved and were often found in seabeds that were covered with fine-grained sediment, promoting early burial of a ship structure. On the other end of the scale, Class 5 sites were the poorest preserved, containing only heavy metal and stone artefacts with ship structure having been ‘smashed up by natural forces’. However, he does note that there are many intermediate environments resulting in various grades of wreck site preservation. The study concluded that more research into the subject is required and proposed that every archaeological report should include assessment of environmental data relevant to a studied wreck site.

However, as Gregory (1992) and Maarleveld (1995) noted, although Muckelroy’s hypothesis was frequently cited in subsequent research (Hardy, 1990; Nash, 1990; Owen, 1991), it was not significantly developed upon (Hunter, 1994; Lenihan et al., 1981; Parker, 1995). Furthermore, the feasibility of statistical analysis was also questioned by researchers that had access to substantially larger datasets of shipwrecks (Gibbins, 1990). Cederlund (1983) attempted to define the factors which determined the preservation of 600 shipwrecks in Swedish waters, including geographical appearance of the coast; the situation of the wreck in relation to the water surface (depth); type of seabed (geographical (hard, sandy) or geological (sedimentation, erosion)), height differences on the bottom (movement of sand), existence of ice, currents, heavy wave movements, and special biological conditions (presence of Teredo spp).

Nevertheless, Gibbins (1990) notes that classification systems can contribute to the formalisation of ideas and theories within a framework, and site assessment studies provide a fundamental role in improving the viability of site classification systems. However, few such studies have been carried out. Parker (1992) outlined that of 1189 sites recorded and reported in the Mediterranean, as many as 25% have been so briefly reported that neither condition nor depth is known. His classification of wrecks from the classical period based on depth and condition (perfect, coherent, scattered, hull only, unknown) served to emphasise that those wrecks which were well preserved lie in deeper waters.

These early advances were predominantly conceptual works and in the 1990s there was renewed interest in site formation processes from a cultural resource management perspective (Oxley, 1998; Oxley & Gregory, 2002). This was due, in part, to international legislation advocating for in situ management of underwater cultural heritage (UCH). Within Europe this was galvanised by changes within developer-led archaeology (CETS-143, 1992) and by the ICOMOS Charter (ICOMOS, 1996). This was followed by a changing interest in multiple case-studies, validating earlier site formation theory and models. Ward et al. (1999) extended site formation models, by adding key components related to changes in net sediment budget. Gibbs (2006) added cultural processes which comprise anthropogenic influence on the integrity of wreck sites. Quinn (2006) highlighted the importance of assessing seabed scour developed around wreck sites and its role in shipwreck preservation. This also came at a time of rapidly progressing marine remote sensing techniques like multibeam echosounder surveying, allowing effective detection, characterization and monitoring of underwater sites with no physical disturbance (Plets et al., 2011; Westley et al., 2019). The remote characterization and monitoring of change at individual sites (Brennan et al., 2016; Quinn & Boland, 2010) was extended to including environmental factors to understand physical formation processes in control of site preservation (Fernández-Montblanc et al., 2016; Geraga et al., 2020; Majcher et al., 2021; Quinn & Smyth, 2018; Smyth & Quinn, 2014). A holistic approach to understanding the formation of an individual wreck site was proposed by Vega-Sánchez and Herrera (2022) who used agent-based modelling. In the study site formation processes were considered as a complex adaptive system where multiple factors lead to a distribution of a wreck at the seabed. The authors consider a variety of pre- and post-depositional factors influencing the site, from shipbuilding plans, shipwrecking event witness accounts to newly acquired in-situ data. Nevertheless, most of the investigations to date focus on highly detailed investigations of individual sites, rather than aiming to establish methods for the preliminary assessment of UCH in wider spatial extents.

UNESCO estimate there are more than 3 million shipwrecks worldwide (Croome, 1999). The UNESCO Convention on the Protection of Underwater Cultural Heritage, like the earlier Valetta and ICOMOS conventions, encourages in-situ preservation as the preferred option for continued preservation of underwater sites that are over 100 years old (UNESCO, 2001). Many metal shipwrecks, including First World War vessels, now fall under this definition. Metal wrecks are important from many other perspectives; for example, they are often referred to as potential pollution sources. Therefore, there is an urgent need to develop desk-based approaches for assessment of the most significant site formation processes for many sites at once.

Gregory et al. (2012) formed a baseline process-based approach to control, assess and achieve in-situ preservation, and one of its central points is the assessment of the most significant chemical, biological and physical threats to the site. Understanding the formation of underwater sites also enables more accurate archaeological interpretation. Assessing environmental factors which may have influenced sites, not only leads to a better understanding of events that have caused the sinking of the vessel but also help to predict pathways of future wreck deterioration. This in turn can enhance understanding of pollution hazards and address UCH management needs.

In this study, baseline post-depositional site formation processes for over 500 shipwreck sites located in the North and Baltic Seas are investigated. Multibeam echosounder bathymetry data are integrated with multiple environmental and anthropogenic variables sourced from open-data portals. The observed seabed distribution of the sites is initially manually classified using a revised Muckelroy (1977) approach, followed by statistical correlation analysis to determine the most significant variables influencing site formation. This study paves the way for big-data driven approaches to assess predominant basin-scale post-depositional site formation factors acting at underwater heritage sites.

Data and Methods

Multibeam Echosounder Data

High-resolution RESON SeaBat 7125 multibeam echosounder (MBES) bathymetry data were collected at 549 shipwreck sites during multiple surveys in the Danish Straits, the North and Baltic Seas and nearby Scapa Flow between 2015 and 2019 (Figure 1). The EIVA software suite was used for raw data processing and the Navimodel Producer v.4.5.6 application for data exports. The data were vertically reduced to chart datum and corrected for refraction using sound velocity profiles (Valeport, MiniSVS) acquired during the respective surveys. Subsequently MBES point clouds were gridded, with excessive noise manually cleaned. The data were then exported to digital elevation models (DEM), with spatial resolution ranging from 0.2 to 0.3 m.

Figure 1. Overview map representing the study sites (Authors).

Figure 2. A step-by-step block diagram of the workflow designed for the study (Authors).

Environmental Data

To investigate the influence of formation processes on the sites, environmental and cultural variables were extracted from online data portals providing basin-scale, open-source oceanographic datasets. A summary of the data mined from the open-source portals is provided in Table 1.

Table 1. Summary of the open-source environmental variables.

Dataset	Symbol	Source	Reference	Resolution	Unit	Averaging period
Sediment data	–	EMODnet	EMODnet (2023)	Multiscale	Folk7 classification	-
Near-seabed temperature	T	BioOracle v.2	Assis et al. (2018)	∼9.2 km at the equator	^oC	2000–2014
Near-seabed salinity	S	BioOracle v.2	Assis et al. (2018)	∼9.2 km at the equator	PSS	2000–2014
Near-seabed dissolved oxygen - EMODnet	DO	EMODnet	EMODnet (2023)	∼6.2 km	µmol/l	1982–2015 (Baltic), 1977–2017 (North Sea)
Seabed kinetic energy levels – waves and currents	KEw and KEc	EMODnet	EMODnet (2023)	4.5 km	J	Integrated over 3 years. Average of annual percentiles 90th.
Max. near-seabed current velocity	Vmax	BioOracle v.2	Assis et al. (2018)	∼9.2 km at the equator	m/s	2000–2014
Zero-crossing wave period – Baltic	Tz	Copernicus Marine Service	Copernicus Marine Service (2021b)	∼1.8 km	s	2010–2020
Significant wave height – Baltic	Hs	Copernicus Marine Service	Copernicus Marine Service (2021b)	∼1.8 km	m	2010–2020
Zero-crossing wave period – Baltic	Tz	Copernicus Marine Service	Copernicus Marine Service (2021a)	∼2 km	s	2010–2020
Significant wave height – Baltic	Hs	Copernicus Marine Service	Copernicus Marine Service (2021a)	∼2 km	m	2010–2020
Fishing intensities	–	EMODnet	EMODnet (2023)	5.5 km	mW fishing hours	2015–2018
Fishing vessel density	–	EMODnet	EMODnet (2023)	1 km	hours/km^2 per month	2017–2021

Environmental variables were sourced from the BioOracle, Copernicus Marine Service and EMODnet portals. Among the data downloaded from the BioOracle portal (Tyberghein et al.; Assis et al., 2018), near-seabed temperature and salinity represent 14-year averages, while maximum near-seabed current velocity shows long-term extremes in the same period. The BioOracle data were originally sourced from pre-processed global ocean re-analyses, namely ARMOR; Global Observed Ocean Physics Reprocessing, and ORAP; Global Ocean Physics Reanalysis provided by ECMWF (European Centre for Medium-Range Weather Forecasts), which integrate satellite and in situ observations across regular two- and three-dimensional spatial grids. The significant wave height and period data mined from the Copernicus Marine Service (2021a, 2021b) represent 10-year averages. The wave models used for the Atlantic Ocean and the Baltic Sea are WAVEWATCH III (Copernicus Marine Service (2021a) and WAM (Copernicus Marine Service (2021b), respectively, with the atmospheric forcing given by the ECMWF ERA-5 Numerical Weather Prediction reanalysis.

Data sourced from the EMODnet portal (EMODnet, 2023) comprise sediment substrate data represented as Folk 7 classes, with class 1 being the finest and class 7 the coarsest sediment type (Folk, 1954). The sediment data were provided by the EMODnet project’s partners in Europe and generalized to a broad, continental-shelf scale (EMODnet, 2023a). Additionally, near-seabed dissolved oxygen datasets are included for the Baltic (EMODnet, 2023b) and North Sea (EMODnet, 2023c) derived from observational data from SeaDataNet/EMODnet Chemistry Data Network. They represent long-term averages calculated over 33 years (Baltic) and 40 years (North Sea) using data-interpolating variational analysis (DIVA). Seabed kinetic energy exerted by currents and waves, modelled using a combination of oceanographic wave and current circulation models (EMODnet, 2023d, 2023e) was integrated over 3 years and average values of annual 90th percentiles were extracted. Some variables are derived from the mined data. The degree of saturation of oxygen (and thus its activity) is calculated as the measured near-seabed dissolved oxygen divided by the maximum solubility of oxygen at a given salinity and temperature (Garcia & Gordon, 1992). All these environmental variables represent basin-scale modelled data which were cross-validated using in-situ quality-controlled data.

Potential anthropogenic influence on wreck site formation was assessed through the inclusion of variables representing fishing effort. EMODnet (2023) provides a variety of quantified fishing activities in the Baltic and North Seas. For the scope of this study, only the fishing methods entailing direct contact of fishing gear with the seabed were used. Bottom seine, beam and otter trawl fishing data were downloaded from the EMODnet portal, representing fishing intensities averaged over 3 years. The fishing intensity data are aggregated in a polygon grid using C-squares geocoding (EMODnet, 2023f). They are only available for vessels more than 12 m long with vessel monitoring systems (VMS) and vary in terms of spatial extent between the selected gear types. Due to this, and to represent total fishing intensity, bottom seine, beam and otter trawl fishing intensities were summed using GIS.

GIS Database

ArcGIS Pro 3.0.1 was used to integrate the datasets used in the study. MBES data-derived wreck site DEMs, environmental raster and vector data were brought into a GIS project and geodatabase structure. A point shapefile was subsequently created with individual points added for every site. The shapefile’s attribute table was first populated with the wreck description data. Subsequently, the environmental data were extracted for individual wreck points using the spatial analyst toolbox’s ‘Extract Multi Values to Points’ tool. In this way, the shapefile’s attribute table inherited the environmental variables (Table 1), allowing for further statistical analysis (section 2.5).

Wreck Descriptions and Preservation Classes

Muckelroy (1977) separated 20 wreck sites around the British Isles into 5 classes of decreasing preservation state. Class 1 comprised sites that were coherent and well preserved; Class 2 – nearly coherent sites where materials of all types were preserved; Class 3 – a fair amount of organic materials were preserved, with wide range of artefacts found, but no ship structure visible; Class 4 – no organic material preserved, but wide range of artefacts were preserved; Class 5 – smashed up by natural forces, but the most durable objects survived.

For this investigation, the classification system of Muckelroy (1977) was simplified. Four classes were used to grade shipwrecks based on their observed seabed distribution: Class 1 – coherent, well preserved; Class 2 – coherent, but incomplete; Class 3 – scattered ordered, recognizable as a ship; Class 4 – scattered disordered. The classes were manually assigned based on visual interpretation of the high-resolution MBES data. The classification was performed three times, each time by a different person. This resulted in different classes being assigned for some wrecks. Additional examination of the data was performed in order to assign a final class. All the sites which did not fit in any of the classes were assigned Class 5 – comprising wrecks which were known to have been culturally influenced (salvage, fishing), or their appearance on the seabed did not conclusively represent the remains of a shipwreck. Furthermore, information regarding construction and sinking date, hull materials, past salvage activities and cause of sinking were initially incorporated in the study for the known, identified shipwreck sites totalling 170 out of 549 positions. Harmonisation of these variables with the assigned preservation classes was attempted in order to reduce ambiguities related to the initial conditions of the wreck on the seabed.

For each wreck site, the data were inspected in two forms to further reduce classification subjectivity: as DEM surfaces and as point-clouds. DEM data were rendered in ArcGIS Pro 3.0.1 with multi-directional hillshading and raster colormap enhancement. These techniques visually enhanced the sites, consequently allowing more objective classification. Processed point-cloud data were brought into CloudCompare v. 2.12.4. The ShadeVis ambient occlusion plugin and the EyeDome Lighting OpenGL shading were applied to improve the rendering of the 3D point-cloud data. This approach for inspecting MBES shipwreck data was proven effective by Majcher et al. (2021) who used it to assess structural deterioration of historic shipwreck sites.

Correlation Analysis

Possible correlations between wreck class and various natural and anthropogenic environmental parameters were investigated through multivariate regression models. The focus was on wreck classes 1-4, i.e. wrecks that were still visible on the seabed. The dependent variable (wreck class 1-4) is not numeric but rather an ordered categorical variable. The categories have a meaningful order, but the intervals between categories are not necessarily equal, which constrains the choice of regression models. An ordered logistic regression model was chosen for the analysis, carried out using the statistical software package R ver 4.2.2, using the polr function from the MASS package. Ordered logistic regression models assume that the probability (or rather the logarithm of the odds: logit) of an observation falling into a specific category is a linear combination of the independent variables. In our case this allowed multiple regression analysis with a number of variables from Table 1 and calculation of the probability of wrecks being in wreck class 1-4 for some given conditions. The performance of the model was evaluated using the confusion Matrix function of the caret package that compares the predicted class (i.e. the class with the highest probability) to the observed class for all the wrecks and calculates the percentage of wrecks that are categorised correctly by the model. Furthermore, the output of the model was compared to a ‘no information rate’ to see if the model makes better predictions than no model or pure chance.

Different combinations of the variables in Table 1 were tested, starting with complex models with a large number of variables, and then making a stepwise removal of variables that didn’t contribute significantly to the explanatory power of the model. The aim was to make a simple and robust model that can calculate wreck classes for the known wrecks based on local environmental conditions. The best model was subsequently used for making a weighted overlay analysis, giving a modelled predicted class for the whole region.

Multi-Criteria (Weighted Overlay) Analysis

Weighted overlay preservation prediction modelling was conducted based on the output of the correlation analysis. The logistic ordered regression provided a linear combination of the most important environmental variables that describes the odds of known wrecks being in class 1-4 (on the logit-scale). Thus, for some given environmental conditions it was possible to calculate the probability of each class, and three intercept values representing the limits where the highest probability goes from class 1 to 2, from class 2 to 3, and from class 3 to 4. This was used to create a weighted overlay model covering the whole region, where the environmental variables were first multiplied by their respective coefficients and then summed. Subsequently, the overlay raster was classified based on intercept values grouping the raster data into expected preservation classes.

Results

Wreck Classification

Of the 549 wrecks classified using the modified Muckelroy scheme, 76 are classed as well preserved and coherent (Class 1), 126 as coherent but incomplete (Class 2), 157 as scattered but ordered (Class 3), 93 as scattered disordered (Class 4), and 97 as being either heavily modified by anthropogenic activity or impacted so much by external forcing that no wreck remains are visible (Class 5) (Table 2; Figures 3 and 4).

Figure 3. Examples of digital elevation models (left hand column) and point clouds (right hand column) of wreck classes 1–4 (Authors).

Figure 4. The spatial distribution and classification of the 549 wrecks classified in the study (Authors).

Table 2. Numbers of sites within their assigned classes

Class	Description	Count	Example
1	Coherent	76	Fig. 3a,b
2	Coherent, incomplete	126	Fig. 3c,d
3	Scattered, ordered	157	Fig. 3e,f
4	Scattered, disordered	93	Fig. 3g,h
5	No wreck/significant human impact	97	–

Environmental Variables

Although the North and Baltic Seas make up an interconnected system, they are often treated as three separate geographic areas; the North Sea, the Danish Straits, and the Baltic Sea (Figure 1). The North Sea is a shelf sea, with an average depth of 100 m, and is open to the North Atlantic between Scotland and Norway in the north, and through the English Channel in the south-west (Figure 5a). The Baltic Sea has shallow and narrow connections to the North Sea through the Danish Straits. Comprising a number of basins connected by relatively deep channels separated by shallow sills, the Baltic Sea averages 55 m, with typical basin depths between 100 and 250 m (Figure 5a).

Marine sediments in the North Sea are dominated by sand, with mud, sandy mud and muddy sand located in low-energy deep areas (Figure 5b). Coarse and mixed sediment is present to the west of Jutland, through the Danish Straits and into areas of the Baltic Sea, where higher energy hydrodynamic regimes dominate (Fig 5b). North-east of Bornholm in the Baltic Sea, the seafloor is dominated by mud and sandy mud, while south of Bornholm sand dominates. Good correlation is noted between depositional environments (Fig 5b), current velocity (Fig 5c) and wave period/height (Figures 5d and e), as the relative kinetic energy of the ocean environment is reflected in depositional environments.

Bottom salinity in the North Sea is close to 35, except near river outlets in the south-east and along the northeast where the Baltic Sea and large rivers lower the surface salinity to 25 (Figure 5f). Approaching the Baltic Sea through the Danish Straits, salinity drops to 15, and in the Baltic Sea salinity ranges from 7-10 in the south to 5-6 in the north (Figure 5f). Bottom temperatures in the Baltic Sea are generally low compared to the North Sea (Figure 5 g), with sea ice forming in the Baltic every year. Anthropogenic nutrient enrichment and rising water temperatures can lead to reduced oxygen concentrations. This is particularly noted for the Baltic Sea, now recognised as one of the largest dead zones in the world (Conley et al., 2009), with nutrient input from agricultural fertilisers leading to oxygen depletion and eutrophication (Figure 5h).

Figure 5. Environmental variables: (a) depth, (b) seabed substrate, (c) near-seabed maximum current velocity, (d) zero-crossing wave period, (e) significant wave height, (f) salinity, (g) temperature, (h) near-seabed dissolved oxygen (Authors).

Bottom Fishing Intensities

Fishing is one of the dominant anthropogenic activities affecting marine systems (Halpern et al., 2008) and global concern about adverse effects of bottom trawling on seabed habitats (Eigaard et al., 2016) and UCH (Brennan et al., 2016; Krumholz & Brennan, 2015) is increasing. Bottom trawling intensities for the North and Baltic Seas are presented as heat maps for beam trawling (Figure 6a), otter trawling (Figure 6b), bottom seines (Figure 6c), and cumulative bottom fishing intensity (Figure 6d). Bottom trawling is widely distributed over the study area, with otter trawling prevalent in both the North and Baltic Seas (Figure 6b), while beam trawling and seine fishing is concentrated in the North Sea (Figures 6a and 6c). Cumulative impact shows high fishing effort in areas in the Norwegian Trench, in large areas of the northern and southern North Sea to the west of Jutland, in the Danish Straits, and sustained effort around Bornholm in the Baltic Sea.

Figure 6. Fishing intensities: (a) bottom beam trawling, (b) bottom otter trawling, (c) bottom seines, (d) total bottom fishing intensity. Grey areas represent lack of data (Authors).

Correlation Analysis

Figure 7 shows a boxplot where different variables are shown grouped by wreck classes. It indicates that wreck class cannot be explained by one single parameter, as there is a pronounced overlap between the different classes – thus for instance the best-preserved class 1 wrecks are found at depths from 6 to 172 m, while the degraded class 4 wrecks are found at depths between 5 and 142 m. Correspondingly, ordered logistic regression models using only one of the variables in Table 1 make very poor predictions that are not significantly better than no model (P > 0.05).

Figure 7. Boxplots for environmental conditions at the wreck sites, grouped by the wreck classification. Numbers of wrecks in each class are given in Table 2. The boxplot to the lower right gives the year of sinking /wrecking which is known for 160 of the classified wrecks (Authors).

Things improve if more parameters are included. A regression model that includes depth, current velocity, sediment type, kinetic energy of waves, temperature, salinity, dissolved oxygen, fishing intensity and the year of wreckage classifies 52% of the wrecks correctly and is significantly better than no model (P = 0.009). However, it is also a complex model and there is a significant risk of overfitting, especially as some of the variables are not independent. In order to simplify the model, one parameter at a time was removed and the model tested for reduced explanatory power. In this stepwise manner, the model was reduced to a simpler and more robust solution based solely on depth, salinity, bottom fishing intensity and year of wreckage. This simpler solution classifies 54% of the wrecks correctly and is significantly better than no model (P = 0.0002).

Removing any of these four variables from the model results in significantly poorer predictions. Still, it is of interest to remove the year of sinking / wrecking and only look at external factors such as the environment and influence from fishing. This increases the number of wrecks in the analysis, as the year of wreckage is only known for 170 out of the 452 classified wrecks. A simple model based on three variables (depth, salinity and fishing activity) classifies 40% of the wrecks correctly, which is actually not much better than no model (P = 0.07). This demonstrates that intrinsic parameters such as wreck age, type, construction etc. obviously play important roles on what is preserved today.

The reduced model gives a linear equation that allows us to calculate the probabilities of wrecks being in class 1-4 under some given conditions: (1) $\begin{array}{rl}& \mathrm{Logit}\left(\mathrm{Clas}{\mathrm{s}}_i\right)={\alpha }_{\mathrm{i}}+0.0089\cdot \mathrm{Depth}\left[\mathrm{m}\right]\\ & +0.0895\cdot \mathrm{Salinity}\left[\mathrm{PSU}\right]+0.00275\\ & \cdot \mathrm{bottom}\mathrm{fishing}\mathrm{activity}\left[\mathrm{mW}\right],\\ & \mathrm{where}\mathrm{intercepts}{\alpha }_{1|2}=0.727,{\alpha }_{2|3}\\ & =2.329\mathrm{and}{\alpha }_{3|4}=3.999\end{array}$(1)

The probability of being in class 1 (well preserved) increases at low numbers i.e. for greater depths (as these are negative), lower salinities, and lower fishing activity. This is illustrated in more detail in Figure 8, showing the probability of different wreck classes under some given conditions. It is based on our wreck database (and ignoring intrinsic parameters) where low salinities (upper line) are only found in the Baltic, whereas in the North Sea all sites have a salinity close to 34. Although the model is not perfect, it is a best fit for our data. Furthermore, the linear equation can be used to make predictions for the whole area, and not only for the individual wreck sites (Figure 8).

Figure 8. Estimated probabilities of different wreck classes under varying conditions, using a reduced model that looks solely at depth (m), salinity S (PSU) and bottom fishing intensity (mW fishing hours) (Authors).

Multi-Criteria (Weighted Overlay) Analysis

Multi-criteria analysis was informed by the results of the statistical analysis, with the three input variables resulting from the simplified model: salinity, depth and bottom fishing intensity. Weighted overlay analyses show highest preservation potential in the southern Baltic Sea and in the deeper water of the Norwegian Trench (Figure 9). Lowest preservation potential is modelled for the shallower areas of the southern North Sea, off the east coast of England, the north coasts of The Netherlands and Germany, and west of Jutland.

Figure 9. Map with the classified wrecks and the predicted preservation potential/class overlay raster. Points coloured with shades of red represent the classes manually assigned for shipwrecks (1 to 4), while the grey points represent class 5 sites, excluded from the classification. Grey areas result from insufficient data to predict preservation potential (Authors).

Discussion

In this study, we demonstrate that the structural integrity of shipwreck sites (n = 549; year of loss 1801-2003) in the North and Baltic Seas is correlated with depth, salinity, bottom fishing intensity and time. In terms of wreck preservation, the southern Baltic Sea and the Norwegian Trench offer the highest potential, while the shallower areas of the North Sea and off England, the Netherlands and Germany are least favourable.

Wreck Classification and Regression Analysis

Undoubtedly MBES bathymetry data provide an informative overview of the distribution of those parts of a shipwreck lying on the seabed and in the water column (Figure 3). Furthermore, MBES data yields information on seabed morphology and dynamics prevalent on and around sites (for example Figures 3 c and 3 g show the effects of scour around the wrecks). However, classifying the state of preservation and the preservation potential of a given wreck site using these data is problematic due to the inherent subjective nature of visual assessment and interpretation.

Although correlation between wreck class and environmental parameters was established through regression analysis (Table 1), the fit was not perfect. This is not so surprising, as the wreck class was determined by decay accumulated over a long time-period and by different processes, while the environmental variables we used represent the current conditions and hence the contemporary decay rate. For instance, it makes sense that including time on the seabed would improve our models significantly. The identified wrecks in the model sank between 1801 and 2003, and obviously a ship that sank over 200 years ago is likely to have undergone more deterioration than a ship that sank only 20 years ago. This parameter, together with other wreck-specific data, should be included in a next iteration of the regression analysis in order to separate more clearly the effects from external stressors and intrinsic variables. Focussing on wreck classes limits statistical correlation as the dependent variable (wreck class 1-4) is not continuous but rather an ordered categorical variable which constrains the choice of regression models. Thus, future work will attempt to improve on objectively quantifying the state of wreck preservation in order to facilitate better modelling of decay.

Wreck Environments

Comparison to the original study by Muckelroy (1977) shows some fundamental differences. He found sediment type and offshore fetch to be best correlated to the preservation of 20 wooden wrecks off Britain and Ireland, whereas we find depth, salinity, bottom fishing intensity and time the most important. It may be surprising that sediment type does not improve our models significantly, but this could be due to relatively poor spatial resolution in the EMODnet-derived substrate layers used in the modelling (see discussion below) and/or the lack of variation in our data, as most of the wrecks in our database are found in sand (Figure 4). The kinetic energy at the sites (including the wave height, wave period and current velocity) did not contribute significantly to the models either, which is also surprising. However, absence of evidence is not evidence of absence, and (some of) the effects from variations in energy may already be included in other model parameters such as depth and sediment type, as preservation is observed to be better at depth where the energy level will generally be lower and substates finer. Variations in dissolved oxygen could also be important, but when we compensated for variations in salinity and temperature to calculate the oxygen saturation the apparent variations in Figure 4 turned out to be much lower and by far most of the wreck locations had close to 100% oxygen saturation.

Data Resolution

A limitation with the approach used here stems from using data sets with different spatial and temporal resolutions (Table 1). For example, the MBES data used to classify the integrity of the wreck sites has a resolution of centimetres (Westley et al., 2019), whereas the BioORACLE, COPERNICUS and EMODnet datasets employed in the correlation analysis and weighted overlay modelling have grid cells in the order of kilometres (Figure 10). For instance, the global-scale BioORACLE model has a spatial resolution of 5 arcmin, which corresponds to 9.2 km at the equator (Assis et al., 2018) and 5.4 km at a latitude of 54°N.

Figure 10. (a) 5 arcmin cells from the BioORACLE model overlain on 4 m resolution multibeam echosounder data. (b) Zoomed view of the hillshade multibeam echosounder data showing the complex erosional and depositional patterns around two First World War shipwrecks (Authors).

This point is illustrated by plotting BioORACLE model output cells at 5 arcmin resolution on top of the results of a multibeam survey carried out over two First World War shipwrecks (Figure 10a). The multibeam-derived DEM has a cell size of 4 m, typical of the resolution of MBES data available for download from national open-data portals (e.g. UK Admiralty Marine Data Portal; Ireland’s INFOMAR). Intricate patterns of erosion and deposition at the two wreck sites are clearly seen in the multibeam data (Figure 10b); complexity in the wreck system that simply cannot be captured by the two BioORACLE grid cells covering the sites. Typically, the resolution of these global- and shelf-scale models are two orders of magnitude poorer than multibeam bathymetry data.

On a related point, the multibeam echosounder data collected from the 549 wreck sites used in this study were collected over different months and years, and the BioORACLE, COPERNICUS and EMODnet data are averaged over long time periods, so temporal offset is also a limitation in the analysis. Not only are model components available at different resolutions, but the spatial and temporal scales at which key environmental and anthropogenic processes take place and interact are also highly variable. Consequently, correlation analysis or multi-criteria analysis conducted using these model ingredients can suffer from oversimplification, and result in smoothing or under-representation (or in the worst-case exclusion) of key variables influencing the system.

MBES Data and Open Data Portals

Ideally MBES and similar high-resolution data types should be used when they are available (e.g. Majcher et al., 2021; Manders 2009) to understand the formation of shipwreck sites. Nevertheless, in the absence of these types of high-resolution data, shelf-scale and global-scale open-data model outputs and their derivatives allow us to model preservation potential at national- and shelf-scales to inform management plans based on real world data. As sensors continue to develop, and national and international seabed mapping initiatives gain momentum (Mayer et al., 2018) and become more widely available to the community, reliance on coarse-resolution data and low-resolution model outputs will naturally decrease.

Future Work

Future work by the authors will build upon this initial study within the recently started five-year (2022–2027) project ENDURE (www.endureerc.com). The overarching aim of ENDURE is to propose a novel framework to cultural heritage management, based on the state of entropy of sites (DeSilvey, 2021). Continued refinement of the model presented here will incorporate more MBES surveys of shipwreck sites over a larger geographic area and the addition of other environmental and anthropogenic factors from open data portals. The refined models will be validated by survey, diving and monitoring of environmental parameters on select sites. The models will be further supplemented with data relating to decay rates of archaeological materials (especially metals and wood) and an improved understanding of the role of artefact, design form and function (Nutley, 2009) on the processes and rates of physical, chemical and biological decay. The hypothesis is that these factors play a significant role in the rate of decay and hence survivability of shipwrecks and artefacts underwater.

Conclusions

• Weighted overlay modelling and regression analysis demonstrates that the structural integrity of wreck sites in the North and Baltic Seas is correlated with depth, salinity, bottom fishing intensity and time.

• The study shows the potential of integrating high-resolution remotely sensed data with environmental information from open data portals in predicting preservation potential of UCH sites across large geographic areas.

Author Contributions

DG secured funding and with RQ formulated the ideas behind the research. MD collected and post processed multibeam echossounder data. DG, GN and JM classified the wrecks based on multibeam echosounder data. RQ and JM collated data from the open access data portals and generated the GIS. JM and HM conducted statistical analyses and produced images for the publication. All authors contributed to the writing and editing of the manuscript

Permissions Statement

Access to the multibeam echosounder data was provided to the authors by Rasmus Normann Andersen, CEO of JD Contractors. The MBES data used within this paper were collected on numerous academic research trips, commercial surveys, and as data exchange with other archives/survey companies. Relevant regulations for collection of data were in line with national and international legislation of the waters where surveys were conducted. Exact positions of wrecks are not given due to issues of commercial interest and national security. Furthermore, data were collected adhering to the principles of the UNESCO convention in as much it was non-destructive and non-invasive and no shipwrecks were disturbed during these surveys. All imagery has been produced by the authors.

Acknowledgments

The authors would like to thank Rasmus Normann Andersen of JD Contractors for granting us permission to use the multibeam echosounder data presented in the paper. The work was supported by the Danish Ministry of Culture under Grant FPK.2017-0037 and the European Union (ERC, ENDURE, 101053993). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The authors report there are no competing interests to declare.

Disclosure Statement

The authors declare no conflict of interest.