Structure-from-Motion-Derived Digital Surface Models from Historical Aerial Photographs: A New 3D Application for Coastal Dune Monitoring

Recent advances in structure-from-motion (SfM) techniques have proliferated the use of unmanned aerial vehicles (UAVs) in the monitoring of coastal landform changes, particularly when applied in the reconstruction of 3D surface models from historical aerial photographs. Here, we explore a number of depth map filtering and point cloud cleaning methods using the commercial software Agisoft Metashape Pro to determine the optimal methodology to build reliable digital surface models (DSMs). Twelve different aerial photography-derived DSMs are validated and compared against light detection and ranging (LiDAR)and UAV-derived DSMs of a vegetated coastal dune system that has undergone several decades of coastline retreat. The different studied methods showed an average vertical error (root mean square error, RMSE) of approximately 1 m, with the best method resulting in an error value of 0.93 m. In our case, the best method resulted from the removal of confidence values in the range of 0–3 from the dense point cloud (DPC), with no filter applied to the depth maps. Differences among the methods examined were associated with the reconstruction of the dune slipface. The application of the modern SfM methodology to the analysis of historical aerial (vertical) photography is a novel (and reliable) new approach that can be used to better quantify coastal dune volume changes. DSMs derived from suitable historical aerial photographs, therefore, represent dependable sources of 3D data that can be used to better analyse long-term geomorphic changes in coastal dune areas that have undergone retreat.


Introduction
The reliability of three-dimensional data is a crucial component in the physical analysis of exogenous geomorphic processes. These types of data can now be gathered much more accurately thanks to recent improvements in and accessibility of modern light detection and ranging (LiDAR), differential global navigation satellite system (DGNSS), and unmanned aerial vehicle (UAV) systems [1]. Despite the wide range of available technologies, the production of high-quality digital elevation models (DEMs) still requires significant investments in personnel time, hardware, and software [2]. In many cases, the high frequency and dense spatial accuracy provided by modern topographic acquisition systems such as UAVs is still not able to provide historical (multi-decadal) datasets due to their relatively recent (previous 10 years) widespread availability. There is a growing need for reliable long-term 3D data across all geomorphic disciplines, especially in studies on rapidly changing dynamic coastal environments [3]. Reliable long-term 3D data would help improve geomorphic predictions such as shoreline variation models [4,5], and would aerial photographs. Twelve methods were used to build and refine the dense point cloud (DPC). Four different DPCs were built using the four options available in the software to filter the depth maps: (i) disabled, (ii) mild, (iii) moderate, and (iv) aggressive. Each DPC was then filtered to remove poor confidence values calculated by the software during the DPC reconstruction. In particular, we chose to clean the DPC in three ways: (i) without any further DPC cleaning (keeping all confidence values calculated by the software); (ii) by removing points from the DPC with confidence values between 0 and 2; (iii) by removing points from the DPC with confidence values between 0 and 3. Therefore, a total of twelve methods were adopted, resulting in twelve final DSMs.
The study site was chosen to compare archival aerial-derived DSMs to recent topographic data acquired from more accurate LiDAR and UAV technologies. Our study, therefore, assesses if the SfM analysis applied to historical aerial photographs can provide a robust method for building reliable, long-term (decadal) 3D models of coastal dune volume loss on retreating coasts.

Study Area
The study site is located in the Murlough beach and dune complex of Dundrum Bay (County Down, Northern Ireland, UK) ( Figure 1). A wide sandy beach covers a large part of the bay, dominated by a unique geomorphic beach system known as multiple inter-tidal bars (MITB) [46], as well as an extensive ebb tidal delta feature, which effectively divides the southwest part of Newcastle and Murlough beaches from the eastern site of Ballykinler. A vegetated dune field dating back to the mid-late Holocene [47,48] extends from Newcastle to the inlet channel and from the latter to Ballykinler. A steep gravel ridge, representing the high tide mark [49], lies along the dune toe. The tidal environment is high mesotidal, with a mean spring tidal range of 3.85 m. Short-period wind waves with a mean significant wave height of 1 m are typical [49,50]. Modal winds are common from south and southwest directions, while storm winds are also common from the southeast quadrant [50,51]. The dune system of the Murlough Special Area of Conservation (SAC) has been managed by the National Trust since the 1960s.
Two areas of interest (AOIs) were chosen to perform comparisons between aerialphotograph-derived DSMs and a DSM derived from LiDAR and UAV surveys. A small AOI measuring approximately 220 m long and 100 m wide is located on a small portion of the coastal dune system (Figure 1). A larger AOI, approximately extending from the Royal County Down Golf Course, Newcastle, to the tidal inlet in the east (Figure 1), was chosen a posteriori based on the aerial photo overlap resulted from SfM process.

