Unveiling buried aeolian landscapes: reconstructing a late Holocene dune environment using 3D ground‐penetrating radar

Across the UK, sandy beaches and dunes protect coastal infrastructure from waves and extreme water levels during large‐scale storms, while providing important habitats and recreational opportunities. Understanding their long‐term evolution is vital in managing their condition in a changing climate. Recently, ground‐penetrating radar (GPR) methods have grown in popularity in geomorphological applications, yielding centimetre‐scale resolution images of near‐surface stratigraphy and structure, thus allowing landscape evolution to be reconstructed. Additionally, abrupt changes in palaeo‐environments can be visualized in three dimensions. Although often complemented by core data, GPR allows interpretations to be extended into areas with minimal ground‐truth control. Nonetheless, GPR data interpretation can be non‐intuitive and ambiguous, and radargrams may not initially resemble the expected subsurface geometry. Interpretation can be made yet more onerous when handling the large 3D data volumes that are facilitated with modern GPR technology. Here we describe the development of novel semi‐automated GPR feature‐extraction tools, based on ‘edge detection’ and ‘thresholding’ methods, which detect regions of increased GPR reflectivity which can be applied to aid in the reconstruction of a range Quaternary landscapes. Since reflectivity can be related to lithological and/or pore fluid changes, the 3D architecture of the palaeo‐landscape can be reconstructed from the features extracted from a geophysical dataset. We present 500 MHz GPR data collected over a buried Holocene coastal dune system in North Wales, UK, now reclaimed for use as an airfield. Core data from the site, reaching a maximum depth 2 m, suggest rapid vertical changes from sand to silty‐organic units, and GPR profiles suggest similar lateral complexity. By applying thresholding methods to GPR depth slices, these lateral complexities are effectively and automatically mapped. Furthermore, automatic extraction of the local reflection power yields a strong correlation with the depth variation of organic content, suggesting it is a cause of reflectivity contrast. GPR‐interpolated analyses away from core control thus offer a powerful proxy for parameters derived from invasive core logging. The GPR data collected at Llanbedr airfield highlight a complex dune system to a depth of 2.8 m, probably deposited in several phases over ~700 years, similar to elsewhere in North Wales.


