Digital Elevation Models and Orthomosaics of the Dutch Noordwest Natuurkern Foredune Restoration Project

1. Summary

Coastal dunes are intrinsically valuable ecosystems along most of the world’s sandy beaches, whose dynamics is determined by plant-geomorphic feedbacks [1,2,3]. They also contribute directly or indirectly to human well-being by protecting against coastal hazards and providing drinking water and space for recreation [4,5]. Human activities related to, for example, urban expansion and development for tourism increasingly damage, degrade, or destroy dune habitats around the world. To reverse this trend of negative transformation, restoration measures are being taken in more and more dune areas to create conditions that allow the dunes to naturally return to a previous biogeomorphic trajectory e.g., [6,7,8,9,10,11]. The number of ways in which dunes can be restored is enormous [12,13,14] and includes efforts to replant sand patches to restore or build dunes e.g., [15,16,17,18], as well as efforts to remove vegetation to reactivate aeolian dynamics and locally reset ecological succession e.g., [6,8,10]. The way in which restored dunes will develop is often highly uncertain due to the long time scales involved (years to decades), the legacy of pre-restoration conditions (e.g., any remaining seeds or rhizomes) [19], the non-linearity of biogeomorphic feedbacks, and changing boundary conditions (e.g., due to global change). Therefore, many restoration projects are adaptively managed based on evidence [20,21]. Monitoring is essential not only to see if the imposed measures are working and if additional actions are needed, but also to "learn by doing". However, since (non-)governmental organisations or consultants are often responsible for monitoring, the resulting data usually remain part of the "grey literature" [12,22] and are not easily accessible to improve scientific understanding of the biogeomorphic feedbacks underlying the changes of restored dunes and to help practitioners develop guidelines to better design future efforts.

The aim of this data descriptor paper is to make available a unique open access data set (CC BY 4.0) of remotely sensed digital elevation models (DEMs) and orthomosaics from the Noordwest Natuurkern (NWNKern) project in Zuid-Kennemerland National Park, located along the west coast of the Netherlands (Figure 1a). The NWNKern project, one of the largest coastal dune restoration projects in northwest Europe, involved the excavation of five artificial trough blowouts (also termed "notches" [8,10]) through the foredune in winter 2012/2013 (Figure 1). The project aimed to reactivate aeolian sand dynamics and geomorphological diversity in order to increase the area of the Natura 2000 habitat type H2120 Shifting dunes along the shoreline with European marram grass (Calamagrostis arenaria, formerly Ammophila arenaria) or "White dunes" and, second, improve conditions for habitat type H2130 Fixed coastal dunes with herbaceous vegetation or "Grey dunes" [23]. The latter is expected, at least in part, from the inland suspension of fine calcareous material; its deposition as a grainfall up to hundreds of metres landward of the foredune [24] can limit soil acidification [25]. Although the project was started to expand and restore coastal dune habitats, it was soon realised that inland transport of windblown sand through the notches would also benefit long-term coastal safety by facilitating elevation gain of the inland dunes and therefore their vertical growth with sea level rise [19]. The data set includes 35 DEMs and 24 orthomosaics collected between 2013 and 2023, and two pre-excavation DEMs from 2008 and 2012. An analysis of geomorphological changes during the first years after notch excavation can be found in [8] and of the data set up to October 2021 in [26]. Preliminary work on the classification of habitats from orthomosaics can be found in [27,28,29].

2. Data Set Description

2.1. Study Site

