Processing and performance of topobathymetric lidar data for geomorphometric and morphological classification in a high-energy tidal environment

. The transition zone between land and water is difﬁcult to map with conventional geophysical systems due to shallow water depth and often challenging environmental conditions. The emerging technology of airborne topobathymetric light detection and ranging (lidar) is capable of providing both topographic and bathymetric elevation information, using only a single green laser, resulting in a seamless coverage of the land–water transition zone. However, there is no transparent and reproducible method for processing green topobathymetric lidar data into a digital elevation model (DEM). The general processing steps involve data ﬁltering, water surface detection and refraction correction. Speciﬁcally, the procedure of water surface detection and modelling, solely using green laser lidar data, has not previously been described in detail for tidal environments. The aim of this study was to ﬁll this gap of knowledge by developing a step-by-step procedure for making a digital water surface model (DWSM) using the laser lidar data.

Abstract. The transition zone between land and water is difficult to map with conventional geophysical systems due to shallow water depth and often challenging environmental conditions. The emerging technology of airborne topobathymetric light detection and ranging (lidar) is capable of providing both topographic and bathymetric elevation information, using only a single green laser, resulting in a seamless coverage of the land-water transition zone. However, there is no transparent and reproducible method for processing green topobathymetric lidar data into a digital elevation model (DEM). The general processing steps involve data filtering, water surface detection and refraction correction. Specifically, the procedure of water surface detection and modelling, solely using green laser lidar data, has not previously been described in detail for tidal environments. The aim of this study was to fill this gap of knowledge by developing a step-by-step procedure for making a digital water surface model (DWSM) using the green laser lidar data. The detailed description of the processing procedure augments its reliability, makes it user-friendly and repeatable. A DEM was obtained from the processed topobathymetric lidar data collected in spring 2014 from the Knudedyb tidal inlet system in the Danish Wadden Sea. The vertical accuracy of the lidar data is determined to ±8 cm at a 95 % confidence level, and the horizontal accuracy is determined as the mean error to ±10 cm. The lidar technique is found capable of detecting features with a size of less than 1 m 2 . The derived high-resolution DEM was applied for detection and classification of geomorphometric and morphological features within the natural environment of the study area. Initially, the bathymetric position index (BPI) and the slope of the DEM were used to make a continuous classification of the geomorphometry. Subsequently, stage (or elevation in relation to tidal range) and a combination of statistical neighbourhood analyses (moving average and standard deviation) with varying window sizes, combined with the DEM slope, were used to classify the study area into six specific types of morphological features (i.e. subtidal channel, intertidal flat, intertidal creek, linear bar, swash bar and beach dune). The developed classification method is adapted and applied to a specific case, but it can also be implemented in other cases and environments. sailing, swimming, hiking, diving and surfing. In addition to human exploitation, climate change also poses a future threat with a predicted rising sea level and increasing storm intensity and frequency, expected to cause erosion and flooding in the coastal zone (Mousavi et al., 2011). All these pressures and different interests underpin the societal need for highresolution mapping, monitoring and sustainable management in the coastal zone.
Historically, the transition zones between land and water have been difficult or even impossible to map and investigate in high spatial resolution due to the difficulties in collecting data in these challenging, high-energy environments. The airborne near-infrared (NIR) light detection and ranging (lidar) is a technique often used for measuring highresolution topography; however, the NIR laser is incapable of measuring bathymetry due to the absorption and reflection of the laser light at the water surface. Traditionally, high-resolution bathymetry is measured with a multibeam echosounder (MBES) system mounted on a vessel, but it does not cover the bathymetry in the shallow water due to the vessel draft limitation.
NIR lidar and MBES are applied in different environments; however, the data are very similar and the processed high-resolution topography/bathymetry is often captured, visualized and analysed in a digital elevation model (DEM). The processed DEM may be applied for various purposes, e.g. for geomorphological mapping. Previous studies classifying morphology in either terrestrial or marine environments have been performed numerous times (Al-Hamdani et al., 2008;Cavalli and Marchi, 2008;Hogg et al., 2016;Höfle and Rutzinger, 2011;Ismail et al., 2015;Kaskela et al., 2012;Lecours et al., 2016;Sacchetti et al., 2011). These classification studies generally focus on either the marine or the terrestrial environment and do not cover the fine-scale morphology in the shallow water, due to the challenging environmental conditions. To overcome this impediment, a new generation of airborne green topobathymetric lidar that enables high-resolution measurements of both topography and shallow bathymetry has been introduced (Guenther, 1985;Jensen, 2009;Pe'eri and Long, 2011). The potential of merging morphological classifications of marine and terrestrial environments enables a holistic approach for managing the coastal zone.
The raw topobathymetric lidar measurements are spatially visualized as points in a point cloud, with each point containing information of its location and elevation. The point cloud must be piped through a series of steps before it can be visualized as a DEM. Most of the processing steps required to process raw topobathymetric lidar data into a DEM are similar to the processing steps of topographic lidar data (Huising and Gomes Pereira, 1998). However, additional processing steps are required for topobathymetric lidar data due to the refraction of the laser beam at the water surface. All submerged lidar points have to be corrected for refraction; therefore, the water depth must be known for each point.
This sets a requirement for making a digital water surface model (DWSM), before the refraction correction can be performed.
Often, the water surface is detected and modelled from simultaneous collection of green and NIR lidar measurements, where the green laser reflects from the seabed and the NIR laser reflects from the air-water interface, and the NIR laser data are then used to detect and model the water surface (Allouis et al., 2010;Collin et al., 2008;Guenther, 2007;Parker and Sinclair, 2012). The use of NIR lidar data for water surface detection has been applied in several studies. For instance, Höfle et al. (2009) proposed a method for mapping water surfaces based on the geometrical and intensity information from NIR lidar data. Su and Gibeaut (2009) classified water points from NIR lidar based on point density, intensity and elevation. They identified the shoreline based on the large sudden decrease in NIR lidar intensity values when going from land to water. Brzank et al. (2008) used the same three variables (point density, intensity and elevation) in a supervised fuzzy classification to detect the water surface in a section of the Wadden Sea. Another study in the Wadden Sea by Schmidt et al. (2012) used a range of geometric characteristics as well as intensity values to classify water points from NIR lidar data.
The capability of NIR lidar data for detecting the water surface is thus well documented. However, deriving all the information from a single green lidar dataset would be a more effective solution for water surface and seabed detection, with respect to the financial expenses and the difficulties of storing and handling often very large amounts of data. However, there is no definitive method for making a DWSM from green topobathymetric lidar data. For this purpose, the Austrian lidar company RIEGL have developed a software, RiHYDRO (RIEGL, 2015), in which it is possible to model the water surface in a two-step approach: (1) classification of water surface points based on areas with two layers (water surface and seabed) and extending the classification to the entire water body, and (2) generation of a geometric gridded DWSM for each flight swath based on the classified water surface points. However, RiHYDRO is commercial software, and thus the algorithms, which form the basis of the classification and water surface modelling, are not publicly available. Other software packages, such as Hy-droFusion (Optech, 2013) and LiDAR Survey Studio (Leica, 2015), also proclaim to have incorporated methods for the entire data processing workflow, but the algorithms in these software packages are also closed and cannot be accessed by public users.
Only few research studies have investigated the potential of water surface detection and modelling from green lidar data. In a relatively recent publication, Guenther et al. (2000) even regarded water surface detection from green lidar data as unacceptable and they justified it with two fundamental issues: (1) no water surface returns are detected in the dead zone, and (2) there is uncertainty of the water surface eleva-tion, because the green water surface returns are actually a mix of returns from the air-water interface and from volume backscatter returns, and they are generally found as a cloud of points below the water surface. Mandlburger et al. (2013) addressed the second issue by comparing the water surface points of NIR and green lidar data, and they concluded that it is possible to derive the water surface elevation from the green lidar data with sub-decimetre vertical precision relative to a reference water surface derived by the NIR lidar data. However, their work addressed only the determination of the water surface elevation without going into detail on the actual procedure of generating a DWSM. An approach for modelling the water surface from green lidar data was presented by Mandlburger et al. (2015), who did their study in a riverine environment with only few return signals from the water surface. Their method was based on manual estimates of the water level in a series of river cross sections, after which interpolation between the cross sections filled out the gaps with no water surface points to derive a continuous water surface model.
The aim of this study was to investigate the potential of improving the procedure of processing green lidar data and generating DEMs in tidal environments, and of improving the classification of morphological units in such environments. More specifically, the objectives were 1. to develop a robust, repeatable and user-friendly processing procedure of raw green lidar data for generating high-resolution DEMs in land-water transition zones; 2. to quantify the accuracy and precision of the green lidar data based on object detection; and 3. to automatically classify morphological units based on geomorphometric analyses of the generated DEM.
The investigations were based on studies undertaken in a section of the Knudedyb tidal inlet system in the Danish Wadden Sea.

Study area
The Knudedyb tidal inlet system is located between the barrier islands of Fanø and Mandø in the Danish Wadden Sea (Fig. 1a). The tidal inlet system is a natural environment without larger influence from human activity. The tides in the area are semi-diurnal, with a mean tidal range of 1.6 m, and the tidal prism is on the order of 175 × 10 6 m 3 (Pedersen and Bartholdy, 2006). The main channel is approximately 1 km wide and with an average water depth of approximately 15 m (Lefebvre et al., 2013).
The study site is an elongated 3.2 km 2 (0.85 km × 4 km) section of the Knudedyb tidal inlet system (Fig. 1b). The section is located perpendicular to the main channel and stretches across both topography and bathymetry. The study site extends northward into an area on Fanø with dispersed cottages (Fig. 1c). The most prominent morphological features within the study site include beach dunes (Fig. 1d), small mounds (Fig. 1e), swash bars ( Fig. 1f and g) and linear bars (Fig. 1h).

Topobathymetric lidar
The topobathymetric lidar technique is based on continuous measurements of the distance between an airplane and the ground/seabed. The distance (or range) is calculated by half the travel time of a laser beam going from the airplane to the surface of the earth and back to the airplane. The wavelength of the laser beam is in the green spectrum, usually 532 nm, since this wavelength is found to attenuate the least in the water column, resulting in the largest penetration depth of the laser (Jensen, 2009). In literature, topobathymetric lidar data are sometimes referred to as either bathymetric lidar or airborne lidar bathymetry. These are different terms with the same meaning, and in this paper, topobathymetric lidar is preferred, since it describes the system's ability to simultaneously measure bathymetry as well as topography.
A single laser beam may encounter many targets of varying nature on its way from the airplane and back again, and different processes are influencing the laser beam propagation through air and water. First, the laser beam may be reflected by targets in the air, such as birds or dust particles, and these can show up as lidar points in the space between the airplane and the surface. When encountering water, the speed of the laser decreases from 3 × 10 8 to, e.g. 2.25 × 10 8 m s −1 in 10 • C freshwater or, e.g. 2.24 × 10 8 m s −1 in 10 • C saltwater of 30 PSU (Millard and Seaver, 1990).
The changing speed of the laser beam also affects the direction of the laser beam when penetrating the water surface with an angle different from nadir ( Fig. 2) (Guenther, 2007;Jensen, 2009). The laser beam will be refracted according to Snell's law: sin α air sin α water = c air c water = n water n air , where α air is the incidence angle of the laser beam relative to the normal vector of the water surface and α water is the refraction angle in water. n water and n air are the refractive indices of water and air, respectively . The penetration depth in water is limited by the attenuation of the laser beam. Water molecules, suspended sediment and dissolved material all act on the laser beam by absorption and scattering, resulting in substantial reduction in power as the signal propagates into the water (Guenther, 2007;Mandlburger et al., 2013;Steinbacher et al., 2012). The laser beam also diverges in the water column, resulting in a wider laser beam footprint (Guenther et al., 2000), and this effect reduces the resolving capability of fine-scale morphology the deeper the laser beam penetrates. The returned signal is represented as a distribution of energy over time, also called the "full-waveform" distribution (Alexander, 2010;Chauve et al., 2007;Mallet and Bretar, 2009). The peaks in the full-waveform distribution are detected as individual targets encountered by the propagating laser beam. If the laser hits two targets with a small verti-cal difference, such as a water surface and seabed in very shallow water, then the two peaks in the full-waveform distribution may merge together, resulting in the detection of only one target (Fig. 2). This results in a detection minimum of successive returns from a single laser pulse, referred to as the "dead zone" (Mandlburger et al., 2011;Nayegandhi et al., 2009). The dead zone is a clear limitation to the lidar measurements, which is an important parameter to consider in very shallow water, such as intertidal environments.

Surveys and instruments
Lidar data and orthophotos were collected by Airborne Hydro Mapping GmbH (AHM) during two surveys on 19 April 2014 and 30 May 2014.
On 19 April 2014, the quality of the lidar data was validated at two sites along Ribe Vesterå River ( Fig. 1i and j): -Validation site 1, with a 2.50 × 1.25 × 0.80 m cement block, is located on land next to the mouth of Ribe Vesterå River (Fig. 1i). The block was covered by seven swaths retaining 227 lidar points from the block surface, which were used for assessing the accuracy and precision of the lidar data.
-Validation site 2, with a 0.92 × 0.92 × 0.30 m steel frame, is located in the Ribe Vesterå River, its top situated just below the water surface (Fig. 1j). The frame was covered by four swaths retaining 46 lidar points from the surface of the frame, which were used for precision assessment, and for testing the feature detection capability of the lidar system. According to the International Hydrographic Organization survey standards, cubic features of at least 1 m 2 should be detectable in special order areas, which are areas with very shallow water as in the study site (IHO, 2008).
Ground control points (GCPs) were measured for the four corners of the block with an accuracy better than 2 cm using a Trimble R8 RTK GPS. Measurements were repeated three times and averaged to minimize errors caused by measurement uncertainties. GCPs were also collected for the frame; however, during the lidar survey the frame experienced an unforeseen intervention by local fishermen using the frame as fishing platform. Therefore, the frame was only used to assess the deviation between the lidar points (the precision) and not to assess the deviation between the lidar points and the GCPs (the accuracy). On 30 May 2014, the study site was covered by 11 swaths, which were used for generating the DWSM and DEM. The overflight was carried out during low tide, and the water level was measured to −1 m DVR90 at Grådyb Barre, approximately 20 km NW of the study site.
The weather conditions were similar during the two surveys, with sunny conditions, average wind velocities of 7-8 m s −1 (DMI, 2014) and significant wave heights, measured west of Fanø at 15 m water, of approximately 0.5 m coming from NW (DCA, 2014). However, the waves in the lessexposed Knudedyb tidal inlet were observed in the 30 May lidar point cloud to be 0.2-0.3 m, which can be explained by the location of the study site in lee of the western-most intertidal flats and the ebb-tidal delta. The wave heights in the rest Figure 2. Conceptual sketch of the laser beam propagation and return signals. The beam refracts upon entering the water body, and it diverges as it propagates through the water column. Return signals are produced in the air, at the water surface, in the water column and at the seabed. The lidar instrument has limited capability in very shallow water (the dead zone in the figure) because the successive peaks from the water surface and the seabed are not individually separated in time and amplitude. Only the largest peak, which is from the seabed, is detected. of the study site (flood channel and intertidal ponds) were in the scale of sub-decimetres. There were no waves at validation site 2 during the 19 April lidar survey.
Lidar data were collected with a RIEGL VQ-820-G topobathymetric airborne laser scanner in both surveys (RIEGL, 2014). The scanner was characterized by emitting green laser pulses with 532 nm wavelength and 1 ns pulse width. It had a very high laser pulse repetition rate of up to 520 000 Hz. The flying altitude was 400 m, which, combined with a beam divergence of 1 mrad, created a laser beam footprint of 40 cm diameter at the ground. The high repetition rate and narrow footprint made it well suited to capture fine-scale landforms (Doneus et al., 2013;Mandlburger et al., 2011). An arc-shaped scan pattern results in a swath width of approximately 400 m, while maintaining an almost constant 20 • (±1 • ) incidence angle of the laser beam when penetrating the water surface (Niemeyer and Soergel, 2013). The typical water depth penetration of the laser scanner is claimed by the manufacturer to be one Secchi disc depth (RIEGL, 2014).
For each returned signal, the collected lidar data contained information of x, y and z, as well as a GPS time stamp and values of the amplitude, reflectance, return number, attribute and laser beam deviation (RIEGL, 2012). The essential processing steps, which are standard procedure when processing topobathymetric lidar data, were followed to produce a DEM in the study area. These steps included 1. determination of flight trajectory; 2. boresight calibration: calculating internal scanner calibration; 3. collecting topobathymetric lidar data; 4. swath alignment based on boresight calibration: the bias between individual swaths was minimized; 5. filtering: the raw data contained noise located both above and below ground, which needed to be filtered from the point cloud; 6. water surface detection: a DWSM had to be established in order to correct for refraction in the following step; 7. refraction correction: all the points below the water surface in the DWSM were corrected for the refraction of the laser beam; and 8. point cloud to DEM: the points were transformed into a gridded elevation model representing the real world terrain in the study area, including cottages and vegetation on Fanø in the northern part.
Step 1 and 2 were performed prior to the lidar survey. The different instruments (lidar, inertial measurement unit (IMU) and GPS) were integrated spatially by measuring their position relative to each other, when mounted on the airplane, and temporally by calibrating their time stamps.
Step 3 was the actual lidar survey and step 4 was the initial processing step after the lidar survey. The bias between the swaths was minimized in the software RiPRO-CESS (RIEGL LMS) by automatically searching for planes in each swath and then matching the planes between the swaths.
Steps 5-8 represent the processing of the point cloud into a DEM. The methods involved in these steps are the main focus in this study and they are described in detail in the following subsections. Each swath was pulled individually through the processing workflow to account for the continually changing water level in the study area due to tides. The broad term DEM is used, rather than the more specific terms "digital terrain model" (DTM) or "digital surface model" (DSM), because the generated model includes natural terrain in the tidal environment, which is the main focus area in this study, as well as vegetation and cottages on Fanø.

Filtering
The raw lidar data contained noise in the air column originating from the laser being scattered by birds, clouds, dust and other particles, and noise was also appearing below the ground/seabed ( Fig. 3a and b). This noise had to be filtered before further processing. The filtering process involved both automatic and manual filtering.
The automatic filtering was carried out in Hydro-Vish (AHM) with the tool "Remove flaw echoes" (Fig. 3c). The filtering tool was controlled by three variable parameters: search radius, distance and density. The search radius parameter specified the radius of a sphere in which the distance and density filters were utilized. The distance parameter rejected a point, if it was too far from any other point within the sphere. The density parameter specified the lower limit of points within the sphere. The automatic filter iterated through all the points in the point cloud. The settings for the automatic filtering were based on a sensitivity analysis of three fragments of the lidar data, and the settings were selected so that a minimum of valid points were removed by the automatic filter. The settings were the following: search radius = 1 m, distance = 0.75 m and density = 4 points.
The automatic filter could not remove two layers of noise points closely above and below ground, but on the other hand, more widely dispersed points in the deeper bathymetry were removed. To account for this, the point cloud went through manual filtering in Fledermaus (QPS) software, where the remaining noise points were rejected and the valid bathymetric points were accepted (Fig. 3d).
The filtered point cloud (with water points) was used in the following step to detect the water surface. Meanwhile, a copy of the data was undergoing additional manual filtering, removing all the water points (Fig. 3e). After this final filtering step, there were only points representing topography, bathymetry, vegetation and man-made structures left in the dataset.

Water surface detection
The water surface detection was based on determining the water surface elevation and the water surface extent, thereby producing a DWSM. The water surface elevation was determined based on the water surface points, and the extent was determined by extrapolating the water surface until it intersected the surface of the topography. Two assumptions were made in the production of the DWSM: 1. The water surface was horizontal. This was a simplification of the real world. Tidal processes and wind and wave setup may cause a sloping water surface, and the water is often topped by more or less significant waves. A linear fit through the water surface lidar points along the main channel showed a changing water level of 0.13 m over a distance of 400 m, corresponding to a 0.325 × 10 −3 (0.019 • ) sloping water surface. A simi- lar fit through the lidar points along the flood channel showed a slope of 0.156 × 10 −3 (0.009 • ). The maximum wave heights observed in the main channel were 20-30 cm. Based on the moderate slope of the water surface and relatively low wave height, the water surface was assumed to be flat. This assumption is deemed error prone, but at the time of this study, it was our best estimate.
2. The study area contained water bodies with two different water levels: one represented the water level in the main channel and the other represented the water level in the flood channel. This was also a simplification, as the tidal flat contained small ponds with potentially different water levels. However, almost all of these ponds contained no lidar points of the water surface, which means that the water depth in the ponds must have been within the limitation of the dead zone. Therefore, it was impossible to detect individual water surfaces in the ponds.
A series of processing steps was performed to produce the DWSM. The first step was to extract a shallow surface and a deep surface from the filtered point cloud (with water points) in Fledermaus (Fig. 3f). Both surfaces consisted of 0.5 m × 0.5 m cells, and the elevation of each cell was equal to the highest point within the cell (shallow surface) and the lowest point within the cell (deep surface), respectively. The shallow surface should then display the topography along with the water surface, whereas the deep surface should dis-play the topography and the seabed (as long as the seabed was detected by the laser). It is worth noting that the extraction of the shallow surface and the deep surface had nothing to do with the final DEM, as they were merely intermediate steps performed for the water surface detection.
The following steps were focused on the shallow surface to determine the elevation of the water surface (Fig. 3g). First, the shallow surface was down sampled to a surface with a cell size of 2 m × 2 m, and the new cells were populated with the maximum elevation of the input cells. The down sampling was done for smoothing the water surface and thereby eliminating most of the outliers. The exact cell size of 2 m × 2 m, as well as populating them with the maximum value, was chosen based on the work by Mandlburger et al. (2013). They compared water surface detection capability between green lidar data, collected with the same RIEGL-VQ-820-G laser scanner, and NIR lidar data, which were assumed to capture the true water surface. They found that the green lidar generally underestimated the water surface elevation, but that reliable results were achieved by increasing the cell size and only taking the top 1-5 % of water points into account. According to their work, it was assumed that placing the water surface on the highest points in 2 m cells provided a good estimate of the true water level. However, based on their results it could be expected that the water surface elevation in this case would be underestimated on the order of 2-4 cm.
The water-covered areas in the main channel and the flood channel were manually extracted from the newly downsampled shallow surface. The average elevation of the 2 m cells within each area was calculated, and these values constituted the elevation of the water surfaces in the main channel and in the flood channel, respectively. Two horizontal water surfaces were created in the flood channel and the main channel with a cell size of 0.5 m × 0.5 m and cell values equal to the calculated water surface elevations. The high spatial resolution of 0.5 m cells was chosen to produce a detailed DWSM along the edges of the land-water transition.
Finally, the extent of the water surfaces was determined by subtracting the deep surface cell elevations from the water surface elevation and discarding all cells with resulting negative values (Fig. 3h), together forming the DWSM.

Refraction correction
The refraction correction of all the points below the DWSM was calculated in HydroVish (AHM). The input parameters were the filtered point cloud (without water points), the derived DWSM and the trajectory data of the airplane. The refraction correction was calculated automatically for each point based on the water depth, the incident angle of the laser beam and the refracted angle according to Snell's law (Eq. 1 and Fig. 3i).

Point cloud to DEM
After iterating through the processes of filtering, water surface detection and refraction correction for all the individual swaths, the lidar points of all swaths were combined. The transformation from point cloud into a DEM was performed with ArcGIS (ESRI) software. The DEM was created as a raster surface with a cell size of 0.5 m × 0.5 m, and each cell was attributed the average elevation of the points within the cell boundaries. It was chosen to make the resolution of the DEM lower than the laser beam footprint size (i.e. 40 cm), due to the inaccuracies arising from attributing smaller cells with measured elevation values spanning across a larger area. Furthermore, the 0.5 m cell size was chosen to get as high resolution as possible without making any significant interpolation between the measurements. In this way, each cell represented actually measured elevations instead of interpolated values. However, there were still very few gaps of individual cells with no data in the resulting raster surface in areas with relatively low point density. Despite the general intention of avoiding interpolation, it was chosen to populate these cells with interpolated values to obtain a full coverage DEM (except for the bathymetric parts beyond the maximum laser penetration depth). The arguments for interpolation were that (1) the interpolated cells were scattered and represented only 1.7 % of all the cells, (2) they were found primarily on the tidal flat where the slope is generally less than 1 • , meaning that the elevation difference from one cell to a neighbouring cell is usually less than 1 cm, and (3) the general point density in most of the study area was so high that the loss of information by lowering the DEM resolution would represent a larger sacrifice than interpolating a few scattered cells. The interpolation was performed by assigning the average value of all neighbouring cells to the empty cells. The final DEM was thereby fully covering the topography, and the bathymetry was covered down to a depth equal to the maximum laser penetration depth.

Accuracy and precision of the topobathymetric lidar data
The term accuracy refers to the difference between a point coordinate (in this case, a lidar point) compared to its "true" coordinate measured with higher accuracy, e.g. by a total station or a differential GPS, while the term precision refers to the difference between successive point coordinates compared to their mean value, i.e. the repeatability of the measurements (Graham, 2012;Jensen, 2009;RIEGL, 2014). Two "best-fit planes" based on the lidar points on the block and the frame surfaces were established with the "Curve fitting tool" in MATLAB (MathWorks). We propose the use of these two planes to give an indication of the relative precision of the lidar measurements.
Another best-fit plane was established based on the block GPS measurements, and this plane was regarded as the true block surface for assessment of the accuracy of the lidar measurements. The established planes were described by the polynomial equation where x, y and z are coordinates and a, b and c are constants. Inserting x and y coordinates for the lidar surface points in Eq.
(2) led to a result of the corresponding elevation (z) as projected on the fitted plane. The difference between the elevation of the lidar point and the corresponding elevation on the fitted plane was used as a measure of the vertical accuracy (for the GCP-fitted plane of the block) and the vertical precision (for the lidar-point-fitted plane of the block and the frame). Statistical measures of the standard deviation (σ ), mean absolute error (E MA ) and root mean square error (E RMS ) were calculated by where z i is the elevation of the measured lidar points, z plane is the corresponding elevation on the best-fit plane and n is the number of lidar points. The vertical accuracy and precision were determined at a 95 % confidence level based on the accuracy standard presented in "Geospatial Position Accuracy Standards Part 3: National Standard for Spatial Data Accuracy" (NSSDA) (FGDC, 1998): The horizontal accuracy was determined as the horizontal mean absolute error (E MA,xy ) based on the horizontal distances between the block corners, measured with RTK GPS, and the best approximation of the block corners derived from the lidar points of the block surface. The minimum distance between a block corner and the perimeter of the lidar points was regarded as the best approximation. Hereafter, E MA,xy was calculated as the average of the four corners.

Geomorphometric and morphological classifications
The processed DEM was applied in two classification analyses: first, a geomorphometric classification and then a morphological classification. Both were based on the DEM and derivatives of the DEM, but they were differentiated by the resulting classification classes, which showed (1) surface geometry and (2) surface morphology. The analysis mode, as defined by Pike et al. (2009), was "general" in the geomorphometric classification where the surface geometry was continuously classified within the study site, while being "specific" in the morphological classification where discrete morphological units were classified. The northern part of the study site with cottages on Fanø was excluded in the classification analyses, as the objective of this work was to classify the natural terrain (geomorphometry and morphology) in the high-energy and dynamic tidal environment.

Geomorphometric classification analysis
The benthic terrain modeler (BTM) tool (Wright et al., 2005) was used for the geomorphometric classification. The tool is an extension to ArcGIS Spatial Analyst, originally used for analysing MBES data (Diesing et al., 2009;Lundblad et al., 2006;Rinehart et al., 2004). The BTM classification tool uses fine-and broad-scale bathymetric position indexes (BPIs) (Verfaillie et al., 2007) in a multiple-scale terrain analysis to classify fine-and broad-scale geometrical features. The BPIs are measures of the elevation of a cell compared to the elevation of the surrounding cells within the determined scale (radius) size. Positive BPI values indicate a higher elevation than the neighbouring cells and negative BPI values indicate a lower elevation than the neighbouring cells. For instance, a BPI value of 100 corresponds to 1 SD (standard deviation) and a value of −100 corresponds to −1 SD of the cell elevation compared to the elevation of the surrounding cells within the determined scale size. BPI values close to zero are derived from flat areas or from constant slopes.
The elevation values of each cell in the DEM were exaggerated by a factor of 10 before the classification to enable the BTM to detect the shapes of the terrain. The fine and broad scales were determined based on the BPI results for different radius sizes. The best results were obtained from a broad-scale BPI of 100 m radius and a fine-scale BPI of 10 m radius, based on visual inspection. The fine-and broadscale BPIs were used together with the slope of the actual DEM (not the exaggerated) to classify the investigated area into the geomorphometric classes: fine-scale crests, broadscale crests, depressions, slopes and flats (Fig. 4). The classification classes were decided based on previous studies using the BTM classification tool with success (Diesing et al., 2009;Lundblad et al., 2006). The thresholds for the fine-and broad-scale BPIs were in previous studies often defined as 1 SD (Lundblad et al., 2006;Verfaillie et al., 2007); however, thresholds of 0.5 SD have also previously been applied (Kaskela et al., 2012). We used a low threshold of 0.5 SD due to the generally very gentle variations in the terrain geometry of the tidal inlet system. We defined the threshold between slopes and flats as 2 • . This definition was a compromise between detecting as many slopes as possible but avoiding too many "false slopes" being detected along the swath edges, which seemed to be a consequence of lower precision at the outer beams of the swath, as well as differences between overlapping swaths.

Morphological classification analysis
A morphological classification was developed for delineating actual morphological features in the study area. The classification was built partly on different neighbourhood analyses and slopes derived from the DEM and partly on the local tidal range. Broad-scale crests from the geomorphometric classification were also incorporated in the analysis. Figure 5 describes the steps performed in ArcGIS, which led to the classification of six morphological classes: swash bars, linear bars, beach dunes, intertidal flats, intertidal creeks and subtidal channels. All the criteria for defining a particular morphological class had to be fulfilled for a cell to be classified into that class. Cells that did not meet the criteria to be classified into any of the morphological classes were assigned the class "unclassified".
A total of 33 years of continuous measurements of the water level at Havneby on Rømø, 25 km south of the study area, showed a mean low water level of −0.94 m (DVR90) and a mean high water of 0.94 m (DVR90) (Klagenberg et al., 2008). Although the tidal range in Knudedyb was probably slightly different, it was the best estimate for the study site. Therefore, these water levels were used to separate between the supratidal, intertidal and subtidal zones.
Subtidal channels were defined as everything below the mean low water, which was −0.94 m. A "smooth DEM" was created, in which each cell of the original DEM was assigned the average elevation value of its surrounding cells in a window size of 100 m × 100 m (actually 199 × 199 cells, i.e. 99.5 m × 99.5 m). The result was subtracted from the original DEM, creating an elevation change model (ECM),   which made it possible to extract information about the deviation of the cells in the DEM compared to its surrounding cells. The principle is similar to the BPI, and again the purpose was to locate cells with a higher/lower elevation than its surrounding cells. Positive values were higher cells and negative values were lower cells. Certain thresholds were found suitable for classifying beach dunes (> 0.8 m) and intertidal creeks (< −0.3 m). These two classes were further classified into their respective tidal zones (supratidal and intertidal) based on the elevation. Intertidal flats were classified by low slope values (< 1 • ) of a down-sampled 2 m DEM (each down-sampled cell was assigned the mean value of its 4 × 4 original cells). Moreover, to be classified as a flat, the ECM had to be within ±10 cm to avoid any incorrect intertidal flat classification of flat crests on top of bars or flat bottoms inside creeks or channels. The BTM classification class "broad-scale crests" was used as an input, since it was found to capture bar features. However, the thresholds used in the BTM classification resulted in capturing features larger than bars in the broad-scale crests class. To distinguish between bars and larger features, the standard deviation of each DEM cell in a moving window size of 250 m × 250 m (actually 249 × 249 cells, i.e. 124.5 m × 124.5 m) was calculated. A suitable threshold to distinguish between bars and larger features was 0.6 SD. Finally, swash bars and linear bars were identified by an area/perimeter ratio, based on the assumption that linear bars have a smaller ratio than swash bars due to the different shapes. Based on visual interpretation, a ratio of 4 was found to be a suitable threshold.

Refraction correction and dead zone extent
The vertical adjustment of the lidar points (z diff ) due to refraction correction is linearly correlated with the water depth (d) (Fig. 6). An empirical formula is derived for this relationship: A lidar point at 1 m water depth is vertically adjusted by approximately 0.23 m (Fig. 6). The variations around the linear trend in Fig. 6 are due to changing incidence angles of the laser beam that varies with the airplane attitude (roll, pitch and yaw). The vertical extent of the dead zone is approximately 28 cm, determined by plotting the vertical difference between the shallowest and the deepest lidar point within 0.5 m cells -i.e. between the shallow surface and the deep surface (Fig. 7). The difference is manifested by an abrupt change at the dead zone, and the highest rate of change is shown to be at a water depth of approximately 28 cm.

Sub-decimetre accuracy and precision
The vertical root mean square error of the lidar data is ±4.1 cm, and the accuracy is ±8.1 cm with a 95 % confidence level (Table 1 and Fig. 8a). The vertical precision of the lidar data with a 95 % confidence level is ±3.8 cm for the points on the frame and ±7.6 cm for the points on the block ( Table 1).
The horizontal accuracy calculated as the horizontal mean absolute error (E MA,xy ) is determined to ±10.4 cm (Fig. 8b).

Point density and resolution
The average point density is 20 points per m 2 , which equals an average point spacing of 20 cm ( Table 2). The point density of the individual swaths varies between 7 and 13 points per m 2 , and the point density of the combined swaths in the study area varies between 0 and 216 points per m 2 , although densities above 50 points per m 2 are rare.

DEM and landforms
The elevations in the studied section of the Knudedyb tidal inlet system range from −4 m DVR90 in the deepest parts of the flood channel and main channel to 21 m DVR90 on top of the beach dunes on Fanø (Fig. 9). Beach dunes and cottages of the village Sønderho are clearly visible in the northern part of the study site ( Fig. 9a and b). The intertidal zones are generally flat, while the most varying morphology is found in the area of the flood channel ( Fig. 9c and d) and in the area close to the main channel ( Fig. 9e and f). The flood channel is approximately 200 m wide in the western part and it divides into two channels towards the east. The bathymetry of the channel bed is clearly captured by the lidar data in the eastern part and also in the western part down to −4 m DVR90, which approximately equals a water depth of 3 m at the time of survey. An intertidal creek joins the flood channel from the north (Fig. 9d). From the flood channel towards south, the tidal flat is vaguely upward sloping, until it reaches two distinct swash bars, which are rising 0.9 m above the surrounding tidal flat, reaching a maximum elevation of 1.5 m DVR90 ( Fig. 9e and f). Further south, the linear bars along the margin of the main channel are clearly captured in the DEM (Fig. 9e).

Geomorphometric and morphological classifications
The geomorphometric and morphological classifications show that most of the study area is located in the intertidal zone and is mostly flat. This is manifested by the dominating two classes; flats and intertidal flats ( Fig. 10a and b). The geomorphometric classification identifies slopes as stripes with NNW-SSE directionality across the flats. These are following the direction of the survey lines, and thus, they are not real morphological features but more an indication of lower precision of the lidar data, especially at the outer beams of the swath. These swath artefacts are smoothed out in the morphological classification by down sampling the DEM to 2 m resolution, and therefore, the intertidal flats appear uniform and seamless. The bar features close to the main channel are well defined in the geomorphometric classification where they are classified as broad-scale crests and fine-scale crests surrounded by slopes. In the morphological classification, these are identified based on neighbourhood analyses and separated by the area/perimeter ratio into two classes, swash bars and linear bars (Fig. 10c). Broad-scale crests are also found on Fanø in the northern part of the area, and most of these are classified as beach dunes in the morphological classification. The geomorphometric classification identifies more broad-scale crests along the banks of the flood channel; however, these are not real bar features but they are identified as crests due to the nearby flood channel and creeks resulting in a positive broad-scale BPI. In the morphological classification, it is possible to distinguish between these and the actual bar features by looking at elevation deviations at an even broader scale than the broad-scale BPI. The intertidal creek in the northwestern part of the area is a mix of depressions, slopes and fine-scale crests in the geomorphometric classification, whereas it is relatively well defined and properly delineated in the morphological classification (Fig. 10d). The geomorphometric classification identifies slopes along the banks of the main channel, flood channel and the intertidal creek, as well as in front of the beach dunes and along the edges of the swash bars and linear bars. The slopes seem particularly reliable at delineating the features in the intertidal zone, i.e. swash bars, linear bars and creeks. Depressions are primarily identified in the deepest detected parts of the main channel and in the flood channel, in the intertidal creek and in the beach dunes. Fine-scale crests are found in the geomorphometric classification in locations which are high compared to its near surroundings. They are primarily seen as parts of the linear bars close to the main channel, in the beach dunes on Fanø and along the banks of the intertidal creeks.
A few small circular mounds of approximately 5 m diameter with patches of Spartina townsendii (common cord grass) located on the intertidal flat are classified as fine-scale crests in the geomorphometric classification (Fig. 11). It clearly shows the capability of capturing fine-scale features in the DEM and in the derived classification.

Performance of the water surface detection method
The water surface in topobathymetric lidar surveys are most often detected from NIR lidar data and are simultaneously collected along with the green lidar data (Collin et al., 2012;Guenther et al., 2000;Parker and Sinclair, 2012;Wang and Philpot, 2007). However, detecting the water surface and generating a DWSM based on the green lidar data alone provides a potential to perform topobathymetric surveys with just one sensor, thus optimizing the survey costs as well as data handling and storage.
The two critical issues risen by Guenther et al. (2000), as mentioned in the introduction, concerning the water surface detection with green lidar, were thoroughly investigated in this study. The first issue, regarding the gap of detected water surface signals in the dead zone, is addressed by detecting the water surface, based on areas which are known to be covered by water, and thereafter extending the water surface until it intersects the topography so that also the dead zone is covered by the modelled water surface. The second issue, regarding uncertainty in the water surface elevation determination, is addressed using the results presented by Mandlburger et al. (2013) who found a statistical relationship between the cloud of water surface points in the green lidar data and the water surface elevation derived from NIR lidar data. Mandlburger et al. (2013), however, did not describe the actual method of modelling a DWSM, which is done in this study. Mandlburger et al. (2015), on the other hand, did propose a method for modelling the water surface; however, it was done in a fluvial environment and the water level was based on manual determinations of cross-sectional water levels. The water surface detection method in this study is thus new in combining the following properties: (1) it is only using green lidar data, (2) it is based on automatic water level determination, (3) it is applied in a tidal environment (can be applied in any coastal environment) and (4) it is transparent and repeatable due to the detailed description of data processing steps given in the text.
The developed water surface detection method is new but it must be pointed out that the assumption of a flat DWSM leaves room for improvements for the future, especially if it is applied in a fluvial environment. Assuming a flat water surface is indeed a simplification of the real world, since the water surface in reality can be inclined, and it can be topped by waves.

Implications of the dead zone
The vertical extent of the dead zone is in this study determined to approximately 28 cm (Fig. 7), which means that no return signal is detected from the water surface when the water depth is less than 28 cm. The implication of the dead zone along the channel edges is minimized by extending the DWSM until it intersects the topography, but the setting is different for the small ponds on the intertidal flats. They may have different water levels than in the large channels, but no detected water surface points, since the water depth in the ponds is generally less than the vertical extent of the dead zone, i.e. approximately 28 cm. The presented method is not capable of detecting a water surface in these ponds, which means that the bottom points of the ponds are not corrected for refraction. Omitting refraction correction of a 28 cm deep pond will result in −6 cm elevation error according to the calculated refraction (Fig. 6).

Evaluation of the topobathymetric lidar data quality
The vertical accuracy of conventional topographic lidar data has previously been determined to ±10-15 cm (Hladik and Alber, 2012;Jensen, 2009;Klemas, 2013;Mallet and Bretar, 2009). Only few previous studies have focused on the accuracy of shallow water topobathymetric lidar data (Mandlburger et al., 2015;Nayegandhi et al., 2009;Steinbacher et al., 2012). Nayegandhi et al. (2009) determined the vertical E RMS of lidar data in 0-2.5 m water depth to ±10-14 cm, which is above the ±4.1 cm E RMS found in this study (Table 1). Steinbacher et al. (2012) compared topobathymetric lidar data from a RIEGL VQ-820-G laser scanner with 70 ground-surveyed river cross sections serving as reference, and found that the system's error range was ±5-10 cm, which is comparable to the ±8.1 cm accuracy found in this study. Mandlburger et al. (2015) compared ground-surveyed points from a river bed with the median of the four nearest 3-D neighbours in the lidar point cloud, and they found a standard deviation of 4.0 cm, which is almost equal to the ±4.1 cm standard deviation found in this study (Table 1). In comparison with these previous findings of lidar accuracy, the assessment of the vertical accuracy in this study indicates a good quality of the lidar data.
Mapping the full coverage of tidal environments, such as the Wadden Sea, requires a combination of topobathymetric lidar to capture topography and shallow bathymetry and MBES to capture the deeper bathymetry. The two technologies make it possible to produce seamless coverage of entire tidal basins; however, merging the two products raises the question whether the quality of the data from the two different sources is comparable. Comparing the lidar accuracy with previous findings of accuracy derived from MBES systems indicates similar or slightly better accuracy from the MBES systems (Dix et al., 2012;Ernstsen et al., 2006). Dix et al. (2012) determined the vertical accuracy of MBES data by testing the system on different objects and in different environments and found the vertical E RMS to be ±4 cm. Furthermore, they tested a lidar system on the same objects and found a similar vertical E RMS of ±4 cm. The vertical E RMS of ±4.1 cm found in this study is very close to both the MBES accuracy and lidar accuracy determined by Dix et al. (2012). Another study by Ernstsen et al. (2006) determined the vertical precision of a high-resolution shallowwater MBES system based on seven measurements of a shipwreck from a single survey carried out in similar settings as the present study, namely in the main tidal channel in the tidal inlet just north of the inlet investigated in this study. They found the vertical precision to be ±2 cm, which is slightly better than the vertical precision of ±3.8 cm (frame) and ±7.6 cm (block) found in this study. Overall, accuracy and precision are within the scale of sub-decimetres for both topobathymetric lidar and MBES systems, which enables the mapping of tidal basins with full coverage and with comparable quality.
Due to technical and logistical reasons, the data validation and the actual survey were carried out on different days and in different locations. Based on this, it is a fair question to ask whether the determined quality actually represents the quality of the data within the study site. Differences between the determined data quality at the validation sites and the data quality at the study site may arise from (1) different environmental conditions on the two surveying days and/or (2) different environments at the validation sites compared to the study site.
The environmental conditions were similar on the two surveying days (as mentioned in Sect. 3.2), meaning that the different days are not affecting the representation of the data quality within the study site.
The environmental differences between validation site 2 and the study site include the presence of up to 0.2-0.3 m waves in the main channel next to the study site. The waves introduce a source of error, because the proposed water surface detection method is not modelling the waves. The precision of the seabed points within the study site are therefore expected to be worse than the ±3.8 cm precision determined at validation site 2.
The water clarity/turbidity impacts the accuracy of the lidar data negatively, due to scattering on particles in the water column, which causes the laser beam to spread (Kunz et al., 1992;Niemeyer and Soergel, 2013). Moreover, part of the light is reflected in the direction of the receiver, and such return signals can be difficult to distinguish from the seabed return (Kunz et al., 1992). The turbidity was measured at validation site 2 and in the flood channel close to the study site during the 19 April survey by collecting water samples and subsequently analysing the samples for suspended sediment concentration (SSC) and organic matter content (OMC). The analyses showed that the average SSC was higher in the flood channel (17.2 mg kg −1 ) than in the river (10.2 mg kg −1 ). In contrast, the average OMC was lower in the flood channel (25.5 %) than in the river (40.0 %). These observations indicate that (1) the underwater precision is assessed in a location with higher turbidity than the environment within the study site; therefore, the turbidity cannot be a cause of lower precision in the study site, and (2) the penetration depth seems to be controlled by the OMC rather than by the SSC. This is new knowledge, since no previous studies (from what we know) have investigated the relative effect of organic matter as opposed to inorganic matter on the laser beam penetration depth. However, in order to determine the relationship with statistical confidence, a more comprehensive study is needed, involving measurements of penetration depth at different SSCs and OMCs, and without disturbance from other environmental parameters.

Spatial variations of topobathymetric lidar data quality
The quality of spatial datasets is often provided as single values, such as ±8.1 cm for the vertical accuracy in this case, and then the determined value represents the accuracy/precision of the whole dataset. However, in reality the value is only a measure of the local quality at the location where the assessment is conducted. The quality of the dataset varies spatially, and one way to illustrate this is to extract the maximum vertical difference between the lidar points of the processed point cloud within every 0.5 m × 0.5 m cell throughout the study site (Fig. 12). In flat areas, without multiple return signals, this shows the spatially varying precision of the dataset. There are large differences on Fanø, which is expected due to vegetation causing multiple lidar returns from both the vegetation canopy and from the bare ground.
In contrast, the differences on the very gently sloping, nonvegetated tidal flat are up to 10 cm, and there is no simple and natural reason for this variation. A range of factors contribute to the observed variations: -The laser beam incidence angle is determined by a combination of the scan angle, the water surface angle and the terrain slope. The shape of the footprint is stretched with larger incidence angles, and this effect can cause pulse timing errors in the detected signal, which lead to a decreasing vertical accuracy (Baltsavias, 1999). The error associated with larger scan angles is generally causing the outer beams, toward the swath edges, to attain a lower accuracy (Guenther, 2007). This is a reason for the observed variations along the swath edges (Fig. 12). Terrain slopes have the same effect of decreasing the vertical accuracy due to the footprint stretching. The measured elevation tends to be biased toward the shallowest point of the slope within the laser beam (Guenther, 2007). However, the influence of slope is not crucial in the Knudedyb tidal inlet system, since it is generally a very flat area.
-Areas covered by more than a single swath tend to show more vertical variation in the lidar point measurements, which means that there is a vertical bias between overlapping swaths. This can be caused by variance/error in the GPS measurements and/or IMU errors (Huising and Gomes Pereira, 1998). The vertical bias between swaths has been observed in the point cloud to be up to 5 cm, but it is varying throughout the study site. In most environments, a bias of 5 cm would be unnoticeable, but because of the large and very flat parts of the Knudedyb tidal inlet system, even a small bias becomes readily evident.
-The water depth is a factor, since the accuracy and precision are expected to be lower as the laser beam penetrates deeper into the water column (Kunz et al., 1992). The laser beam footprint is diverging as it moves through the water column, resulting in a larger footprint on the seabed. The elevation of the detected point is thus derived from the measurement on a larger area on the seabed, which will decrease the vertical accuracy, as well as decrease the capability of detecting small objects. With this in mind, the higher precision at the frame compared to the block is opposite of what would be expected, since the frame is below water and the block is on land. In this case, other factors, such as overlapping swaths and/or scan angle deviations, have more influence on the precision than the water depth. Also, it should be remembered that the frame surface was close to the water surface, and the effect of the water depth on the precision would most likely be more evident if it was located in deeper water.
Additional factors, beside the ones mentioned above, may influence the quality of lidar datasets. For instance, a dense vegetation cover of the seabed or breaking waves that makes the laser detection of the seabed almost impossible. However, these factors do not have a great influence in the studied part of the Knudedyb tidal inlet system, and thus they are not further elaborated.

Evaluation of the morphological classification
The morphological classification presented in this study is based on the studied section of the Knudedyb tidal inlet system. The overall concept of using tidal range, slope and variations of the elevation at different spatial scales proves to be a reliable method for delineating the morphological features in this tidal environment. The concept, however, can be applied in other environments. The specific thresholds in the classification determined in this study may deviate in other areas. Morphological features of different sizes require steps of other spatial scales in the neighbourhood analyses to produce a successful classification. In the future, the classification method will be improved by implementing an objective method for determining the scales, which can make it applicable in areas with different morphological characteristics. Such an objective scale determination method is presented by Ismail et al. (2015), who determined the scales based on the variance of the DEM at progressively larger window sizes. In this way, the sizes of the morphological features are determining the scales for the classification.
5.6 Using topobathymetric lidar data to map morphology in a highly dynamic tidal environment The study demonstrates the capability of green topobathymetric lidar to resolve fine-scale features, while covering a broad-scale tidal inlet system. Collecting topobathymetric li-dar data with a high point density of 20 points per m 2 on average enables detailed seamless mapping of large tidal environments, and the lidar data have further proved to maintain a high accuracy. The combined characteristics of mapping with high resolution and high accuracy in a traditionally challenging environment provide many potential applications, such as mapping for purposes of spatial planning and management, safety of navigation, nature conservation, or morphological classification, as demonstrated in this study. The developed lidar data processing method is tailored to a morphological analysis application. The best representation of the morphology is mapped by gridding the average value of the lidar points into a DEM with a 0.5 × 0.5 resolution. Other applications would require different gridding techniques. For instance, hydrographers, who are generally interested in mapping for navigational safety, would use the shallowest point for gridding. However, the overall method for processing the point cloud can be used regardless of the application. Only the last and least-challenging/time-consuming step of gridding the point cloud into a DEM may vary depending on the application. Applying topobathymetric lidar data for morphological analyses in tidal environments enables a holistic approach of seamlessly merging marine and terrestrial morphologies in a single dataset. However, a combination of topobathymetric lidar and MBES data is required, in order to map the morphology of tidal environments in full coverage. The comparable quality and resolution of lidar and MBES data gives a potential to map broad-scale tidal environments, such as the Wadden Sea, in full coverage and with high resolution and high accuracy.

Conclusions
A method was developed for processing raw topobathymetric lidar data into a DEM with seamless coverage across the land-water transition zone. The method relies on basic principles, and the entire processing method is described with a high level of detail, which makes it transparent and easy to implement for future studies. Specifically, a new procedure was developed for water surface detection in a tidal environment utilizing automatic water level determination solely based on green lidar data. The water surface detection method presented in this work did not take into account the variation in wave heights and surface slopes, which therefore constitutes a challenge to be addressed in future studies.
The vertical accuracy of the lidar data was determined by object detection of a cement block on land to ±8.1 cm with a 95 % confidence level. The vertical precision was determined at the cement block to ±7.6 cm, and ±3.8 cm at a steel frame placed just below the water surface. The horizontal mean error was determined at the block to ±10.4 cm. Overall, vertical and horizontal precision are within sub-decimetre scale.