Introduction
Climate projections suggest that coastal geomorphology will change dramatically by 2100 due to sea-level rise and an increased frequency of intense storms (Russell et al., 2018). Coastal dunes are important landforms in that they provide protection from wind gusts and storm surges, reducing coastal erosion and providing habitats for unique ecosystems (García-Mora et al., 2001;Judge et al., 2003;Sigren et al., 2014). Although some aspects of coastal dunefields can rapidly adapt to climatic shifts, anthropogenic activities have reduced their resilience (Williams and Davies, 2001;Provoost et al., 2011); for example, coastal infrastructure has restricted the landward migration of dunes (Jones et al., 2013). Projections of sea-level rise and erosion suggest a net 36% loss (~36 800 ha) in the total dune area across the UK by 2060 (Dargie, 1993;Radley, 1994;Dargie, 1995;Saye and Pye, 2007), with implications for dunefield geomorphology, regional habitats and coastal infrastructure.
Most coastal dunefields in northwest Europe formed during the Little Ice Age (LIA), between the early 13th and mid-17th centuries (Pye and Neal, 1993;van Vliet-Lanoë et al., 2016;Kelley et al., 2018;Jackson et al., 2019). Climatic cooling during the LIA led to an increase in the frequency, magnitude and duration of strong winds and storms, causing widespread coastal erosion and inland migration of sand (Wilson et al., 2001(Wilson et al., , 2004Szkornik et al., 2008;Anthony et al., 2010). A subsequent lowering of the wave base (the maximum depth at which a wave passage causes substantial water motion) also contributed regionally to the increased landward movement of sand (Reed et al, 2009). Thereafter, a period of stable climate and sea level allowed dense dune vegetation to form, leading in turn to the preservation of many UK dunefields; prominent examples remain in northeast England, north Wales, southern England and northwest Scotland (e.g. Wilson, 2002;Bateman and Godby, 2004;Orme et al., 2015;Pye and Blott, 2017).
The study of coastal dune systems has historically applied a multifaceted approach, involving (i) examination of sedimentological sequences, exposed at the surface (Bernhardson and Alexanderson, 2018), (ii) auger, gouge and piston drilling, supplemented by excavation pits (Matthews and Seppälä, 2014;Nielsen et al., 2016), and (iii) re-examination of existing borehole data, historical maps and aerial photos (Forman et al., 2008;Lancaster and McCarley-Holder, 2013). Although this approach provides detailed information on the internal structure, composition and stratigraphy of coastal dune fields, it may neglect spatial subsurface variations which are critical to understanding dune evolution. The application of photogrammetry, LiDAR and geophysical methods e.g. groundpenetrating radar (GPR) has furthered the understanding of dune systems (Millington et al., 2008;Hugenholtz et al., 2012;Rockett et al., 2016;Fu et al., 2019). GPR has been widely applied to characterize palaeo-dunefields, providing visualizations of their internal sedimentary structure (e.g. Bristow et al., 2000aBristow et al., , 2010Pedersen and Clemmensen, 2005;Girardi and Davis, 2010;Trivedi et al., 2012;Bernhardson and Alexanderson, 2017;Zhao et al., 2018). The three-dimensional (3D) characterization of near-surface facies architecture, made possible with GPR surveying, can help reconstruct palaeo-wind patterns, depositional histories of Holocene dunefields, and the evolution of the dune system due to environmental forcing(s) (Huddart, 1992;Clemmensen et al., 2009). In this paper, GPR methods are incorporated into an analysis of a late Holocene dune system in north Wales, UK; in particular, this study aims to demonstrate the use of a novel, semi-automated GPR technique, as part of a broader approach to investigate the evolution of the field site.

GPR application in coastal dune research
For GPR methods to be applicable, a target feature must offer a contrast in electromagnetic properties with its host material. Although coastal dunes may seem homogenous, they comprise both sand and humic layers, the latter present in palaeosols formed by decomposed dune vegetation (Fig. 1). vegetation (Flemming and Hansom, 2011;Paranjoti, 2016). (b) Example of a buried dune palaeosol layered beneath dune sands. Typical core through a buried dune palaeosol, displaying the sharp transition between organicrich and quartz-rich dune sand (Nichol et al., 2004;Monaghan et al., 2013). [Color figure can be viewed at wileyonlinelibrary.com].
The potential for the organic material to retain moisture increases the water content of organic-rich layers, thereby reducing their dielectric permittivity. Thus, GPR can be readily applied to map palaeosols (Van Dam et al., 2002;Gusmeroli et al., 2015) and reconstruct the evolution of late Holocene dune fields.
A single traverse with a GPR instrument produces a two-dimensional (2D) cross-sectional slice through the ground, termed a 'radargram', but the efficiency of modern GPR systems means that it is commonplace to acquire multiple parallel survey lines to assess 3D variation of the target (Schennen et al., 2016;van Schoor et al., 2017;Allroggen et al., 2019). Certainly, the detailed architecture of a dune system is more efficiently revealed with a 3D GPR survey than by conventional coring methods (Żuk et al., 2017). A powerful advantage of the 3D dataset is that it can be displayed as horizontal 'timeslices' (Nuzzo et al., 2002;Trinks et al., 2010;Zhao et al., 2015) (Fig. 2), which are plan-view maps of the GPR response at a certain position beneath the ground surface. Such slices are particularly beneficial when delineating the spatial and vertical extent of features in a complex subsurface. With an estimate of the GPR propagation velocity through the subsurface, a 'timeslice' may be converted to a 'depth slice'.
The interpretation of features in GPR slices involves recognizing variations in reflection continuity, shape, amplitude, internal configuration and external geometries (Beres et al., 1999;González-Villanueva et al., 2011;Franke et al., 2015). Such variations can be categorized into socalled radar facies (Beres and Haeni, 1991;Jol and Smith, 1991). This analysis parallels the approach applied to interpret marine seismic datasets (e.g. Emery et al., 2019), but the recognition of unique patterns within GPR data can be more difficult (Moysey et al., 2006;Tronicke and Allroggen, 2015). Recently, inspired by practice in the medical field (e.g. X-ray, computed tomography and magentic resonance imaging scanning) (Hendriks et al., 2011;Christodoulou et al., 2012;Erdt et al., 2012;Mishra et al., 2015;Ali et al., 2017), image segmentation methods have been applied to GPR feature extraction, and include tools such as Thresholding and Edge Detection. Given an image file, such methods can highlight features of interest and reduce the volume of superfluous data. These methods have therefore been explored for GPR applications (Tronicke and Allroggen, 2015;Dou et al., 2017;Dinh et al., 2018), with the rationale of objectifying and automating facies analysis.
The application of a semi-automated algorithm to appraise a large GPR dataset would offer a useful development to GPR interpretation (e.g. Francke, 2007;Krause et al., 2007;Pasolli et al., 2009;Maas and Schmalzl, 2013;Tronicke and Allroggen, 2015;Uus et al., 2016). Although full automation may be desirable in theory, and practical in certain cases, such as the detection of diffraction hyperbolae (Shihab et al., 2002(Shihab et al., , 2003Neal, 2004;Nagrale and Bagde, 2013), the variability of GPR responses means that a single algorithm is unlikely to deliver a suitable interpretation in all cases, and hence user-supervision of the algorithm will typically remain valuable.

Field site and data acquisition
The field site is a flat grassed area within the grounds of Llanbedr Airfield, Gwynedd, North Wales, UK (British National Grid SH 57079 26608) (Fig. 3). Llanbedr Airfield was constructed in the 1940s, built on a low-lying coastal plain that protrudes into Cardigan Bay (Pye et al., 2007a). It is underlain by a combination of superficial material, including blown sand and tidal flat deposits (British Geological Society, 2020), though the extent of these units is unknown. The current airfield is 1 km north of the 10.6 km 2 transgressive parabolic Morfa Dyffryn dunefield, the northern limit of which is constrained by the southern perimeter of the airfield. The modern coastal dune morphology follows the classification outlined in Fig. 4 with a foredune ridge, followed by two primary parabolic dune ridges, dominated by extensive dune slack in the depressions formed behind these ridges, terminating at the southern perimeter of the airfield (Fig. 4). The Morfa Dyffryn dunefield lies parallel to the prevailing south-westerly winds, the rate of dune migration reaching up to 3.6 m per year, with a mean migration rate of 1.0 m, similar to other areas of the North Wales coast (Pye et al., 2007a(Pye et al., , 2007b(Pye et al., , 2014. GPR data were collected at Llanbedr Airfield in 2016 within a 110 × 52-m area, using a Sensors & Software PulseEKKO PRO system with 500 MHz monostatic antennas at a common-offset of 0.25 m. A GPR profile was recorded every 0.5 m along the grid, with a GPR trace recorded every 0.04 m along the profile. The time sampling interval in the survey was 0.2 ns. Coring was conducted in 2019 using a 1-m gouge corer to groundtruth the GPR data. Changes in groundwater conditions between the two survey periods may produce small mismatches between the horizons observed in the GPR and core data, although the accuracy of any such comparison is in any case governed by the accuracy of GPR time-to-depth conversions. Following the approach of Booth et al. (2010), a velocity of 0.075 m ns −1 was measured, with a precision of ± 8% for both the velocity analysis and the subsequent depth conversion. Considering how water content may change the GPR velocity and depth conversion, Robinson et al. (2013) suggest 0.08 m ns −1 for a dry sandy soil, representing a 7% increase in velocity and an equivalent decrease in depth. By contrast, saturated sandy soil could have a velocity as low as 0.055 m ns −1 , decreasing velocity and increasing depth by 20%. This range of velocity variation could, at most, add 0.5 m to our depth estimates. The potential range of water content change is therefore more significant than the underlying precision of the  Application of exponential gain that compensates for attenuation in the GPR data, boosting deeper reflection events and improving image quality.

Time-to-depth conversion
Conversion of two-way travel time (ns) to depth (m), using a wavelet velocity estimate.
Conversion Velocity: 0.075 m ns −1 inferred from hyperbolic diffraction fitting. depth conversion. However, such extreme variations are not anticipated at this site as the lithology has not changed over timehence we consider ± 8% to be a reasonable uncertainty for our depth conversion. Core locations targeted key reflections but also had to avoid the buried infrastructure visible within the shallow subsurface. The cores reached a maximum depth of 2 m, limited by the stiffness of sand, and were logged in the field following the methodology of Troels-Smith (1955).

GPR processing and image segmentation analysis
GPR image processing was conducted in Sandmeier Reflex-Win software (processing steps in Table 1). Conversion of GPR two-way travel time to depth used two velocity control methods: (i) matching the curvature of diffraction hyperbolae in the data, and (ii) a common-offset semblance approach (Booth et al., 2010) which provided a more objective velocity appraisal (see Supporting Information Data 1). Segmentation of the GPR data was conducted in Mathworks MATLAB 2017a, using conventional image segmentation tools coded into a semi-automated algorithm (Sharma et al., 2010;Dahab et al., 2012). This was applied to data in the timeslice domain, to optimize the appraisal of the spatial dimensions of targets within any one depth interval. The algorithm was tuned to maximize its sensitivity to the lithological changes and buried structures that appeared in the GPR data. The segmentation methodology consists of four phases given in Fig. 5 (see Supporting Information Data 2).
Phase 1: ground-penetrating radar data conversion Normalization of GPR data (Fig. 5a). Normalization is required to prevent datapoints with anomalously high or low amplitude dominating the segmentation process. The normalized amplitude, A norm (x,t), of a sample at time t in trace x, is given by: where A(x,t) is the observed amplitude at coordinate (x,t), and A max (x) and A min (x) are the maximum and minimum amplitudes observed in trace x. A norm (x,t) is then restricted to values between −1 and 1.
Phase 2: data segmentation and feature extraction Thresholding filter (Fig. 5b). Thresholding performs the most extensive segmentation of GPR data, enhancing reflections of interest and suppressing superfluous data noise. Each A norm (x,t) is compared to a threshold T, where 0 < T < 1, and allocated a binary value A bw (x,t) of either 0 or 1 according to (Zakeri et al., 2017): The definition of an appropriate threshold can be conducted via automated or user-defined means (Huang et al., 2015). User-defined thresholds give flexibility, but several must be trialled for each image, which can be time-consuming. Fig 6(a) shows the impact of different thresholds for segmenting the Figure 5. Schematic GPR data segmentation framework. Phase 1: ground-penetrating radar data conversion. Normalization of the GPR data (a). Phase 2: data segmentation and feature extraction. Application of the image segmentation filters, including (b) thresholding and (c) edge detection. Phase 3: data visualization and analysis. Analysis of the image segmentation filters performance and key reflection events (d,e). Phase 4: statistical analysis of segmented data. Groundtruth data (f) are compared with the variation of amplitudes within a segmented region; the causes of the underlying reflectivity can be suggestede.g. whether reflections correspond to stratigraphic, lithological or hydrological changes (g), and this can then be used to reconstruct reflectivity logs from GPR reflectivity data (Supporting Information Data 3). [Color figure can be viewed at wileyonlinelibrary.com]. same depth slice. However, automated thresholding is applied for the analysis in this paper: following normalization, the standard deviations of the data are calculated, removing many of the outlying data points. This process creates a measure of the variability in the data and assigns a single value for each trace in the timeslice. The mean of the standard deviation is calculated, and a unique threshold value is created for each timeslice (Fig. 6b).
Edge detection filter (Fig. 5c). Edge detection defines the edges of features in images, while simultaneously suppressing the superfluous noise around and between them. A 3 × 3 window of A norm (x,t) is multiplied by a 3 × 3 gradient matrix, 1 0 1 2 0 2 1 0 1 to give A y (x,t). If the sum of A y (x,t), A edge , is above the edge detection threshold then the central element of A norm (x,t) is deemed to be an edge. As an example, the following A norm (x,t): The sum of A y (x,t) is non-zero, =0.38.
Two of the most common edge detection filters are the Canny and Sobel filters (Öztürk and Akdemir, 2015). The Sobel filter works by calculating intensity gradients at each trace. It finds the direction of the largest increase from low to high and the rate of change in that direction (Lamba and Raj, 2017). The resulting plot shows how quickly or smoothly the timeslice data change at each trace, and how likely it is that the trace represents an edge. The Canny filter removes noise with a low-pass filter first and then performs a nonmaximum suppression to pick out the best data points for edges when there are multiple possibilities. The Sobel and Canny filter were both trialled and the Sobel filter has been judged as the most appropriate, due to its computational simplicity.

Phase 3: data visualization and analysis
Reflection pattern visualization (Fig. 5d,e). The segmented patterns are then overlaid on the original GPR data. This allows the user to assess the match between the segmented data and the GPR facies and assess whether a different threshold may be required.
Phase 4: quantitative analysis of segmented data GPR and core data integration (Fig. 5f,g). Although the segmentation algorithms initially provide qualitative data interpretations, they can be used to yield qualitative image attributes. By considering the variation of amplitudes within a segmented region, the causes of the underlying reflectivity can be suggested, for example whether reflections correspond to stratigraphic or hydrological changes (Fig. 5f). In this analysis, we compare 'reflectivity logs' extracted from the segmented GPR volumes to our core data, to allow reflectivity to be explained and therefore interpreted across the 3D data volume.

Results
Our recorded GPR data volume comprises 379 timeslices, each of which is passed through the segmentation routine. Using a velocity of 0.075 m ns −1 (Table 1), the 0.2-ns interval between each slice corresponds to a depth interval of 0.075 m; the maximum depth sampled is~2.7 m. Signal-to-noise ratios remain high throughout the depth range of the dataset. The algorithm completed the segmentation of the dataset in 1 min 50 s on a desktop PC with 16 Gb of RAM. The segmentation algorithm categorizes the timeslice amplitude data into two categories: (i) A bw (x,t) = 1 is coloured white, and (ii) A bw (x,t) = 0 is coloured grey (Fig. 7). Figure 7 shows six representative depth slices, at depth increments of 0.5 m, to show the variation in reflectivity with depth (the full output is shown in the Supporting Information). The segmented depth slices show an evolving pattern of reflectivity, indicating patterns of subsurface inhomogeneity. These features are less easily perceived in the depth slices, showing the efficacy of the segmentation approach. The depth interval 0-1 m is dominated by linear features, interpreted as shallow airfield infrastructure; the sensitivity of segmentation to such features highlights the potential of the method for (e.g.) civil engineering applications, but further investigation of these responses is beyond the scope of the current analysis. By contrast, the depth interval 1.5-2.5 m reveals a more complex pattern of reflectivity; the slice at 1.5 m is dominated by a channel-like feature, which intersects a rather circular feature. In the 2-m slice, the morphology is dominated by thin, meandering, features, similar in morphology to small channels. In the deeper intervals, the channel-like features become thinner and sparser.
The variation of features with depth can be interpreted as evidence of the evolving dune system through time. To identify the origins of GPR reflectivity at this site, the segmented profiles are compared to observations from cores acquired at four locations (Fig. 8). A similar stratigraphy is observed in each core, classified into four main units: • Unit 1dark brown-red soil humus, • Unit 2yellowy-brown fine sand, • Unit 3dark grey silty sand, and • Unit 4mid-dark grey organic silty sand.
The transition between Units 1 and 2 does not appear to provide a substantial GPR reflection; additionally, it could be too shallow to be resolved by the 500 MHz wave or it is obscured by the arrival of the direct air-wave (an airborne arrival that is unavoidable in common-offset GPR records). Unit 3 is rather homogenous; its depth extent typically correlates with a region of segmented responses in which A bw (x,t) = 0, consistent with the lack of substantial dielectric contrast within the unit. By contrast, Unit 4 shows inclusions of organic matter, which represent a dielectric contrast. Given the greater internal heterogeneity of this layer, the segmented GPR slices show a greater tendency towards A bw (x,t) = 1. By linking these observations, it appears that the preserved GPR responses correspond to greater concentrations of organic matter, hence mapping where the segmentation A bw (x,t) = 1 yields a proxy for the depth and spatial variation of organic matter.
Reflectivity-logs of GPR reflection patterns are created and compared to core logs (Fig. 9). For each depth slice, these logs describe the segmented percentage of a trial window around the core location; the red curves in Fig. 9 consider the entire slice as the effective window (defining an overall 'site trend'), whereas the blue curves assume a window of 25 × 3 m at the analysis location (defining a 'local trend' at the core location). The dimensions of this window encompass the spatial resolution of the GPR wavelet and the positional precision of the GPR grid and therefore the core location. A higher percentage within the window at any depth suggests a greater degree of scatter within a unit, here interpreted to be caused by organic deposits.
The local trend of GPR reflectivity-log 1 (blue, Fig. 9a) shows little variation in segmentation percentage with depth: values peak at 22%, at 1.65 m, but are < 10% in the depth range sampled by the core. The core stratigraphy is dominated by the dark grey sand of Unit 3, and the organic-rich Unit 4 is absent in the depth range of the downcore (unit 3). By contrast, logs 2, 3 and 4 ( Fig. 9b-d) show two discrete increases in the percentage of segmented values, that typically exceed 30% and occasionally exceed 40-50%. The shallower of these, between 1.1 and 1.4 m depth, correlates well with the location of Unit 4 and, therefore, the presence of organic matter in the corresponding cores. Only core 3 reveals the deeper appearance of Unit 3, terminating in the unit at~2 m depth.
Comparing local and site-wide reflectivity-logs gives an indication of the depositional variability across the field area. Although the site-wide trends (red, Fig. 9) support there being two periods of organic deposition at the site, these are rather more diffuse in depth than the local trends. It is therefore unlikely that the same depositional setting was active across the entire study site at any one time, which is consistent with the spatial variation seen in segmented depth slices in Fig. 7. However, the depth span of the upper Unit 4 (1.1-1.4 m) is less than its deeper appearance (1.7-2.5 m).

Llanbedr dunefield reconstruction
Dunefields are common landforms along many coastlines, and their stratigraphy often reveals important indications about the sensitivity of aeolian landforms to climatic perturbations (Arbogast et al., 2004;Miot da Silva et al., 2013;Miot da Silva and Shulmeister, 2016;Mu et al., 2016). GPR images have often aided the study of landscape change, particularly where older periods of aeolian activity are unknown (e.g. Bristow et al., 2000aBristow et al., ,b, 2010Girardi and Davis, 2010;Trivedi et al., 2012;Zhao et al., 2018). However, the study of buried aeolian dunefields has proven to be much more difficult, even more so in areas where subsequent land use has altered the surface landscape (i.e. urban developments). The application of 3D GPR, coupled with an image segmentation-based analysis algorithm, has provided a methodology with which to study buried palaeo-environments in greater detail. The results of this investigation have revealed a late Holocene dunefield that has undergone periods of stability and migration (Fig. 10). The segmented GPR and core data suggest two phases of dune stability, interspersed with episodic dune migration. Phase 1 is inferred between 1.7 and 2.5 m depth, interpreted as the earlier of the two periods of dune vegetation growth, shown by thin, sparsely distributed dendritic reflection morphologies, slowly beginning to establish across the wider dunefield. Evidence of the termination of the first phase of dune stabilization is observed at 1.6 m depth, where organic matter content (and thus GPR reflectivity) is reduced; dune migration destroyed the earlier vegetation and halted organic deposition, but it is reasonable to assume that the unit would have been encountered in cores 2 and 4 had they sampled deeper. This could imply two phases of organic-rich deposition at the site. Phase 2 is inferred between 1.0 and 1.4 m depth, where GPR responses indicate a different morphology that involves irregular circular deposits connected via a large 'channel' and suggesting that the channel sediments may be more permeable (Mountney and Russell, 2009). Taking depth as a proxy for time since deposition, this may suggest that the more recent deposition was more synchronous across the site and certainly between the four core locations. The greater spatial concentration of reflections in this phase suggests a wider coverage of dune vegetation, potentially consistent with rapid and synchronous colonization across a greater area. At shallower depths, there is no further evidence of the evolving dune system owing to the airfield development and the installation of infrastructure.
The main periods of sand deposition across North Wales took place during the Little Ice Age (LIA) between 1300 and 1800 (Lamb, 2005). Bailey (2003) summarizes seven main periods of sand deposition across the coastline of Wales over the past 1000 years: AD 1260AD -12801300-13401450-15501560-16101630-17501820-1870and1890-1980. Although no chronological data are available from Llanbedr, age control has been acquired from the 10-km-long parabolic dune system at Morfa Harlech, which is situated approximately 8 km north of the study area. Bailey et al. (2001), Bailey (2003) and Doody (2013) use optically stimulated luminscence to establish that sand dune deposition in this area started from AD 1385, which is commensurate with the development of transgressive dune fields around Europe during the LIA (Jackson et al., 2019). Furthermore, Bailey et al. (2001) dated the deepest sand unit at Aberffraw, Anglessy (45 km away) to AD 1180-1460. This, within the uncertainty range, overlaps with the period of initial sand establishment documented at Morfa Harlech. It is therefore reasonable to assume the establishment of sand dunes at Llanbedr also occurred from the 14th century onwards (Doody, 2013) with the~2.8 m of sand buried under the airfield recording one or more phases of sand deposition over the past 700 years, similar to that recorded elsewhere in North Wales. The area surveyed within Llanbedr airfield was part of a large, low-relief hummocky dunefield (Pye et al., 2007a,b). The buried hummocky dunefield formed behind the parabolic dunes of Morfa Dyffryn, still extant, which acted as a barrier, protecting the area from extreme weather events. This relative stability allowed for dune vegetation to expand, resulting in stabilization of the wider dunefield (Pye et al., 2014).

Methodological improvements
Integration of core data into geophysical data alongside image segmentation techniques is a powerful way to develop detailed Figure 9. Core logs (left-hand side) compared to percentage above threshold GPR responses versus depth (right-hand side of each sub-panel). The 'local' log (blue) represents windows of segmented data at four core locations. The 'site' log (red) assumes the whole area of the depth slice as the analysis window. Logs 2, 3 and 4 display similar relationships with high log values appearing between 1-1.4 and 1.8-2.2 m, which correlate with the appearance of organic material in the gouge core. Core 1 shows no organic content within the range of the GPR reflectivity-log. [Color figure can be viewed at wileyonlinelibrary.com].
reconstructions of buried landscapes. Traditional coring approaches are limited by 2D data coverage (Neal and Roberts, 2000); they require logs to be interpolated between sparse control points, thus limiting the accuracy of the interpretation. Although GPR provides improved spatial coverage, interpretations are limited in the absence of core data. Here, GPR and core data are integrated: GPR allows a coring strategy to target key horizons and to avoid damaging vulnerable infrastructure; the combination of the two allows temporal and spatial variations in deposition (e.g. horizon depth, thickness and continuity) to be assessed. Additionally, the two datasets are compared quantitatively via the unique application of GPR segmentation approaches. The segmentation approach is rapidly and objectively implemented, facilitating a first-pass appraisal of data and/or further comparison to core logs. Having calibrated GPR reflectivitylogs against core observations, inferences can be made in regions without core control, extending interpretations either deeper than the real core sample depth or into areas where coring was impossible (e.g. due to modern infrastructure or an Figure 10. Schematic reconstruction of the evolution of the Llanbedr dunefield. Six key timeslices have been collated from the GPR volume to reflect the evolution of the Llanbedr dunefield. The segmented interpretation (i) has been divided into two key evolutionary phases. The resulting interpretation and modelling of the environment (ii) has utilized segmented reflection patterns and core data to create a schematic evolution of the coastal dunefield. [Color figure can be viewed at wileyonlinelibrary.com]. impenetrable ground surface) (see Supporting Information Data 3 for reflectivity GPR core logs reconstructed from segmented GPR data).
The ability to characterize deep reflections has led to an improved understanding of the Llanbedr dune system, with our results supporting the observations of Pye et al. (2007aPye et al. ( ,b, 2014. The proposed methodology has characterized the spatial variation in deposition of organic-rich sands, revealing different stages of dune stability, from a now-buried coastal dunefield, most probably formed over the past~700 years (Bailey et al., 2001;Bailey, 2003;Bailey & Bristow, 2004), the extent of which was previously not known to extend beneath the northern perimeter of the airfield. The results of this investigation demonstrate the value of an integrated sedimentological, geophysical and image segmentation approach to visualizing, categorizing and groundtruthing complex buried coastal landscapes, and could also be applied to other geomorphological study areas.

Conclusions
We present a novel methodology utilizing image segmentation technology, to automatically classify GPR reflection patterns using a limited core dataset, aiding the reconstruction of a buried aeolian landscape. Our results demonstrate that complex GPR reflection patterns can be extracted in a largely automated manner, highlighting to the user instantaneously the key reflection anomaly within their data, allowing for phases in subsurface palaeoenvironments to be modelled, as well as the spatial and temporal variations in their distribution.
GPR is well-suited to characterizing the 3D architecture and depositional history of dune palaeosols, providing a noninvasive means for mapping systems that are otherwise inaccessible. At sites such as Llanbedr, which have undergone development and the installation of infrastructure, sites where coring is possible are restricted. Hence the non-invasive nature of GPR is particularly powerful and offers a key to understanding the palaeo landscape. Importantly, they provide the ability to visualize the data in 3D, allowing investigation into the evolution of the landscape (e.g. buried or urban environment), something that is much more difficult when utilizing a core-based approach.
The image segmentation method for semi-automated GPR analysis presented here has been instrumental in advancing the understanding of the Holocene coastal change on the coast of North Wales. The segmented GPR reflections coupled with quantitative analysis have provided a hypothesized evolution of this previously highly active aeolian system. The extent of a former hummocky dune region was unknown, as was its evolutionary history. The identified periods of dune stability and migration highlight an active aeolian system, which responded dynamically. The framework set out in this investigation utilizing GPR, image segmentation tools and sedimentological analysis can assist and advance the traditional workflow, characterizing buried coastal environments in greater detail.

Supporting information
Additional supporting information can be found in the online version of this article. Figure S1. Velocity analysis of the Llanbedr GPR dataset utilizing semblance analysis outlined in Booth et al. (2010). V ST refers to the best-fit velocity, with precision range. t(x 0 ) refers to the travel-time, in nanoseconds. To the apex of the hyperbola. A velocity of 0.075 m ns −1 ± 8% was used to depthconvert the Llanbedr dataset. Figure S2. (a) Simplified flow chart describing the flow of steps within the data segmentation algorithm; (b) implified flow chart describing the GPR data conversion procedure. Figure S3. (a) Core logs (left-hand side) compared to percentage above threshold GPR responses versus depth (right-hand side of each sub-panel). (b) Site map highlighting the location of the theoretical groundcores. (c) Pseudo GPR sedimentary structures reconstructed utilizing the responses from the segmented timeslice data. The example shows how the method can be used to possibly reconstruct ground conditions inaccessible to conventional coring. Data S4. Image segmentation scripts https://doi.org/10.6084/ m9.figshare.c.5071436.