The five NWNKern notches, labelled N1 to N5 from south to north, are located in an approximately 1 km stretch of the 20 m high established foredune between Bloemendaal aan Zee in the south and IJmuiden in the north (Figure 1a). In recent decades, the beach sediment budget has been moderately positive [30,31], mainly due to the impact of the 2.5 km long IJmuiden harbour moles on waves and currents near the coast and probably also related to various sand nourishments on the beach of Zandvoort aan Zee, 6 km to the south. The foredune borders a low-gradient intertidal beach and an upper beach with incipient foredunes up to 2.5 m high (embryo dunes) on its sea side, and an ≈5 km wide dune system with vegetated parabolic dunes and dune slacks on its land side. Vegetation was removed from several parabolic dunes as part of the NWNKern project [23,32] (e.g., P1-P3 in Figure 1b,c). Most of the national park is part of the European Natura 2000 network of protected nature areas. The most common habitat types in the vicinity of the notches are the aforementioned White (H2120) and Grey (H2130) dunes, together with "Dunes with sea buckthorn Hippophae rhamnoides" (H2160, often together with wild privet Ligustrum vulgare) and "Humid dune slacks" (H2190) [33]. The sediment consists of medium sand, with the median grain size fining from 300–350 $\mu$m on the beach to 250-300 $\mu$m in the notches.

A total of 170.000 m³ of sand was excavated in winter 2012/2013 to create the five notches. As detailed in [8], the notches had a V shape (viewing landward) through the foredune, with the notch floor starting at approximately 6 m above Mean Sea Level (MSL) at the seaward side and sloping linearly to about 9 m +MSL at the landward side. N2 and N1 were the narrowest (50 m) and widest (100 m) notches, respectively, while the cross-dune length of the notch floor was the smallest in N1, N2 and N3 (100 m) and largest in N5 (200 m). The maximum excavation depth ranged between 9 m in N3 to 12.5 m in N5. The three southern notches (N1-N3) are approximately normal to the coast, while the two northern notches (N4 and N5) have a more southwest-northeast orientation (Figure 1b). The notches deliberately had different characteristics (width, depth, orientation, perimeter) in order to investigate which combination of characteristics would result in the largest landward aeolian sand transport.

The Netherlands experiences a temperate oceanic climate. The annual mean wind speed in the North Sea in front of the Dutch coast is about 8 m/s. In winter, the mean wind speed is substantially higher (≈ 9.5 m/s) than in summer (≈ 6.5 m/s). The highest 10-minute average wind speed varies from year to year but is typically in the range of 21 to 25 m/s. The predominant direction of the wind is from the southwest (Figure 1a). Especially storms from the northwest can raise the tidal water level by 1 m or more (storm surge), but dune erosion at the site is rare due to the wide intertidal beach and the wave sheltering by the IJmuiden harbour moles during these conditions. Annual rainfall amounts to about 800 mm.

The site experienced a variety of human activities after notch excavation. Some of these measures were planned in advance given the experience in nearby restoration measures [19,34], other activities happened to be in the study area but were not related to the NWNKern restoration measures. In 2013 and 2014, marram grass was manually removed from the mouth of the notches. Several times, the rubble from the notch floors was mechanically removed (sieving). The rubble is primarily made up of small concrete blocks or bricks that come from bunkers on the Atlantic wall of World War II or from an abandoned bicycle track and beach entrance. Part of the sand blown on the current bicycle track landward of the southernmost notch (N1; Figure 1b) was bulldozed to the beach (≈ 900-1500 m³ per year between 2013 and 2016) or elsewhere on the depositional lobe (e.g., in December 2020), but this practice was subsequently abandoned and the track has since been buried under several metres of sand. The N1 notch is a public beach entrance ("Strandslag Kattendel"), but the other notches are closed to the public. Beachgoers could walk freely between the incipient foredunes on the upper beach, but this region was fenced off in spring 2023 to create roosting and breeding places for coastal birds (project "Groene Strand"). A club for catamaran sailors and surfers was located at the base of the foredune south of N1 until summer 2017. Finally, in summer 2019, a section of incipient dunes north of the northernmost notch (N5) was cleared to detonate several World War II bombs found near Schiphol airport. The absence of these dunes caused a significant (windblown) supply of sand on and beyond the foredune in subsequent autumn-winter periods. The deposits were so thick in places that marram grass and small sea buckthorn shrubs were completely buried. In subsequent years, vegetation has recolonised most of the foredune, but the embryo dunes have not fully recovered by the end of 2023.

2.2. Morphological Data