Historical Aerial Photographic Dataset
Six aerial photographs dated from 1963 were provided by the Ordnance Survey of Northern Ireland (OSNI). The aerial photographs (APs) were black and white with a 1:10,000 scale. Each photograph was scanned at 1200 dpi to obtain a high pixel resolution and to simplify ground control point (GCP) positioning and the surface matching during the SfM process. According to [27], the optimal scanning resolution to apply the SfM process on historical aerial photos is between 800 and 1600 dpi. The typical black border of archive aerial photographs was removed using a batch process in Adobe Photoshop [36]. The batch process was launched using the same selection area and taking care to select the largest area possible to maximise the overlap between images (Figure 2a). The final size achieved for each image (.TIFF format) was 9184 × 10352 pixels, with an average weight of 278 Mb per image.

Historical Aerial Photographic Dataset
Six aerial photographs dated from 1963 were provided by the Ordnance Survey of Northern Ireland (OSNI). The aerial photographs (APs) were black and white with a 1:10,000 scale. Each photograph was scanned at 1200 dpi to obtain a high pixel resolution and to simplify ground control point (GCP) positioning and the surface matching during the SfM process. According to [27], the optimal scanning resolution to apply the SfM process on historical aerial photos is between 800 and 1600 dpi. The typical black border of archive aerial photographs was removed using a batch process in Adobe Photoshop [36]. The batch process was launched using the same selection area and taking care to select the largest area possible to maximise the overlap between images (Figure 2a). The final size achieved for each image (.TIFF format) was 9184 × 10352 pixels, with an average weight of 278 Mb per image. Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 24 5 Figure 2. Preparation workflow of aerial photographs (APs) from 1963 prior to the structure-from-motion process (a) and main steps (b) of the SfM process adopted for the APs from 1963 to produce the 12 digital surface models.

SfM Processing of Historical Aerial Photographs
The SfM process was performed by means of the commercial software Agisoft Metashape Pro (version 1.6.3; 64 bit). Prior to beginning the SfM process, a quality check performed by the software was made of each image, resulting in none having a coefficient quality lower than 0.63. Software specifications recommend to only use images with a Figure 2. Preparation workflow of aerial photographs (APs) from 1963 prior to the structure-from-motion process (a) and main steps (b) of the SfM process adopted for the APs from 1963 to produce the 12 digital surface models.