The core of the data set consists of 37 digital elevation models (DEMs) and 24 orthomosaics collected between May 10, 2008 and September 25, 2023. The survey dates are provided in the Open Document Spreadsheet file NWNKern_Surveys_ErrorStatistics.ods, Sheet NWNK, Column B as YYYY-MM-DD, where YYYY, MM, and DD are the year, month, and day of the survey, respectively. The first two DEMs (2008-05-10 and 2012-02-20) show the situation before the notch excavation, the third DEM (2013-01-14) during the excavation, and the other 34 DEMs, starting from 2013-05-01, show the morphological evolution after the excavation work was completed. DEMs were computed from survey data collected with different remote sensing techniques (NWNKern_Surveys_ErrorStatistics.ods, column C): airborne laser scanning (ALS, 14 surveys), lidar from an uncrewed aerial vehicle (UAVLidar, 2 surveys), and UAV-acquired Structure-from-Motion photogrammetry (UAVSfM, 21 surveys). Orthomosaics are available for all 23 UAV surveys and a single ALS survey, as indicated with Y in NWNKern_Surveys_ErrorStatistics.ods, column E; otherwise, an N is provided in this column.

The DEMs are provided as GeoTIFF files named NWNKern_YYYYMMDD_METHOD_Zf.tif, where NWNKern is the abbreviation for the research site, METHOD is ALS, UAVLidar or UAVSfM, and Zf indicates that the file contains gap-filled (f, see Section 3) elevation Z. The horizontal coordinate system is Amersfoort/RD New (EPSG:28992). The horizontal axes are easting (X) and northing (Y), and are orientated east and north, respectively. The unit of measure is metres, m. Elevation is in m with respect to Normaal Amsterdams Peil (NAP), with a positive value above NAP. The value $Z=0$ m NAP corresponds approximately to the mean sea level. The DEMs have a rectangular grid with a resolution of 1 m. The lower left corner coordinates ($X_{\mathrm{ll}},Y_{\mathrm{ll}}$) are (98300, 492900) m, and the upper right corner coordinates ($X_{\mathrm{ur}},Y_{\mathrm{ur}}$) are (99200, 494100) m. Each grid thus measures 900×1200 m; its outline is shown in Figure 1a with the red box and corresponds precisely to Figure 1b. The 900×1200 Z values in each DEM are provided for the centres of the grid cells. The NoData value for Z is $-9999$. Additional cartographic information and grid properties can be read from the metadata in the GeoTIFF files.

The 2013-05-01 DEM is provided as an example in Figure 2a and clearly illustrates the five excavated notches. The morphological change in the available data set is demonstrated with the DEM of Difference (DoD) for 2013-05-01 – 2023-02-08 in Figure 2b, with positive and negative difference values corresponding to deposition and erosion, respectively. The DoD shows that depositional lobes have formed landward of the notches, locally almost 10 m thick and extending to 300 m inland. Sand was also deposited on the foredune remnants between the notches (rim deposition) and on the seaward side of the foredune to the south and north of the five notches. Erosion is concentrated in the five notches, especially the lateral walls (locally up to almost 9 m), and at the reactivated parabola dunes, most notably P2 and P3. Time series of deposition volume $V_{\mathrm{d}}$, erosion volume $V_{\mathrm{e}}$ and net volume change $V_{\mathrm{n}}$ (Figure 2c) were calculated for a polygon (Figure 2b) with its seaward side on the upper beach to exclude changes in beach volume from the calculations and its landward side to capture the evolution of the notches and lobes throughout the interval of available DEMs. The polygon is provided in the repository with the GeoJSON file notches.geojson. Volumes are taken relative to the 2013-05-01 DEM, and therefore the time series show the cumulative morphological evolution after the excavation of the notches. At the end of the time series in September 2023, $V_{\mathrm{d}}$ and $V_{\mathrm{e}}$ were approximately 360,000 and −170,000 m³, respectively, implying a net sand gain since May 2013 of about 190,000 m³ (Figure 2c). The three time series show an approximately linear trend over time, with best fit slopes m of 32,785±1,552, -16,405±1,354 and 16,380±2,283 m³/yr ( ± values are the 95% confidence interval) and correlation coefficients r of 0.99, -0.98 and 0.94 for $V_{\mathrm{d}}$, $V_{\mathrm{e}}$ and $V_{\mathrm{n}}$, respectively.