SfM Processing of Historical Aerial Photographs
The SfM process was performed by means of the commercial software Agisoft Metashape Pro (version 1.6.3; 64 bit). Prior to beginning the SfM process, a quality check performed by the software was made of each image, resulting in none having a coefficient quality lower than 0.63. Software specifications recommend to only use images with a coefficient quality higher than 0.5, since poor-quality input images can adversely affect alignment results, and therefore can also affect the following steps in the photogrammetric process [52]. Image alignment was performed by selecting "highest" as the accuracy in order to generate a denser sparse point cloud (SPC, a.k.a. tie points) at the end of this step (Figure 2b).
GCPs were identified by comparing the 2014 LiDAR-derived orthophotos and the 1963 aerial photographs by identifying locations where changes were minimal over the five decades. GCP coordinates (X,Y,Z) were retrieved from the 2014 LiDAR-derived DSM by carefully selecting points where the Z value was not influenced by nearby walls, buildings, or tall vegetation (Figure 3a-i). Highly dependent on local landscape changes (natural and anthropic) in the area during the 50 years in between, this was the most laborious but important step, since poor GCP recognition and distribution can adversely influence the final DSM output [34,35]. For the entire dataset, a total of 28 GCPs were identified, with a good spatial spread across the area represented by the six aerial photographs (Figure 4a).    Since each GCP must be visible in a minimum of three photographs (i.e., three projections for each GCP) to make it reliable [36], 17 of the original 28 GCPs were used in the photogrammetric process (Figure 4a). This resulted in a reduced reliable area produced by the SfM process approximately corresponding to the large AOI. Considering GCPs identified from the 1963 APs (spatial resolution 0.47 m) and from the LiDAR-derived DSM and orthophotos (spatial resolutions of 0.5 and 0.05 m, respectively), the XY accuracy of GCPs assigned in the SfM process was 0.5 m, whereas the Z accuracy was 0.069 m, as reported in the LiDAR metadata by the original contractors. After removing poorly projected GCPs, images were optimised and the SPC was manually cleaned, removing points located outside the whole GCP configuration (Figure 2b). GCPs returned mean errors on their locations of 2.3, 2.5, and 0.06 m for X, Y, and Z, respectively.
At this stage of the SfM process, various methods were used to build and refine the dense point cloud (DPC). Four different DPCs were built using the four options available in the software to filter the depth maps: (i) disabled, (ii) mild, (iii) moderate, and (iv) aggressive. For each method, we kept the DPC quality at "ultra high" to take advantage of the entire pixel number of APs. We also let the software estimate the confidence values of the four DPCs ( Figure 2b).
Each DPC was then filtered to remove poor confidence values calculated by the software during the DPC reconstruction. In particular, we chose to clean the DPC in three ways: (i) without any further DPC cleaning (so keeping all the confidence values calculated by the software); (ii) by removing points from the DPC with confidence values between 0 and 2; (iii) by removing points from the DPC with confidence values between 0 and 3. Therefore, a total of twelve methods were adopted, resulting in twelve final DSMs (Figure 2b).

LiDAR Dataset
An airborne LiDAR dataset was acquired on the 4th of March 2014 for the coastal area that extends from the Newcastle town center to Ballykinler ( Figure 1). The dataset was provided by the National Trust of the Murlough National Nature Reserve and the Ministry of Defence of the Ballykinler army camp. The LiDAR sensor acquired multiple returns with a point density greater than 30 ppm. Horizontal and vertical accuracy (root mean square error) values of 0.100 and 0.069 m, respectively, were reported in the metadata. The dataset was organised into tiles measuring 250 × 250 m provided by the original contractors, with a spatial resolution of 0.5 m. The flight was undertaken during high tide, meaning the topographic dataset is only reliable up to the gravel ridge, which lies at the dune toe.

UAV Dataset
A UAV survey was undertaken on the 2nd of October 2019 over the small AOI in Murlough beach ( Figure 1) in ideal weather conditions, with wind speeds of less than 5 m/s. A flight plan was performed using Emotion flight planning software and an eBee fixed-wing UAV equipped with an RGB camera with 18.2 MP resolution at a flight altitude of 120 m and with 70% and 80% forward and lateral overlap, respectively. These settings enabled a ground image resolution of 0.030 m. Before the flight, nine checkerboard targets (GCPs, Figure 3l,m) were placed in order to include the foredune and the facing subaerial beach in the survey area ( Figure 4b). Each GCP was measured using a DGNSS device (Trimble R10), with mean accuracy values of +/−0.016 m (horizontal) and 0.017 m (vertical). A total of 210 images were taken during a flight time of 17 min (divided into two subsequent flights), covering an area of 0.224 km.

SfM Processing of UAV Images
As was done for the APs, the SfM processing of UAV images was performed using Agisoft Metashape Pro (version 1.6.3; 64 bit) according to the following subsequent steps: (i) image alignment at the highest accuracy level; (ii) GCP positioning; (iii) manual cleaning of the SPC in the area outside the GCP spatial configuration; (iv) optimisation of camera alignment; (v) generation of the DPC with "ultra high" accuracy and the filter "disabled" for the depth maps; (vi) DPC cleaning by removing points with confidence values between 0 and 3, namely points where the software was unsure about the position; (vii) DSM and orthophoto generation from the dense point cloud. The final UAV-derived DSM reached a resolution of 0.030 m. GCPs returned mean errors on their locations of 0.026, 0.033, and 0.025 m for X, Y, and Z, respectively. A summary of the main characteristics of each SfM product is presented in Table 1.

1963 DSM Accuracy Assessment and Validation
The quality of the DSMs produced by the twelve methods adopted for the processing of the 1963 APs was investigated by comparing Z values derived from the 2014 LiDARderived DSM, which were more reliable. The accuracy assessment was only possible for the large AOI (Figure 1), since areas with minimal changes such as roads, parking areas, and footpaths had to be identified between 1963 and 2014. In total, 512 independent points (a.k.a. check points) where topographic changes between 1963 and 2014 were supposed to be minimal were used in the validation. The root mean square error (RMSE), mean absolute error (MAE), standard deviation error (SDE), and the coefficient of determination (R 2 ) were estimated on the same point population to check the quality of the twelve 1963 DSMs.

Geomorphic Changes and Volumetric Variation
Using the best DSM from the accuracy assessment, quantitative analyses were performed within the small AOI to estimate volume variations from 1963 to 2019. DEMs of differences (DoDs) and shorelines (i.e., stable vegetation lines) between each period were produced. Shorelines were digitised on the orthophoto of each dataset. DoDs were cut accordingly to the most recent shoreline for each period of comparison in order to avoid the influence of the dune vegetation as much as possible, which can complicate surface expression in the DSM derivation process. For each shoreline, an error was calculated as the root sum of the squares of the individual uncertainties associated with each dataset where the shoreline was digitised. Volume calculations were performed from the lowest elevation among the three DSMs (1.12 m above Ordnance Datum Belfast). Areas that experienced elevation changes small enough to be within the limit of detection (LoD) of each DSM comparison were removed from volumetric calculations. The LoD is a commonly used threshold in geomorphology when DoDs are analysed [44,53,54]. The LoD was considered equal to the standard deviation of the RMSE of each validated dataset (UAV-, LiDAR-, and AP-derived DSMs). Volume uncertainty was, therefore, estimated by multiplying the LoD for the area that experienced significant changes (i.e., above or below LoD).

Accuracy Assessment of Aerial-Photograph-Derived DSMs
Among the twelve DSMs produced by adopting different depth map filtering and DPC cleaning modes (see Section 3.1.1), the best results were shown by the DSM produced with a "disabled" filter for the depth maps and with DPC cleaning performed removing points with confidence values ranging from 0 to 3. This method gave the smallest error values in terms of RMSE, MAE, and SDE and the strongest relationship (R 2 ) between the 2014 LiDAR elevations and the elevations extracted from the AP-derived DSM ( Table 2). The twelve DSMs showed RMSE values ranging between 0.93 and 1.12 m; with MAE values showing a larger variability, ranging from 0.86 to 1.25 m; and with SDE values that were restricted between 0.042 and 0.050 m. The variability was also very low among the R 2 values, which were between 0.991 and 0.994 (Table 2 and Figure 5). This also means that the elevation residuals between the LiDAR and the AP-derived DSMs were quite consistent (following a linear trend) over the entire study site (i.e., 512 check points) and for each elevation range ( Figure 5). The frequency distributions of elevation residual were normal, with a slightly positive skew in the majority of cases. The majority of residuals were included between +1 and −1 m ( Figure 5).
The elevation residuals from the best DSM resulting from the SfM process of APs showed a spatial pattern-the lowest values were concentrated where the best overlap between APs was possible, namely where all six aerial photographs overlap, in the central part of the large AOI ( Figure 6). This was also the area that included the small AOI where the UAV survey was performed. Outside this area of maximum overlap, the largest differences were found in the dune area-cold colours (positive values) were more common moving towards the SW near the golf course, while warm colours (negative values) were more frequent moving towards NE near the tidal inlet and the sand spit ( Figure 6). Cold colours (positive values) mean that the LiDAR elevations are higher than the elevations from the 1963 APs derived from DSM, and vice versa.
The tie points variance in the SPC from the image alignment and camera optimisation (common to all twelve methods; Figure 2b) was good (cold colors in Figure 7a) in both the large AOI and in the small AOI. The dense point cloud confidence values varied between the three DPC cleaning modes. With no re-motion for poor DPC confidence values, the points used by the software to build the DSM covered a wider area, however with high uncertainty regarding their identification (warm colors in Figure 7b). With confidence values ranging between 0 and 2 removed from the DPC, the software then relied on more solid points, despite a restriction in the overall area (Figure 7c). With confidence values ranging between 0 and 3 removed, the area containing reliable 3D points became smaller, corresponding to the area of the maximum overlap between APs, however with more certainty in their identification and location (Figure 7d). A further qualitative check on how the twelve methods performed the DSM reconstruction from 1963 APs is given in Figure 8. A reference transect was chosen within the small AOI and each 1963 DSM was compared to the recent and more reliable topography of the 2014 LiDAR and 2019 UAV-derived DSMs. The dune topography was well-outlined until the foredune crest using each of the 12 methods; a greater variation is visible in the stoss dune faces, with decreasing noise from the methods with no DPC cleaning to the methods where the DPC cleaning was performed, with confidence values between 0 and 3 being removed (Figure 8). Despite the variation between the twelve methods used to reconstruct the topography of the dune slipface, a clear subsequent erosion of the 1963 foredune was observed.  The elevation residuals from the best DSM resulting from the SfM process of APs showed a spatial pattern-the lowest values were concentrated where the best overlap between APs was possible, namely where all six aerial photographs overlap, in the central part of the large AOI ( Figure 6). This was also the area that included the small AOI where the UAV survey was performed. Outside this area of maximum overlap, the largest differences were found in the dune area-cold colours (positive values) were more common moving towards the SW near the golf course, while warm colours (negative values) were more frequent moving towards NE near the tidal inlet and the sand spit ( Figure 6). Cold colours (positive values) mean that the LiDAR elevations are higher than the elevations from the 1963 APs derived from DSM, and vice versa.