Orthomosaics are available as 8-bit GeoTIFF files with 3 bands (band 1: red; band 2: green; band 3: blue) and are named NWNKern_YYYMMDD_RES.tif. Here, RES is 005 or 100 for 0.05 or 1.00 m square grid resolution. The corner coordinates in all orthomosaic files are identical to those in the DEM files. With the 1 m resolution the orthomosaic grids are thus identical to those of the DEMs. The 0.05 m resolution files have 18000×24000 pixels and are provided for improved visualisation purposes, as well as for future studies on vegetation dynamics. NoData in the orthomosaics is expressed by the colour white. Examples of orthomosaics are provided in Figure 3.

2.3. Meteorological Data

To advance future studies on the biogeomorphic evolution of the study region, the morphological data set is supplemented with high-frequency times series of wind characteristics and rainfall. Both series span the time period between January 1, 2008 and December 31, 2023 with a resolution of 10 minutes and were obtained from the Royal Netherlands Meteorological Institute Data Platform (https://dataplatform.knmi.nl/).

Wind   Wind observations (Figure 4a) from the weather station (IJmuiden WMO 06225, see Figure 1a; 52°27.733 N, 004°33.300 E) are available in the data repository with the ASCII space-delimited file windIJmuiden.txt. The file has seven columns: year, month, day, hour, minute, wind speed, and wind direction. The time is in UTC +01:00 (Mean European Time), the wind speed is in m/s, and the direction in °N. Both wind variables are at a height of 10 m, and NoData values are indicated with NaN. Also shown in Figure 4a are the 24 storm events in the period 2013-2023. Notable are the four storms in January-February 2022 (named Corrie, Dudley, Eunice, and Franklin), the last three of which occurred within a week (16-20 February). The highest 10-minute average wind speed (31.4 m/s) was measured during the storm on 5 July 2023, named Poly. The wind rose of the observations is provided in Figure 1a.

Rain   Rainfall observations (Figure 4b) are provided for the Wijk aan Zee weather station (WMO 06257; 52°30.316 N, 004°36.182 E, which is in the dunes northeast of the northern IJmuiden harbour mole), with the ASCII space-delimited file rainfallWijkAanZee.txt. The file has six columns: year, month, day, hour, minute, and rainfall intensity in millimetres per hour. As in the wind file, the time is in UTC +01:00 and the NoData value is NaN. The median annual rainfall in the period 2008-2023 was about 810 mm; the highest annual rainfall was recorded in 2023 (1113 mm).

The IJmuiden and Wijk aan Zee stations were chosen for two reasons: they are the two stations closest to the study site where wind and rain are measured by the Royal Netherlands Meteorological Institute, respectively, and they are located near the foredune. Therefore, they can be considered the most representative stations for the study region.

3. Methods

3.1. Uncrewed Aerial Vehicle Structure-from-Motion (UAVSfM) photogrammetry

UAV surveys to acquire true-colour aerial images (RGB) suitable for SfM processing have been conducted on average 2 to 3 times per year since 2013-05-01. The first image data set was collected with a fixed-wing Trimble X100 by a commercial contractor (Kragten). All subsequent surveys were conducted by drone experts from Utrecht University, first with a fixed-wing Easystar I equipped with a 12.1 Mpix Canon Powershot D10 and from July 2019 with a DJI Phantom 4 RTK quadcopter with a 20 Mpix camera. The Phantom 4 has centimeter-level accuracy Real-Time Kinematic (RTK) positioning on-board; spatially interpolated RTK-correction data were provided during surveying by a commercial service provider, 06-GPS.nl. The coordinate scheme for the Phantom 4 camera locations was ETRS89 (EPSG:4258). Before SfM processing, the coordinates were converted to Amersfoort/RD New (EPSG:28992) and NAP elevation with Trimble^® Business Center software. The images collected with the Trimble X100 and Easystar I did not have location information. Two surveys (2020-03-09 and 2021-04-26) had to be abandoned prematurely for logistical reasons and contained part of N1 only. The 2018-04-19 survey was considered completed in the field, but it turned out to miss part of N1 during processing. These three surveys are indicated with N in NWNKern_Surveys_ErrorStatistics.ods, column D (`uncompleted’); otherwise, Y is provided in this column (`completed’). The latter does not imply that the entire grid has elevation values; instead, it implies that the DEM contains all the notches and depositional lobes.

Ground control point (GCP) targets were deployed during each survey to be used in the SfM workflow and to georeference the images. Although Phantom 4 images were directly georeferenced, we found, consistent with previous studies e.g., [35,36], that a small number of GCPs considerably improved the accuracy of image locations, especially elevation. Up to 40 GCPs were placed throughout the study region, with a smaller number (typically 10) during Phantom 4 surveys (compare Figure 3a with Figure 3b). The actual number of placed GCPs is provided in NWNKern_Surveys_ErrorStatistics.ods, column G. The GCPs had circular or checked black-and-white patterns. Their centres were measured with a Trimble^® RTK-GPS receiver, connected to 06-GPS.nl, with a horizontal and vertical accuracy of approximately 0.02 and 0.04 m, respectively. The GCP coordinate scheme was Amersfoort/RD New with elevation with respect to NAP.

SfM processing was performed in Agisoft Metashape Pro (v. 1.8.4; https://www.agisoft.com/), except for the May 2013 survey, for which we do not have details of the workflow. The workflow for all other UAV image data sets was based on the guidelines in [37] and consisted of the following main steps, all performed with an Ubuntu 2.04 Virtual Machine with a single A10 GPU in the HPC SURF Research Cloud:

1.

All images were imported into a new Agisoft project, as well as camera locations and orientations with their associated error estimates, if available. Image quality was calculated and poor quality images, often due to motion blur, were disabled. When applicable, the sea surface was masked in the images. The images were aligned in the high-quality setting, with a keypoint limit of 40,000 and unlimited tiepoints.

2.

Low-quality tiepoints were selected and removed in several consecutive steps. First, obvious outliers were manually selected and removed. Second, tiepoints were selected and removed based on reconstruction uncertainty, projection accuracy, and reprojection error, for which threshold values of 10, 3 and 0.3 were applied, respectively. After each selection and removal step, the internal and external camera parameters were optimised. The thresholds mentioned are target values. Depending on how many tiepoints would actually be removed, the threshold values were sometimes relaxed [37].

3.

The GCPs were imported into the project and subsequently used (1) to georeference the point cloud (or, in Phantom 4 surveys, to improve the georeferencing based on camera locations and orientations) and (2) to optimise the internal camera parameters. Columns H–K in NWNKern_Surveys_ErrorStatistics.ods provide statistics of the $XYZ$ residuals of the GCPs as the X, Y, Z and 3D (`total’) root mean square (rms) errors. The maximum 3D rms error for complete surveys was 0.05 m (2015-04-21). For most Phantom 4 surveys, the total rms error was less than 0.02 m. The error statistics for the first five EasyStar I surveys (2014-04-10 – 2016-04-01) provided here are lower than in [8] due to improvements in the selection and removal of low-quality tiepoints (step 2).

4.

A dense point cloud was calculated using medium-quality and aggressive filtering settings. Points with a confidence of less than 3 were removed, as well as off-terrain points due to cars, wooden posts, information signs, garbage bins, animals, or people. No attempt was made to remove off-terrain points due to vegetation; this is discussed in Section 4. Then, a polygonal mesh model was calculated, which was finally used to generate an orthomosaic with ground-pixel resolution ($<0.05$  m) using the mosaic blend mode (default).

The orthomosaics provided in the repository were exported from the Agisoft project with 0.05 and 1.0 m resolution. Each dense point cloud was exported in LAZ format and subsequently processed into a DEM using MATLAB^® [38] (https://nl.mathworks.com/). The LAZ file was read into MATLAB^® using the lasFileReader function of the Lidar toolbox. Elevation Z at each location was taken as the average Z of all points within a search radius equal to the grid resolution (here, 1 m). These points were found using a $k-d$ tree search method, implemented in the rangesearch function of the Statistics and Machine Learning toolbox. Regions without Z data that were entirely surrounded by grid points with data were subsequently filled with Laplace interpolation, implemented in the TopoToolbox (v. 2) [39].

3.2. Lidar from an Uncrewed Aerial Vehicle (UAVLidar)

The company Shore Monitoring & Research performed lidar surveys on 2018-09-17 and 2019-02-22 using a DJI Matrice 600 Pro UAV equipped with an AL3-32 mobile laser scanner from Phoenix LiDAR Systems. For each survey, the collected point cloud was filtered to retain only ground points and provided as a LAZ file, with horizontal coordinates in Amersfoort/RD New (EPSG:28992) and elevation with respect to NAP. The LAZ file was processed into a DEM as described in the UAVSfM workflow (Section 3.1). The UAV was also equipped with an RGB camera. The aerial images collected were processed by Shore Monitoring & Research with AgiSoft PhotoScan Professional to an orthomosaic with a resolution of 0.02 and 0.10 m. These orthomosaics were resampled at a resolution of 0.05 and 1.0 m using the gdal software package (https://gdal.org/) for inclusion in the data set presented in this paper.

3.3. Airborne laser scanning (ALS)

All ALS point clouds except 2015-02-13 were collected as part of the annual Coastal Lidar campaign of the Dutch government organisation Rijkswaterstaat e.g., [40,41,42] and are publicly available in LAZ format from https://downloads.rijkswaterstaatdata.nl/. The 2015-02-13 survey was commissioned by PWN, the nature manager of the study region, and the acquired data were provided as ASC Arc/Info ASCII grids with 0.5 m resolution. For the purpose of the present work, these files were merged and converted to LAZ format. The PWN-commissioned survey also resulted in an orthomosaic of the NWNKern region. ALS surveys were generally conducted in January, February or March, prior to the start of the growing season. The point clouds were corrected for vegetation and additional off-terrain objects by the commercial contractors who performed the data collection. The horizontal coordinates are in the Amersfoort/RD New system (EPSG:28992) and the elevation is with respect to NAP. Data quality documents illustrate that elevation data meet the requirements of Rijkswaterstaat, having a bias less than 0.05 m and a standard deviation less than 0.1 m. Because the focus of the annual campaign until 2022 was on the beach and the foredune, the most landward part of the depositional lobes was missing from the point clouds. Inquiries with Rijkswaterstaat revealed that Lidar data points have been collected landward of the foredune since 2016, but have not been processed as standard into the publicly available point clouds. For the present data set, these additional points were made available and merged with the point clouds of the beach and the foredune. Rijkswaterstaat indicated that the quality of these additional data may deviate from the requirements. The 2008 and 2014 ALS point clouds do contain high-quality data points landward of the foredune, because these surveys were carried out simultaneously with the second and third national-scale AHN campaigns (https://www.pdok.nl/introductie/-/article/actueel-hoogtebestand-nederland-ahn). From 2023, Rijkswaterstaat expanded the area for Lidar data collection and processing in the NWNKern region far landward of the depositional lobes. Each ALS point cloud was processed into a DEM in the same way as described for the UAVSfM point clouds.

4. User notes

The DEMs in the present data set are based on two remote sensing techniques: laser scanning (UAVLidar and ALS) and Structure-from-Motion photogrammetry (UAVSfM). Laser pulses can, at least partially, penetrate through the canopy, and the resulting point cloud comprises vegetation and ground points. With a suitable filter method, vegetation points can be selected and removed (as was done here), implying that the resulting DEM is a Digital Terrain Model (DTM), also in vegetated areas. If the vegetation is very dense, ground points cannot be found; this might be the case on the marram-grass-covered foredune. Structure-from-Motion does not result in ground points beneath the canopy, and hence a DEM based on UAVSfM is a Digital Surface Model (DSM), including the elevation of the canopy e.g., [43,44,45]. The difference between a DTM and a DSM is illustrated in Figure 5a, which shows the DoD between the ALS DEM of 2023-02-08 and the UAVSfM DEM of 2023-02-13. Because the two surveys were conducted only five days apart (with mainly low wind conditions), we can safely assume that the calculated elevation change actually represents the height of the canopy. The DoD shows a clear cross-shore zonation. A canopy height of up to 0.6 m is mostly found on the seaward side of the foredune south and north of the notches, on the foredune between the notches, and along the edges of some depositional lobes. At these locations, the vegetation is dominated by marram grass (habitat H2120, White dunes) [27]. Further inland, canopy heights to over 2 m can be seen, most notably on the seaward side of the foredune south of N1 and north of N5. These higher canopy elevations correspond to dark green vegetation patches in the orthomosaic of 2023-09-25 (that is, taken in the growing season of the same year, when the vegetation is more clearly visible; Figure 5b). These dark green patches are sea buckthorn and wild privet shrubs (habitat H2160) [27]. Figure 5 also shows that substantial parts of vegetated inland dunes are not associated with notable differences between the two DEMs, for example, to the east and south of the depositional lobes of N1 and N2. Here, the dominant habitat is H2130 (Grey dunes) with centimetre-high lichens and mosses, and other short plants [27].

The presence of DTMs (UAVLidar and ALS) and DSMs (UAVSfM) in the data set raises the question to what extent the volume computations in Figure 2c were affected by vegetation-induced bias in the DSMs. Therefore, an analysis of covariance (ANOCOVA) was performed to evaluate whether there is a significant difference between the regression slopes m for $V_d$, $V_e$, and $V_t$ based on DTMs and DSMs. Table 1 shows that the slopes for the three variables are different, but none of these differences is significant at the 99% confidence level $(p>0.01)$. Only the slopes for $V_d$ are different at the 95% confidence level $(p<0.05)$. In conclusion, it can be stated that the presence of canopy height in the DSMs does not have a major influence on the multi-year volume computations. Most likely, the true elevation changes completely outweigh the vegetation-induced bias, both in magnitude (between -9 and +10 m in 2023) and in spatial extent (Figure 2b). Furthermore, the highest vegetation (sea buckthorn and wild privet) is located mostly outside the polygon used in the volume computations, compare Figure 5a to Figure 2b.

Although the previous paragraphs considered vegetation-induced bias in DSMs as a potential limiting factor for their use, a canopy height map, as presented in Figure 5a, is also an important derivative of the current data set and provides opportunities for new research. A map with similar spatial extent and resolution will be virtually impossible to make on the basis of vegetation ground surveys. Potentially, a canopy height map can be used to train and test a machine learning algorithm, such as that presented in [45], to predict and remove vegetation-induced bias in a DSM; in other words, to correct it to a DTM. Another use may be to aid in the prediction of habitat maps from orthomosaics using deep learning. [27] used a Convolutional Neural Network (CNN) with U-net architecture to predict a habitat map from the orthomosaic of 2022-03-22 (with 75% accuracy), with training and testing data derived from vegetation surveys in the field and subsequent visual interpretation of the orthomosaic. Other studies [46,47] have shown that including elevation in the machine learning algorithm can improve the accuracy of prediction. [27] did not obtain any improvement by including elevation, possibly because the habitats in the study region are not limited to unique elevation ranges. Instead, she suggested including cross-shore distance due to the zonation of coastal dune habitats with the distance to the beach [2,48,49,50]. Figure 5a suggests that including the height of the canopy in the machine learning algorithm can also improve the accuracy of habitat predictions, as habitats have different canopy heights.

Finally, we note that monitoring of the foredune notches and their associated depositional lobes in the NWNKern region will continue in the years to come. We are particularly interested in (i) seeing for how long the notches remain corridors for the inland transport of windblown sand and how this depends on the characteristics of the notches, (ii) quantifying spatio-temporal changes in habitats, and (iii) understanding when and under what environmental conditions vegetation will eventually colonise the notches and depositional lobes, as expected from the biogeomorphic succession of natural notches (trough blowouts, [51,52]). We intend to include newly acquired data in the Zenodo repository by DOI versioning on an annual basis. The DOI provided here (https://doi.org/10.5281/zenodo.7010095) will always resolve to the latest version.