Geomorphic Changes and Volumetric Variation
Further details of the 1963 foredune erosion in the small AOI can be seen in the DoDs and the shorelines shown in Figure 9. From 1963 to 2014, the shorelines (i.e., stable vegetation lines) experienced retreat distances of 5 ± 2 m (min) to 21 ± 2 m (max). These values were also representative of locations where the minimum and maximum elevation differences occurred (Figure 9a). From 1963 to 2014, total amounts of 11,139 ± 3429 m 3 were eroded over a total length of 220 m. From 2014 to 2019, the shorelines remained relatively stable, with a maximum retreat distance of 7 ± 1 m and maximum seaward advance distance of 5 ± 1 m. General accretion volumes of 1921 ± 323 m 3 affected the gravel ridge, with some erosion spots corresponding to a recent line of beach cusps detected during the 2019 UAV survey (Figure 9b). The overall period from 1963 to 2019 showed minimum retreat distances of 6 ± 2 m up to maximum values of 25 ± 2 m (Figure 9c), with eroded volumes of 9423 ± 3360 m 3 . To give an extra qualitative visualisation of the shoreline retreat over the last 56 years, the three shorelines were superimposed on the best 1963 DSM that resulted from the SfM analyses ( Figure 10).
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 24 13 Figure 6. Elevation difference between the digital surface model (DSM) derived from the light detection and ranging (LiDAR) survey from 2014 and the aerial photograph (AP) derived from the DSM from 1963 using the best structure-frommotion method resulting from the accuracy assessment. Cold colors (positive values) mean that the LiDAR-derived DSM is higher than the AP-derived DSM, and vice versa. The background 2014 orthophoto is derived from the LiDAR survey. The graticule is in TM75/Irish Grid (transverse mercator) projection.
The tie points variance in the SPC from the image alignment and camera optimisation (common to all twelve methods; Figure 2b) was good (cold colors in Figure 7a) in both the large AOI and in the small AOI. The dense point cloud confidence values varied between the three DPC cleaning modes. With no re-motion for poor DPC confidence values, the points used by the software to build the DSM covered a wider area, however with high uncertainty regarding their identification (warm colors in Figure 7b). With confidence values ranging between 0 and 2 removed from the DPC, the software then relied on more solid points, despite a restriction in the overall area (Figure 7c). With confidence values ranging between 0 and 3 removed, the area containing reliable 3D points became smaller, corresponding to the area of the maximum overlap between APs, however with more certainty in their identification and location (Figure 7d).

Discussion and Conclusions
All twelve adopted methods for depth map filtering and DPC cleaning gave an average error in Z of approximately 1 m. The more DPC cleaning performed and the less aggressive depth map filters adopted, the higher the accuracy of the topographic reconstruction. In our case, performing a good DPC cleaning process (0-3 confidence values removed from the DPC) was far more important than the filter adopted for the depth map reconstruction. Despite small differences in the computed Z errors were found among the adopted methods (RMSE values: 0.93-1.17 m), the largest disadvantage emerged from the removal of points from the DPC with confidence values ranging between 0 and 3, which resulted in a consistent restriction of the modeled area (Figure 7d). This aspect is highly important when approaching SfM techniques used on historical aerial photos, since the AOI cannot be entirely decided a priori compared to UAV surveys, where an ad hoc image capture must be planned according to the AOI [55]. Therefore, the DPC cleaning always has to be performed according to the objectives of the study, the quality of the aerial photographs, and the characteristics of the study site. The best method that emerged from our analyses cannot be considered valid in every topographic setting-robust validation of the SfM outputs always has to be performed depending on the objectives to be achieved and the environmental constraints of the study site. As argued by [35], the adoption of self-calibration routines provided by commercial software does not always guarantee good results, because such datasets are very dependent on local variables that cannot be extrapolated to other areas.
Despite all twelve methods demonstrating good performance in the reconstruction of the inner part of the large AOI (Table 2 and Figure 5) and recreating a similar topography up to the foredune crest, it is not clear why they differed so much in the slipface of the foredune (Figure 8). The small AOI was present on each of the six APs, therefore poor overlap between images was not a reason. Nevertheless, poor topographic reconstruction is also possible with good image overlap if the photographs are taken from similar viewpoints and with greater distances between the ground and cameras [30]. On the other hand, differences in camera orientation between images can influence the illumination of imaged areas and consequently alter the sharpness of homogeneous surfaces, thereby affecting the number and quality of identified tie points [1,[56][57][58]. Having no information on flight conditions for the aerial survey, no precise considerations were carried out on the 1963 APs. Nevertheless, significant changes to camera orientation axes, even though estimated by the software, can be largely excluded ( Table 3). The only exception was the yaw axis, which revealed a 180 • direction change performed by the aircraft to take the second strip of photos ( Table 3), but that should not have changed the illumination between the two strips. Variations in camera elevation can also adversely influence the image matching [1,59]. In our case, the maximum estimated difference between the images was 41 ± 25 m (Table 3), which could have disturbed the matching of homogeneous surfaces in the foredune slipface. Furthermore, that area, being very steep and narrow if compared to its rear, could have disturbed the image matching and could consequently have been interpretated by each depth map filter in different ways (Figure 8). One study [60] using APs with similar survey acquisition elevations (3000-4000 m) and similar image resolutions (0.42-0.49 m) found absolute Z errors of 0.6-0.8 m, even though more recent APs (from 1980s to 2000s) were used. The same authors using another set of APs found an error of 1.2 m with coeval photos (1962), a similar scale (1:10,000), and more than double resolution (0.20 m) if compared to that used in our study. Having the actual camera specifications and flight characteristics could slightly improve the error, as demonstrated by [60], even though it is very rare that APs come with all of these specifications [36].
The main cause of the different reconstruction results for the foredune slipface was probably due to the fact that the area is outside the GCP configuration. Indeed, many authors have argued that for SfM-derived DSMs, better reliability occurs in the area delimited by GCPs [16,59,61,62]. Despite the GCPs being well-spread across the large AOI ( Figure 4a) and one GCP being in the small AOI (Figure 8), the lack of GCPs on the seaward side makes the area outside the entire GCP configuration less reliable in the seaward area. For this reason, coastal areas sometimes containing artificial structures such as groynes or jetties could help extend GCP configurations further seaward and provide a more reliable topographic reconstruction of the coastal fringe. Similar issues are well known in the reconstruction of former glacial topographies, which are often unlikely to be fixed over long periods of time due to the degrading permafrost and glacial recession [63]. In these environments, ground control points are limited to stable (non-glacial) locations within the image scene, which may leave large areas of the surface model unconstrained and subject to greater errors [63]. The DSMs derived from the SfM process of aerial photographs are usually a partial coverage dataset [36,37], since good aerial coverage of the final DSM depends on a good number and spatially uniform distribution of GCPs and a good overlap between photographs. Furthermore, aerial photographs typically come with a black or white border with indentations at the four corners, which reduce the overlap area between photos and the chance to place extra GCPs in this area. As stated by [25], image overlap or configuration and image texture are crucial drivers in the tie point generation and in the bundle adjustment phases. The level of uncertainty provided by the AP-derived DSMs is not comparable to the LoD achievable today with modern spatial acquisition technologies such as LiDAR, DGNSS, UAVs, or even recent satellite data. Nevertheless, aerial photographs seem to be the most common documents available for identifying 3D geomorphic features in the period ranging from the 1930s to 1990s (depending on local photos available) [45,64,65]; therefore, a vertical uncertainty of barely 1 m is acceptable. One study [36] found Z error values (RMSE) varying between 0.50 to 0.91 m over two study sites and three different AP datasets. Another study [2] measured the RMSE in Z of 1.049 m, whereas [27] measured RMSE values in Z ranging between 2.11 and 1.98 m, analysing APs with resolutions of 800 and 1600 dpi over a complex topography. The uncertainty could be improved if a reliable and automatic method could be developed to remove the influence of vegetation from the DSM. When dealing with complex topography such as the vegetated dune area of this study site or glacial topographies in steep mountains, vertical differences may increase, mainly due to the influence of poorly constrained topographies because of a limited viewshed, shading, and steep terrain [63]. The relationship between slope and topographic error is indeed well known [63], and it has been found that in very steep terrain, horizontal errors are closely coupled to the vertical ones [66][67][68]. We chose to validate the AP-derived DSMs against other more accurate DSMs, knowing that the vegetation would have been included. Therefore, during the validation stage, a crucial aspect is to identify areas where minimal changes occur.
Another option would be to measure the check points with a DGNSS, even though in our case, given the large area of AOI, not every measuring site would have been physically reachable. If the LiDAR data are accurate enough (here the values were 0.100 and 0.069 m for X-Y and Z, respectively) to detect long-term geomorphic processes (e.g., dune erosion over several decades), the SfM approach for historical aerial photos can also be adopted elsewhere without the use of a DGNSS and with a reliable GCP configuration, as demonstrated by [34,36]. The option of validating AP-derived DSMs against LiDAR data is a common approach, especially when the AOIs are large [34,36] or where the field sites are not easily accessible [2,38]. In some cases, DEMs derived from SfM process of historical APs can also be validated against more recent and more accurate DEMs [26,29]. Furthermore, finding sufficient check points that are easily and precisely identifiable to be measured with a DGNSS can be problematic in actively changing landscapes [28]. As already highlighted by [36], SfM approaches using historical aerial photos can be used to approximately quantify volume variations on sandy coastlines, albeit with consideration of the LoD in any volume calculations [44,53]. Furthermore, photos taken at low tide are preferred if volume calculations of the supratidal and intertidal beach are performed.
Despite no information being available for 1963 APs in terms of flight conditions and camera characteristics, the twelve applied methods gave a very satisfactory vertical error of approximately 1 m, with a maximum RMSE of 0.93 m. In our case, the best combination of steps in the SfM process resulted from removing confidence values ranging between 0 and 3 from the DPC and disabling the filter for the depth map reconstruction. Some issues arose in the reconstruction of the slipface of the foredune, likely due to the paucity of GCPs on the seaward side of the AOIs. This represents a typical limitation in the use of SfM process for historical aerial photographs and possibly explains the lack of adoption of this method for coastal geomorphological studies. AP-derived DSMs can be used to quantify approximate volumetric changes if reliable estimations of the Z error and LoD of the DoDs are provided. Despite these limitations, the SfM process of historical aerial photographs is promising and represents an opportunity to build long-term 3D datasets of coastal dune volumes to add to other existing datasets at sites. Reliable long-term 3D datasets could be used to better refine the numerical modeling of shoreline variation, dune erosion, or even hydrodynamic models of the subaerial beach if the GCP configuration can be improved in seaward areas. Funding: This work is part of the MarPAMM project, funded by the European Union's Interreg V-A Programme (https://www.mpa-management.eu). This work is a contribution to the Natural Environment Research Council project NE/F019483/1.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.