Topo-Bathymetric LiDAR for Monitoring River Morphodynamics and Instream Habitats—A Case Study at the Pielach River

Abstract

Airborne LiDAR Bathymetry (ALB) has been rapidly evolving in recent years and now allows fluvial topography to be mapped in high resolution (>20 points/m²) and height accuracy (<10 cm) for both the aquatic and the riparian area. This article presents methods for enhanced modeling and monitoring of instream meso- and microhabitats based on multitemporal data acquisition. This is demonstrated for a near natural reach of the Pielach River, with data acquired from April 2013 to October 2014, covering two flood events. In comparison with topographic laser scanning, ALB requires a number of specific processing steps. We present, firstly, a novel approach for modeling the water surface in the case of sparse water surface echoes and, secondly, a strategy for improved filtering and modeling of the Digital Terrain Model of the Watercourse (DTM-W). Based on the multitemporal DTM-W we discuss the massive changes of the fluvial topography exhibiting deposition/erosion of 10³ m³ caused by the 30-years flood event in May 2014. Furthermore, for the first time, such a high-resolution data source is used for monitoring of hydro-morphological units (mesohabitat scale) including the consequences for the target fish species nase (Chondrostoma nasus, microhabitat scale). The flood events caused a spatial displacement of the hydro-morphological units but did not effect their overall frequency distribution, which is considered an important habitat feature as it documents resilience against disturbances.

1. Introduction

Remote sensing has been recognized as a viable method for habitat monitoring, i.e., the repeated assessment of the state and its changes [1,2]. A periodic assessment in a six year cycle of the state of habitats is required by the European Union’s Fanua-Flora-Habitat (FFH) directive [3,4], including aquatic habitats. Concerning rivers also other European directives, the water framework directive [5] and the flood directive [6], formulate requirements on monitoring.

River morphodynamics are important driving agents for renewing instream habitats and thus improving the quality of the aquatic environment. However, a major limitation of instream habitat studies so far is that numerical approaches used to evaluate residual flow and/or river training measures are based on the assumption of stable morphological boundaries [7,8]. In most cases sedimentological and morphodynamic processes tend to be neglected. In some smaller studies the importance of river morphology and morphodynamic processes on instream habitats is discussed, e.g., [9]. In addition, as noted by [10], the dynamics in channel morphology influence the occurrence and diversity of functional habitats (mesohabitats) in semi-natural and degraded rivers in many ways.

Morphodynamic processes in instream habitat evaluation can be considered in two different ways. On the one hand, sediment transport models with integrative adaption of morphological changes may be used to consider the impacts of channel forming flows [11]. On the other hand, repeated topographical surveys allow detailed studies of morphological (habitat) changes after periods of high flows, e.g., [12]. Within our work, we concentrate on the second approach of repeated surveying and change detection based on Digital Elevation Model of Differences (DoD).

Beside the importance of channel forming flows for the aquatic environment, the related inundation of floodplains has to be seen as a natural hazard especially in densely populated areas (e.g., determination of endangered zones). To cope with flood hazards, standardized methodologies are required and partially established by water management authorities. In this context, the European Commission released the Floods Directive [6] which requires flood risk maps for all river basins and sub-basins with significant potential of risk of flooding for Europe [13]. In particular, the cost effectiveness of mitigation measures in comparison to the flood hazard has to be studied based on detailed risk maps [14,15]. Examples of flood risk analysis are available on various scales like the municipal level [16] or the catchment level [17]. Depending on the scale of the hazard and risk analysis, the complexity of evaluation tools (e.g., models) range from simple interpolation methods to spatially detailed models using depth-averaged hydraulics on unstructured grids [13,18]. However, the latter can be regarded as state of the art where usually infrared LiDAR for the floodplains and banks above water level are fused with (cross sectional) terrestrial survey or echosounding for the wetted channel. While surveying with two different methods is feasible, it also raises the question of the common reference frame for both campaigns, and one technique for observing both the dry and the wet area, would be an advantage.

Surveying of riparian areas and water bodies by remote sensing is challenging. Vegetation is often overhanging the river channel, but also two media need to be considered for the ray propagation: air and water. Sensing of water color with multi- or hyper-spectral remote sensing is one field of study, providing optical properties and near shore bathymetry [19,20,21,22,23] or water quality [24]. Airborne LiDAR, exploiting near infrared pulses for acquisition of point clouds over land surfaces, has proved to be a suitable method for the assessment of the 3-dimensional structure of vegetation [25,26,27] and terrain surface [28]. The LiDAR principle is to emit a laser pulse in a known, controlled direction and measure the time from emission until receiving of its echo, scattered back from surfaces within the instantaneous field of view. Direct georeferencing provides the position and orientation of the measurement platform. Together, this is used for 3-dimensional location of the echoes, thus providing a georeferenced 3D point cloud [29].

Airborne laser bathymetry (ALB), also referred to as airborne laser hydrography (ALH) or LiDAR bathymetry, is a technique for measuring the depths of relatively shallow, coastal waters from the air using such a scanning, pulsed laser beam [30]. With the recent advent of topo-bathymetric LiDAR, simultaneous monitoring of land surfaces and water bodies became possible. A green laser is used for the water body, as this wavelength propagates into the water and is reflected from the bottom surface, while an infrared laser is used for the detection of the water surface (air/water interface), the land surface, and the vegetation above it. Knowledge of the water surface is required as the green light is refracted at this interface. The point cloud from the river bottom surface, therefore, can only be computed after the river top surface is modeled. Using airborne acquisition, thus moderate flying heights (few hundreds up to approx 2000 m), and high pulse repetition frequency, high resolution of multiple points per square meter can be achieved.

The green light, however, does not only propagate into the water, but it is also reflected at the land surface and vegetation. Thus, also LiDAR systems for topo-bathymetric surveys using the green wavelength only are built. An obvious advantage is that only one LiDAR is required, leading to simpler measurement devices.

Using only the green wavelength, on the other hand, also poses challenges for the determination of the point cloud in the water area as especially smooth water surfaces lead to very few reflections from the water surface itself. As detailed in [31] using only the green signal poses a difficult challenge in accurately and reliably determining the air/water interface for each pulse, whereas [32] demonstrated the feasibility to determine the water surface based on statistical methods from a sufficient number of green near water surface returns.

Furthermore, as the riparian area itself has a complex morphology, the modeling of the land surface is additionally challenging in areas of natural rivers. It is characterized by steep banks and the alluvial forest with complex understory.

The objective of this paper is to present and evaluate new methods (i) for determination of the river water surface using only few reflections available from the green LiDAR; (ii) for modeling the topography of the river bed and the riparian area; (iii) derive spatially variable DoD error estimates for significantly distinguishing between real changes and noise; (iv) to study the dynamic river morphology from the aforementioned models and verify the findings by field visits; and (v) to demonstrate how these models can be used for modeling aquatic habitats.

We will first review related work, then present the study area and the acquired data (Section 2). The new methods are detailed in Section 3 and evaluated and discussed in Section 4.

Related Work

With the increasing availability of cost-effective methods for the acquisition of topographic data at high spatial resolution and extent Digital Terrain Models (DTM) are widely used for understanding geomorphic processes in general [33], and river channel systems in particular [34]. In this context the surface of interest includes the topography of the river bed within the inundated area and the bare ground adjacent to the wetted channel. This specific geometry is referred to as the Digital Terrain Model of the Watercourse (DTM-W) [35] but the more common term DTM is used synonymously throughout this paper.

For capturing riverine environments the entire range of modern acquisition techniques is applied. Terrestrial surveys using Global Navigation Satellite System (GNSS) rovers or total stations are, e.g., reported by [12,36,37]. Although modern instruments allow rapid acquisition, ground based surveys are restricted to smaller rivers and the density of river bed points is low and inhomogeneous. Topographic LiDAR can directly be used for capturing the riparian area [38,39,40] but the near-infrared laser signal is absorbed at the water surface. Measuring water depths based on RGB and/or hyperspectral images became feasible by estimating the correlation between water depth and image color [21,22] and is frequently fused with LiDAR and/or ground based data [41,42,43,44]. Passive optical methods require terrestrial calibration measures and are applicable to clear rivers with water depths below the Secchi depth. Acoustic methods, e.g., [45], are not reviewed in detail here as their application is not feasible within the unnavigable study reach due to occasional riffles and shallow shoals. In most of the aforementioned studies the high resolution geomorphic analysis is performed using DEM of Differences (DoD) models [46] which is also used in our study.

LiDAR bathymetry was used in coastal zones for mapping and classification of the habitats (benthic substrata and macroalgal communities) [47]. They used the LADS Mk II instrument and underwater videos for training data and ground truth. Derivatives of the bottom surface (curvature, etc.) were employed in an automatic classification. In a comparable approach, [48] used the SHOALS bathymetric LiDAR to study the complexity of benthic habitats. A 4 m grid of the bottom surface was used to derive variables of substrate rugosity (complexity) and correlate those to in-situ measurements. In [49] it is shown that the additional use of LiDAR bathymetric bottom surface data from the Hawk Eye MKII system, again in the form of terrain variables, improves the classification of coastal habitats compared to using only passive optical imagery. Surveying of the coastal zone after storm events is described in [50], applying the same instrument, the Riegl VQ-820-G, also used in our study. A comparison of different instruments over one site is presented in [51].

In [52] airborne LiDAR bathymetry for the wet area and standard airborne topographic LiDAR from two different campaigns are combined for modeling the riverine area. A numerical hydrodynamic 2D model is used to derive mesohabitats (pool, riffle, glide) and microhabitats on the basis of Froude numbers (esp spawning habitats). This is performed for one epoch and verification was performed by field investigations, going along the river with a boat. This reference also includes an accuracy assessment, indicating a mean elevation difference of 15–25 cm between terrestrial and LiDAR bathymetric measurements and a standard deviation of 11–35 cm. Also [53] presents an accuracy analysis of a topo-bathymetric small-footprint LiDAR by comparison to terrestrial measurements, confirming the quality measures of the previous reference. In [54] the SHOALS instrument is investigated for a river, deriving a bias of 10–20 cm and a standard deviation of 15–40 cm. The Aquarius sensor was tested by [55] both in coastal and river environments. For the latter, RMSE values in the range of 10–25 cm are reported comparing the bathymetric LiDAR with RTK GPS reference points. The effects of bathymetric LiDAR errors on the flow properties derived from a multi-dimensional fluid dynamics model were assessed in [56] for the Experimental Advanced Airborne Research LiDAR (EAARL). The authors conclude that ALB can map channel topography with sufficient accuracy to support numerical flow models.

This literature review shows that monitoring of instream habitats via repeat surveys has not been performed by airborne topo-bathymetric LiDAR so far, especially not at simultaneously high temporal and geometric resolution.

2. Study Area and Data Sets

Multitemporal topo-bathymetric surveys were conducted at the pre-alpine Pielach River, a medium-sized right side tributary of the Danube (Figure 1). The modeling reach Neubacher Au is located about 100 km west of Vienna downstream of the village of Loosdorf in Lower Austria (48 °12 ’50 ”N, 15 °22 ’30 ”E; WGS 84). The Pielach River originates at an elevation of about 1000 m above sea level and flows 68 km northwards down to the Danube (mouth at 208m a.s.l.), reaching stream order 4 according to Strahler [57]. The catchment size of this pre-alpine river is approximately 590 km² with a recorded mean annual precipitation of 875 mm. The mean discharge for the study site is about 7 m³·s⁻¹ considering the gauging stations Hofstetten and the main upstream tributary Sierningbach (cf. Table 1, Figure 2). The river morphology is classified as widening with a variability of gravel bars ranging from mid-channel bars to point bars. The river is classified as riffle-pool type [58] with an average gradient in the study reach of 0.39% . Bed-load sediments are dominated by coarse gravel (2–6.3 cm) along the active channel and bars, and partially cohesive sediments in areas of bank erosion leading to steep bank slopes. The longitudinal continuum for fish migration in the Pielach River is disrupted by weirs built for hydropower use and engineering measures, but it is possible for fish to enter the first two kilometers from the Danube (Figure 1) [57]. Nevertheless, the Pielach River has retained some of its natural self-forming morphological characteristics in the studied alluvial forests of the Neubacher Au, namely periodically inundated side arms, dynamic gravel bars, large woody debris, small oxbows, inundation areas and various succession stages of floodplain forests [59]. The Neubaucher Au is part of the Natura2000 conservation area Niederösterreichische Alpenvorlandflüsse (Area code: AT1219000). The fish community in the studied reach is dominated by nase (Chondrostoma nasus, 20% ) and barbel (Barbus barbus, 20% ). Moreover, other abundant species out of a total of 29 species are Alburnoides bipunctatus, Barbatula barbatula (L.), Cottus gobio (L.), Gobio gobio (L.), Hucho hucho (L.), Leuciscus cephalus (L.), Leuciscus leuciscus (L.), Phoxinus phoxinus (L.) and Salmo trutta (L.) [60].

Figure 1. Neubacher Au, Pielach River, Lower Austria; (a) Flight block Loosdorf, shaded relief of digital surface model superimposed with color-coded elevation map and additional vector data (flight trajectory, river axis, freshwater ponds); (b) Overview map of Austria with marked location of study area; (c) Pielach River catchment and location of gauging stations; (d) Study area Neubacher Au, mean flow and annual flood area highlighted; (e–g) Terrestrial photos; (e) meander bend; gravel bank and bar, steep bank, backwater; (f) main channel with gravel bar, deadwood, and alluvial forest; (g) meander bend from top of the steep bank during flood event on 16 May 2014.

Figure 1. Neubacher Au, Pielach River, Lower Austria; (a) Flight block Loosdorf, shaded relief of digital surface model superimposed with color-coded elevation map and additional vector data (flight trajectory, river axis, freshwater ponds); (b) Overview map of Austria with marked location of study area; (c) Pielach River catchment and location of gauging stations; (d) Study area Neubacher Au, mean flow and annual flood area highlighted; (e–g) Terrestrial photos; (e) meander bend; gravel bank and bar, steep bank, backwater; (f) main channel with gravel bar, deadwood, and alluvial forest; (g) meander bend from top of the steep bank during flood event on 16 May 2014.

Table 1. Annuality discharge statistics, MQ = mean discharge, HQ_n = flood discharge for recurrence interval of n years.

Annuality	MQ	HQ1	HQ5	HQ10	HQ30	HQ100
Discharge [m3 ·s−1]	7.16	136.6	203	261	347	440

Table 1. Annuality discharge statistics, MQ = mean discharge, HQ_n = flood discharge for recurrence interval of n years.

Annuality	MQ	HQ1	HQ5	HQ10	HQ30	HQ100
Discharge [m3 ·s−1]	7.16	136.6	203	261	347	440

Figure 2. Discharge of Pielach River at Loosdorf 2013–2014 obtained from combining the discharges at gauging stations Hofstetten (Pielach) and Sierningbach; vertical red lines: flight dates; horizontal black lines: characteristic discharges.

Figure 2. Discharge of Pielach River at Loosdorf 2013–2014 obtained from combining the discharges at gauging stations Hofstetten (Pielach) and Sierningbach; vertical red lines: flight dates; horizontal black lines: characteristic discharges.

Since March 2013 the study area has been repeatedly captured with the topo-bathymetric sensors VQ-820-G [61] and—since October 2014—VQ-880-G [61] primarily for the purpose of sensor calibration. The scanner is mounted in the nose pod of a Diamond DA42 light aircraft flying at 600 m above ground level which corresponds to the scanners extended nominal ocular hazard distance (eNOHD). At this flying height the diameter of the laser footprint on the ground is 60 cm. To achieve a high point density for detailed mapping of both the river bed and the riparian area, the scanners were operated with the highest possible pulse repetition rate of 510 kHz (VQ-820-G, effective measurement rate: 220 kHz) and 550 kHz (VQ-880-G, effective measurement rate: 550 kHz), respectively. Both scanners record the echo waveforms and perform online waveform processing resulting in additional attributes (amplitude, pulse shape deviation, reflectance) per echo [62].

Within the April 2013–October 2014 period two flood events occurred at the Pielach River: An annual flood (HQ₁) in June 2013 and a major event in May 2014 with a discharge beyond a 30-years flood (HQ₃₀). Both events had an impact on the fluvial topography of the channel and the riparian area. From the existing flight data, four epochs have been chosen for this study. The initial dataset (April 2013) was captured before the June 2013 floodHQ₁, followed by two intermediate epochs between the June 2013 and May 2014HQ₁ and the HQ₃₀ events (October 2013 and February 2014). The last flight in October 2014, finally, was carried out after the HQ₃₀ flood. Table 2 summarizes the main data capturing properties. Whereas the flights in winter and spring (April 2013 and February 2014) were conducted under leaf-off conditions enabling a good penetration of the vegetation, the defoliation was still in progress for both October flights. This especially applies to the complex and dense understory within the alluvial forest and affects the achievable DTM quality. The mean last echo point density is in the range of 20 points per square meter, with a slightly higher density in February 2014 due to additional flight strips compared to the standard flight planning. In general, the entire project area was captured with four longitudinal and three cross strips (cf. Figure 1) which allows on-the-job calibration of the mounting calibration (boresight angles).

Preprocessing of the scan data included the following steps: (i) derivation of the raw 3D point cloud via direct georeferencing based on the LiDAR sensor model [63]; (ii) strip adjustment/boresight alignment [64,65]; (iii) transformation of the WGS84 coordinates to the target spatial reference system (ETRS89, UTM zone 33, ellipsoidal heights); (iv) assessment of the strip fitting accuracy and point density [66]; and (v) datum transformation of the October 2013, February 2014, and October 2014 flight blocks using the initial April 2013 epoch as reference. Only corresponding smooth and tilted planes were used to estimate a best fitting rigid body transformation (3 shifts and 3 rotations) for each block via Least Squares Matching [66].

Preprocessing of the LiDAR point cloud was performed with the Riegl ALS software suite RiProcess [67] (waveform processing, range determination, direct georeferencing, strip adjustment) and the scientific laser scanning software OPALS [68] (quality control, datum transformation). The precision and accuracy of the flight block registration are documented in Table 2. The standard deviation of the height discrepancies in overlapping strip areas expressed as the median of absolute differences (sigma_MAD ) is less than 2 cm for all flight dates (see Table 2, column: precision) as is the absolute fitting accuracy of the individual flight blocks compared to the reference epoch (April 2013, Table 2, column: accuracy). These measures are below the nominal ranging precision of the scanners (25 mm) as stated by the manufacturer. It is noted that a good conformance of the flight block orientations is a precondition for reliable detection of topographic changes as a result of fluvial activity.

Table 2. Summary of data capturing.

Flight Date	Sensor	Q [m3 ·s−1]	Point Density [points/m2]	Precision [m]	Accuracy [m]	Foliage
15 April 2013	VQ-820-G	12.2	12	0.016	0.000	leaf-off
26 October 2014	VQ-820-G	2.2	18	0.019	0.015	leaf-on
21 February 2014	VQ-820-G	3.4	25	0.018	0.017	leaf-off
16 October 2014	VQ-880-G	3.7	22	0.019	0.019	leaf-on

Table 2. Summary of data capturing.

Flight Date	Sensor	Q [m3 ·s−1]	Point Density [points/m2]	Precision [m]	Accuracy [m]	Foliage
15 April 2013	VQ-820-G	12.2	12	0.016	0.000	leaf-off
26 October 2014	VQ-820-G	2.2	18	0.019	0.015	leaf-on
21 February 2014	VQ-820-G	3.4	25	0.018	0.017	leaf-off
16 October 2014	VQ-880-G	3.7	22	0.019	0.019	leaf-on

3. Methods

In this Section, specific methods for processing ALB point clouds are presented. In general, the ALB workflow is very similar to the processing chain of traditional topographic LiDAR data described in, e.g., [68]. One crucial difference is the correction of the water echoes due to the refraction of the laser beam at the air/water-interface. A novel approach for modeling the water surface and classification of water echoes is presented. Furthermore, the complex structure of the river bed and riparian area poses specific problems for DTM filtering and requires specialized strategies which will be described. Finally, our approach for deriving spatially variable DoD error estimates based on the methods of [46,69] are outlined. First, an overview on the assembly of methods is given.

3.1. Modeling Terrain and River Surfaces

Figure 3 gives an overview of the methods applied to derive the Digital Water surface Model (DWM), the Digital Terrain Model of theWatercourse (DTM-W), and finally a “DEM of Differences” (DoD) providing reliably determined deposition and erosion areas between epochs. Starting from the full waveform LiDAR point cloud, preprocessing is applied as described above. Water surface modeling requires additionally the river axis and is described in Section 3.2. Points below this surface are subject to the refraction correction, resulting not only in the corrected position of those points, but also in an additional attribute per point (wet/dry-indicator). The wet points are further classified into river surface, river body, and river bottom points, described in Section 3.3, utilizing full waveform features and water depth information. The river bottom echoes, the dry last echo points, and automatically derived break lines are used for deriving a terrain model describing the watercourse (Section 3.4). Also in this step full waveform features are exploited. These models are input for the pairwise DoD computation. Spatially variable DTM error measures are calculated for each epoch and propagated into the DoD for separating reliable changes from noise (Section 3.5). Finally, based on a two-dimensional depth-averaged numerical model [70] the abiotic characteristics (flow velocity, water depth and bottom shear stress) are simulated and analyzed for different epochs and different discharges (low flow, mean flow and the magnitude of an annual flood event). Considering the variability in discharge dependent changes of flow variables, a mesohabitat evaluation approach (MEM) was selected to determine the impact of morphological changes on the habitat distribution [9]. The MEM approach enables a differentiation of six hydro-morphological units, namely (1) riffle-; (2) fast run-; (3) run-; (4) pool-; (5) backwater-; and (6) shallow water-habitats.

Figure 3. Overview on modeling.

Figure 3. Overview on modeling.

3.2. Water Surface Modeling

For modeling the river bottom surface, the water top surface needs to be known as well. At this interface the incident laser beam is refracted, thus requiring elevation and gradient of this surface for modeling the beam propagation. In the case echoes are recorded from both the water surface and the river bottom, a classification of these echoes needs to be performed. This is hindered by spurious false detections and by echoes from vegetation. Vegetation occurs above, below, and horizontally near to the water surface.

In [32] it was demonstrated that a reliable reconstruction of the air/water-interface based on the green laser signal is feasible using statistical methods. The vertical location of the echoes from the air/water-interface show a distribution not centered on the surface level, but systematically below it.

If, however, the water surface generates only few detectable echoes, this results in a sparse water surface point set. Therefore, the basis for application of statistical methods is not given. For ALB acquired under low flow conditions, a sparse water surface representation indeed occurs. The February 2014 data acquisition, for instance, exhibited a water surface echoes coverage of only 25% counting cells of 1 m² size with at least a single surface echo as valid. The stretch of the data voids along the river axis amounted to more than 20 m.

In this case, the smooth water surface results in a high amount of specular reflections and, consequently, surface echo drop outs. Thus, the required higher point density is only available in areas with rough water surface. As no major gradient changes can be expected in smooth areas it is, therefore, still possible to derive a continuous water surface model by axis oriented interpolation.

The suggested method targets an efficient reconstruction of the air/water-interface in the case of a very low water surface echo density. It is a combination of manual and automatic steps. In general, it is based on the concept for generating a river bed model from sparsely captured cross sections [71]. The procedure requires the 3D LiDAR point cloud and the 2D river axis as input.

In a second step, the point cloud is extracted in cross sectional slices perpendicular to the river axis . The operator, then, manually defines the water level height (assuming a constant water level within the section) and the lateral extent (slightly cutting into the terrain) in a simple graphical editor.

The definition of water level height and/or lateral extent is repeated for all sections along the axis. The editor assists the operator by automatically extrapolating the water surface height of the consecutive section and by showing already measured sections in a ground plan and longitudinal section, cf. Figure 4. A continuous water surface model is derived in a final interpolation step by employing an axis oriented cross section densification approach [71]. The final result of this processing step is a Digital Water surface Model (DWM) which can be stored as a regular grid, e.g., grid width 0.5 m.

Figure 4. Semi-automatic derivation of water surface model; (a) cross section with LiDAR echoes colored by reflectance and manually defined water level and extent (red line); (b) longitudinal section with measured water levels (green circles) and/or extents (gray dots); (c) plan view, measured sections marked in green/orange.

Figure 4. Semi-automatic derivation of water surface model; (a) cross section with LiDAR echoes colored by reflectance and manually defined water level and extent (red line); (b) longitudinal section with measured water levels (green circles) and/or extents (gray dots); (c) plan view, measured sections marked in green/orange.

3.3. Water Body Classification

The DWM serves as the basis for range and refraction correction of the water echoes based on Snells law [72]. All points on or below this surface are classified as water (i.e., wet points). The strategy for classifying the wet points further into (i) surface; (ii) body; and (iii) bottom points (Figure 3) relies on full waveform features and features computed from the local neighborhood, as well as spatial reasoning, thus general point cloud processing [29].

As a consequence of the interaction of the green laser signal with the medium water, the echoes either stem from a reflection at or near the water surface, from scattering at water particles (volume backscatter) or from the river bed. In and next to the river bed, the correct classification is difficult, especially in deep water areas. There, the progressive attenuation of the laser signal in the water column causes the density of the river ground points to drop sharply. In the case of very low river bottom point density, volume backscatter can easily be interpreted as ground in automatic DTM filtering algorithms relying on coordinates alone. However, echoes from volume backscatter typically exhibit a low reflectance (cf. Figure 5) and show a sparse spatial occurrence. Finally, using the water surface model (Section 3.2), the depth of each echo below the air/water-interface can be computed and aids classification.

In an initial step a first set of surface echoes is identified by their low reflectance and small water depth. Occasional near water surface echoes with higher reflectance values are included in a second step if the majority of points in a spherical 3D neighborhood belongs to the water surface class. This will not affect river bed points of the littoral area as the ground point density is typically high in shallow water. The volume backscatter points are identified in a third step by their low reflectance and sparse point spacing. To avoid misclassification of river bed points at the maximum penetration depth—where the point density is generally low—volume backscatter echoes are restricted to a certain water depth. This strategy requires three types of thresholds related to: (a) water depth; (b) object reflectance; and (c) neighborhood definition. The multiple depth related thresholds can be scaled according to the maximum penetration depth of the sensor depending on the Secchi depth [73] of the water body. The same applies for the dimension of the search spheres and/or cylinders (c), which are scaled depending on the overall point density (i.e., mission parameters). The reflectance threshold mainly depends on the environmental condition (attenuation within the atmosphere and the water column). Mission independence can be achieved by performing radiometric calibration using external reference targets [74].

Figure 5. Classification of water echoes based on echo reflectance, distance from water surface and spatial distribution; principal sketch.

Figure 5. Classification of water echoes based on echo reflectance, distance from water surface and spatial distribution; principal sketch.

3.4. Terrain Modeling

The fluvial topography of near natural rivers is typically very complex. The river area is often overgrown with riverside vegetation, the embankments exhibit steep slopes and even overhangs, and the alluvial forest is covered by dense understory. Obtaining the ground surface echoes from airborne laser scanning was studied extensively, with comparative evaluations given, e.g., by [75,76]. In these surveys step edges (overhangs), low vegetation, and areas of missing ground points are identified as problems in ground detection. However, these problems are abundant in areas of natural rivers. The suggested method tackles these problems by (i) explicit modeling of the break lines; (ii) exploitation of full waveform echo shape information for identifying low vegetation; and (iii) using only the river bottom points identified as above.

The suggested method for coping with the immediate gradient changes is first modeling the 3D break lines following the approach of [77]. This method identifies potential segments of break lines in the unfiltered point cloud by curvature analysis. Longer break lines are established by tracing the line starting at such a start segment. The break line is modeled locally as a pair of intersecting planes, each plane being supported by measured points on either side. This modeling is performed in the 3D point cloud, thus no implicit ground projection is performed, and therefore also the rim above a vertical face can be modeled.

An additional ground point filtering problem arises from the dense understory in the alluvial forest. It is addressed by an adaptation of the strategy proposed by [78] based on full waveform echo attributes. While topo-bathymetric LiDAR sensors often record the full waveform, it is not necessarily that the echo width used in the approach of [78] is provided. The VQ-820-G sensor, e.g., stores for each point next to its location the so-called pulse shape deviation . This parameter is determined during online waveform processing [79]. It is used for eliminating unreliably determined echoes before the actual DTM filtering.

The ground filtering is performed using the hierarchic robust interpolation approach [80] implemented in the program system SCOP++ [81]. In this approach, the terrain surface is determined iteratively. Points obtain higher weights, if they are more likely to represent the ground surface and lower weights, if they are in the vegetation. The break line points are not subject to this reweighting scheme, as they can be assumed to represent the ground. In a final step the DTM model is interpolated based on linear prediction, also referred to as simple kriging [82], using all topo-bathymetric LiDAR ground points and break lines as input.

3.5. DEM of Differences

The sediment budget (deposition and erosion rates) are estimated based on gridded DEM of Differences (DoD) models as:

z_DoD = z_new − z_old

(1)

with z_DoD = deposition/erosion value of a single DoD-cell, and z_new/z_old the respective DTM cell elevations of the studied epochs. The basic question if a certain DoD value represents a real change or just measurement noise requires the estimation of the DTM errors which are propagated into the DoD according to [83]:

$${\sigma }_{DoD}=\sqrt{{\sigma }_{DTM,new}^2+{\sigma }_{DTM,old}^2}$$

(2)

σ_{DT M,new/old} hereby denote the spatially variable DTM accuracy. In contrast to Wheaton et al. [46] who estimate the DTM error using a probabilistic approach by applying a Fuzzy Inference System (FIS) based on point quality, point density and surface slope as input variables, our approach for deriving spatially variable DTM error values follows the deterministic approach of Kraus et al. [69] based on the LiDAR point precision, the local point density and surface curvature. This method makes use of the fact that the LiDAR point density is generally higher than the DTM grid point density (20 points/m² vs. 4 points/m² in our case). As the DTM interpolation (linear prediction) is not forcing the DTM surface through the data points, the height deviations of the LiDAR ground points are used for the estimation of the point precision. The achievable DTM precision furthermore depends on the local LiDAR point density and the relief (the higher the point density and the smoother the surface, the higher the DTM precision).

The absolute DTM accuracy is further influenced by the variability in elevation of natural targets, which cannot be captured by airborne LiDAR, e.g., elevation changes within and between footprints. This value is termed σ_rough and mainly depends on land cover and land use. It is referred to as “uncertainty of definition” in photogrammetry. To achieve realistic DTM error estimates this variability is added to the calculated DTM precision.

$${\sigma }_{DTM,acc}=\sqrt{{\sigma }_{DTM,prec}^2+{\sigma }_{rough}^2}$$

(3)

with σ_DTM,acc representing the (absolute) DTM accuracy, σ_DTM,prec the DTM precision according to [69] and σ_rough the surface “roughness” of natural targets. In this paper the latter is assumed to be constant within homogeneous areas of land cover (i.e., river bed, forest, grassland, gravel, cf. Section 4.1).

For deciding if a certain DoD cell should be considered in the sediment budget we follow the statistical approach of Wheaton et al. [46]. A student’s t-score is derived for each cell by comparing the absolute DoD value and the local DoD error and the 95% confidence interval is used as the critical threshold for the derivation of a binary grid mask. Wheaton et al. furthermore propose to preserve additional depositional or erosional cells in the sediment budget based on their spatial coherence even if their t-score would mark them as noise. We do not apply this approach as it is only applicable for DTM data free of systematic errors. In Section 4 we will exemplary show that this does not hold for airborne LiDAR data flown at different phenological stages (leaf on, leaf off). We, therefore, stick to the more conservative estimate based on the 95% confidence interval.

3.6. Habitat Modeling

Habitat modeling was performed using the MEM approach [9] based on hydraulic parameters as input. The applied two-dimensional depth-averaged hydrodynamic-numerical (HN) model Hydro_as-2D [70] calculates hydraulics by applying a finite volume approach based on an unstructured Triangular Irregular Network (TIN). The convective flow of the 2D model is based on the upwind-scheme by Pironneau [84] and the discretization of time is done by an explicit Runge-Kutta method in second order. The simulations were performed with the parabolic eddy viscosity model and a constant turbulent viscosity coefficient of 0.6. The total computational time of the model was set to 9.25 h. The models performance showed that steady state conditions could be achieved in 6.25 h (low-flow) and 3.5 h (high-flow) for the various HN models of the study reach. The depth criterion was set to 0.01 m for wetting/drying computations. Roughness (Manning n-values) was calibrated for low-flow. For high-flow, existing calibrated roughness coefficients have been used. For the main channel Manning n-values vary from 0.043 (low-flow) to 0.033 (high-flow). For the inundation areas Manning n-values exhibit a range of 0.067 (dense floodplain vegetation) to 0.017 (sealed road in the overbank areas). The input discharge data were obtained from long-term recording gauging stations in the Pielach River catchment. In addition to gauging station Hofstetten, located 26 km upstream of the study reach, the hydrological data of the major tributary Sierningbach was additionally considered to obtain representative discharges for the modeling site. The following steady state discharges were used for the detailed morphological and habitat investigations: mean flow (7.16 m³·s⁻¹) and annual flood (136.6 m³·s⁻¹).

Mesohabitat description and quantification was carried out based on the Mesohabitat Evaluation Model (MEM). The conceptual MEM-model was developed by [9] and allows the differentiation of six different mesohabitats according to their abiotic characteristics. Three abiotic parameters (flow velocity, water depth, and bottom shear stress) are used in the MEM-classification in which the energetic potential of the mesohabitats due to bottom shear stress is one of the criteria for differentiation. Riffles and fast runs are classified as high energetic (high bottom shear stress), run and pool as moderate energetic and backwaters and/or shallow water habitats as low energetic hydro-morphological units [85].

4. Results and Discussion

In this section accuracy assessments of the LiDAR point cloud, the multitemporal DTMs, their morphodynamic changes, and the results of flood modeling as well as habitat modeling are presented and discussed.

4.1. Point Cloud Accuracy Assessment

For validating the accuracy of the captured topo-bathymetric point cloud a terrestrial survey was conducted in February 2014 (cf. Figure 6). More than 2000 points were measured in sub-centimeter accuracy with a total station (Leica TPS1200). The surveyed reference points were classified into five categories: (i) sealed road; (ii) river bed; (iii) alluvial forest; (iv) grassland; and (v) gravel. Those five categories represent all relevant land cover categories within the study area. The spatial distribution of the acquired points is shown in Figure 6. The figure additionally shows histograms of the height discrepancies between the LiDAR point cloud and the reference points. The numeric results of this accuracy assessment are presented in Table 3. For each reference point the elevation difference was calculated based on the median of the four nearest 3D-neighbors in the LiDAR point cloud.

The statistics in Table 3 summarize all parameters contributing to the error budget of both the LiDAR point cloud and the reference points. Reference points were measured on artificial, man-made surfaces (sealed road) and on natural surfaces (river bed, forest, grassland, and gravel). The discrepancies can be split up in bias (mean/median) and precision (standard deviation, sigma_MAD). In the following only the robust statistics (median, sigma_MAD) are used for the further interpretations and calculations.

Table 3. Distribution of height discrepancies between topo-bathymetric point cloud and reference points; last row: surface roughness of natural targets.

Parameter	Sealed Road	River Bed	Alluvial Forest	Grassland	Gravel Bank
Mean [m]	0.001	0.000	0.069	0.038	0.025
Median [m]	0.001	–0.006	0.055	0.031	0.022
Std.dev. [m]	0.011	0.040	0.058	0.042	0.037
sigmaMAD [m]	0.007	0.025	0.050	0.031	0.022
Surface roughness [m]	0.000	0.024	0.050	0.030	0.021

Table 3. Distribution of height discrepancies between topo-bathymetric point cloud and reference points; last row: surface roughness of natural targets.

Parameter	Sealed Road	River Bed	Alluvial Forest	Grassland	Gravel Bank
Mean [m]	0.001	0.000	0.069	0.038	0.025
Median [m]	0.001	–0.006	0.055	0.031	0.022
Std.dev. [m]	0.011	0.040	0.058	0.042	0.037
sigmaMAD [m]	0.007	0.025	0.050	0.031	0.022
Surface roughness [m]	0.000	0.024	0.050	0.030	0.021

Figure 6. Reference survey for validation of ALB point cloud accuracy.

Figure 6. Reference survey for validation of ALB point cloud accuracy.

The height deviations for reference points measured on sealed roads are unbiased and the sigma_MAD of 7 mm document a very good precision/accuracy of the topo-bathymetric point cloud on artificial surfaces, especially considering the medium size of the laser footprint of 60 cm. This, once again, underlines the ranging performance of the sensors, the proper calibration of the multi-sensor system and the absence of geodetic datum errors. All other reference points were measured on natural targets. The respective statistics are biased and show a lower precision compared to the sealed road results depending on the data category. According to Kraus [82], the height accuracy of a natural target can be computed as:

$${\sigma }_{z,nat}=\sqrt{{\sigma }_{target}^2+{\sigma }_{rough}^2}$$

(4)

with σ_z,nat as the height accuracy of a natural target, σ_target as the height accuracy of the target point constrained by the data acquisition method, and σ_rough as the surface roughness. Using sigma_MAD of the sealed road category as a reasonable estimate of the potential LiDAR height accuracy, the surface roughness can be calculated for all categories as:

$${\sigma }_{rough,ca{t}_n}=\sqrt{{\sigma }_{z,ca{t}_n}^2-{\sigma }_{z,sealed}^2}$$

(5)

with σ_rough,catn as the surface roughness for the category n, σ_z,catn as thesigma_MADaccording to Table 3 for category n, and σ_z,sealed as the sigma_MAD for the sealed road category. The final values for the surface roughness are documented in the last row of Table 3. For the DoD error estimation they are assumed constant for the respective area within the study reach (cf. colored polygons in Figure 6).

4.2. Classification and Surface Modeling

The methods for classifying the point cloud and modeling the surface of the water course described above were applied to the topo-bathymetric point clouds of the four epochs from April 2013 to October 2014. Figure 7a shows a color map of the water surface derived for the October 2014 dataset using the semi-automatic modeling approach described in Section 3.2. As indicated by the colors the ellipsoidal water level heights range from 261.0 m (dark blue) to 259.0 m (white).

Figure 7. Water surface modeling: Results (a); accuracy assessment (a+b); and refraction correction effect (c).

Figure 7. Water surface modeling: Results (a); accuracy assessment (a+b); and refraction correction effect (c).

A validation of the DWM was performed for a 50 m section of the meander based on simultaneous Terrestrial Laser Scan (TLS) which was co-registered with the simultaneous ALB flight using circular retro-reflecting targets. The 3D coordinates of the reference targets were measured with a total station (Leica TPS1200). The TLS (Riegl, VZ-1000) delivered 44908 measurements of the water surface and for each point the deviation from the modeled water surface was plotted (Figure 7a) and statistically analyzed (Figure 7b). The RMSE of the height discrepancies is 3.3 cm which results in a height error of the water echos of ca. 1 cm (assuming the standard refractive index of n = 1.33 for water). This measure is well below the ranging accuracy of the topo-bathymetric sensors (25 mm) according to the datasheets. It should be noted that the mean deviation of –1.7 cm is significantly different from zero (one-sample T-test with 3823 samples, confidence level: 95% , confidence interval: ±1.9606, test statistic: –36.8279). The possible error sources include (i) geodetic datum effects; (ii) water surface modeling errors; (iii) inclined water surface levels; and (iv) penetration of the green laser signal into the water column. From all these effects only the geodetic datum effect (i.e., misregistration between TLS and ALB scan) could be quantified based on reflections on dry land (gravel banks, median: 2.2 cm, sigma_MAD : 2.2 cm). The color-coded height discrepancies presented in Figure 7a show a systematic pattern. As the water surface modeled from green LiDAR echoes is rather rigid, the systematic effect can be attributed to the long-wave behavior of the DWM. However, it is noted that the reported values are smaller i.e., lower bias and std.dev.) than those presented in the related work (Section 1).

The validation results confirm that the reconstruction of the water surface based on sparse green laser echoes is feasible with an accuracy in the order of the ranging accuracy of the laser sensor. Figure 7c highlights the effect of the laser signal refraction in water by showing the raw and the corrected echoes.

Figure 8. Classified LiDAR points; (a) pre-classification of water echoes (water surface, volume backscatter, river bed); (b) DTM filtering (ground, non-ground).

Figure 8. Classified LiDAR points; (a) pre-classification of water echoes (water surface, volume backscatter, river bed); (b) DTM filtering (ground, non-ground).

As detailed in Section 3.4 ground point filtering is carried out in a two-step-process by first roughly separating river bed, water surface and volume backscatter (Figure 8a) and subsequently distinguishing ground from off-terrain points (Figure 8b) based on hierarchic robust interpolation [28]. As can be seen in the detail of Figure 8a not all the volume backscatter echoes were identified correctly in the pre-classification, but some were erroneously classified as river bed. However, most of the non-river bed points could be removed in the pre-classification step. This was verified visually for the investigated reach for all epochs. As our final goal is not a perfect classification of the water echoes, but a good DTM for subsequent habitat modeling, the occasional misclassifications are considered acceptable as the remaining off-terrain points within the river bed are correctly reclassified in the robust interpolation filtering procedure. Figure 8b shows the final classification results where all points are classified in terrain (river bed echoes and ground points of riparian area) and off-terrain (water surface, volume backscatter, vegetation). To preserve the sharp surface discontinuities in steep bank areas, where standard DTM filtering approaches tend to flatten the surface, additional 3D break lines were introduced in the DTM filtering.

4.3. DTM Quality and DoD Masking

In the introduction, the importance of river morphodynamics for the quality of the aquatic environment was highlighted. So far, repeated DTM surveys of the river bed had to be conducted by time consuming terrestrial survey, by optical remote sensing and/or, in mid-size to large rivers, by echosounding. For the investigated Pielach River the impacts of high flows (annual flood and a 30-years flood event) on river morphology and consequently habitat distribution have been evaluated and quantified based on topo-bathymetric LiDAR. Regardless of the employed data acquisition method the estimation of the uncertainty of the resulting DTMs is crucial for interpreting the detected changes. The method for deriving spatially variable DoD error estimates was detailed in Section 3.5 and the intermediate and final results of the procedure are shown in Figure 9 for the central meander bend. It is noted that the biases calculated from the terrestrial reference measurements are not considered in the error budget as they are eliminated in the DoD difference as long as they are constant in both epochs. This holds for all open areas (river bed, gravel banks, grassland) but the stage of the vegetation may well influence the bias within the alluvial forest. Thus, in the alluvial forest the DoD for different phenological states requires careful interpretation.

The DTM precision maps [69] of the February 2014 and October 2014 epoch are displayed in the first row of Figure 9. The green color tones indicating errors in the range of 0–3 cm are predominant, especially in the February 2014 dataset (leaf-off), whereas more yellow to red color tones representing errors of 5–10 cm appear in the October 2014 dataset. The higher error level is mainly caused by the lower ground point density due to incomplete defoliation in mid-October. The sub-optimal vegetation penetration resulted in data voids and, as a consequence, in unusable (white) DTM areas. The second row shows the DTM accuracy maps where the surface roughness derived from the ground based reference measurements was added based on a GIS-based delineation of homogeneous land cover areas and the unmasked DoD for the February–October 2014 period. Negative differences (erosion) are hereby visualized in red and positive differences (deposition) in blue color tones. The last row in Figure 9, finally, shows maps of the propagated DoD accuracy, the student’s t-score, and the masked DoD thresholded at the 95% confidence interval. These changes are highly significant. The DoD accuracy map documents higher accuracies within the channel area compared to the adjacent riparian forest. This is remarkable and mainly results from the high LiDAR bathymetry point density of ca. 20 points per square meter which drops only for the pool areas (depth: ca. 3 m) at the maximum depth performance of the sensor.

Figure 9. DoD accuracy estimation; first row: DTM precision maps and color bars; second row: DTM accuracy maps and unmasked DoD; last row: combined DoD accuracy, student’s t-score map, masked DoD thresholded at 95% confidence interval.

Figure 9. DoD accuracy estimation; first row: DTM precision maps and color bars; second row: DTM accuracy maps and unmasked DoD; last row: combined DoD accuracy, student’s t-score map, masked DoD thresholded at 95% confidence interval.

Figure 10 and Figure 11, and Table 4 show the DoD and the erosion and deposition, respectively, both masked and unmasked by the 95% confidence test. The unmasked DoD in Figure 10a show light green color tones indicating moderate deposition outside the main channel. However, these areal positive differences do not denote topographic changes but can be attributed to the slightly higher DTM surface in October 2013 compared to April 2013 due to incomplete defoliation. Especially the understory is still in partial leaf and prevents the laser signal from fully penetrating the vegetation layer. In Figure 10b showing the DoD of the October 2013–February 2014 period the same areas appear in light red color tones indicating artificial erosion resulting from the full vegetation penetration in February 2014. The mean height deviations within the alluvial forest are in the range of 5 cm (April 2013/October 2013: –7.4 cm, October 2013/February 2014: 4.9 cm) which is in good conformance with [86]. Figure 10 deliberately shows the raw DoD not considering the mask indicating significant changes based on the DoD error estimation. In Table 4 the raw and the masked volumetric changes are presented. The figures underline that most of the changes within the alluvial area (HQ₁ are statistically insignificant (e.g., April 2013/October 2013; raw deposition: 10.631 m³ textitvs. masked deposition: 877 m³). The same conclusions can also be drawn by comparing the raw volumetric changes of the April–October 2013 and the October 2013–February 2014 period. Whereas the former is dominated by (artificial) deposition (10.631 m³), the latter is predominantly erosional (–8.964 m³). Especially the high erosion rate of the latter is implausible, as no flooding of the riparian area took place in that period. The respective masked DoD values are much more balanced. From these findings we conclude that (i) the derivation of significant changes based on spatially variable DoD accuracy measures is crucial for the interpretation of geomorphological dynamics and (ii) reliable estimates of volumetric changes in the alluvial area can only be derived based on leaf-off datasets. Figure 10c shows the changes between the leaf-off datasets April 2013 and February 2014. Within the floodplain the differences are predominantly close to zero indicated by the light gray color tones as a result of the full vegetation penetration in both datasets.

Within the channel area (MQ), however, the majority of the reported raw changes are significant. As seen in Table 4 (columns: MQ raw, MQ masked) the values deviate less from each other compared to the values of the overbank area (HQ₁). This holds true especially for the April 2013 to October 2014 period (raw deposition: 6.906 m³ vs. masked deposition: 5.470 m³).

Table 4. Volumetric changes derived from DoD. Separate values are reported for: deposition and erosion, mean flow (MQ) and annual flood (HQ₁) areas, and raw and statistically significant (masked) changes.

Volumetric Change [m3]		Raw		Masked (95% Confidence Interval)
Period	Change	MQ	HQ1	MQ	HQ1	Total	Sum
Apr13–Oct13	deposition	4.749	10.631	2.565	877	3.442	435
	erosion	–3.353	–2.134	–2.089	–918	–3.007
Oct13–Feb14	deposition	1.548	1.307	554	175	729
	erosion	–2.981	–8.964	–1.056	–174	–1.230	–501
Apr13–Feb14	deposition	4.040	3.664	2.282	642	2.924
	erosion	–4.078	–2.823	–2.562	–885	–3.447	–523
Feb14–Oct14	deposition	5.793	10.634	4.161	1.006	5.167
	erosion	–5.021	–8.093	–3.495	–5.754	–9.249	–4.082
Apr13–Oct14	deposition	6.906	11.690	5.470	1.610	7.080
	erosion	–6.171	–8.308	–4.527	–5.958	–10.485	–3.405

Table 4. Volumetric changes derived from DoD. Separate values are reported for: deposition and erosion, mean flow (MQ) and annual flood (HQ₁) areas, and raw and statistically significant (masked) changes.

Volumetric Change [m3]		Raw		Masked (95% Confidence Interval)
Period	Change	MQ	HQ1	MQ	HQ1	Total	Sum
Apr13–Oct13	deposition	4.749	10.631	2.565	877	3.442	435
	erosion	–3.353	–2.134	–2.089	–918	–3.007
Oct13–Feb14	deposition	1.548	1.307	554	175	729
	erosion	–2.981	–8.964	–1.056	–174	–1.230	–501
Apr13–Feb14	deposition	4.040	3.664	2.282	642	2.924
	erosion	–4.078	–2.823	–2.562	–885	–3.447	–523
Feb14–Oct14	deposition	5.793	10.634	4.161	1.006	5.167
	erosion	–5.021	–8.093	–3.495	–5.754	–9.249	–4.082
Apr13–Oct14	deposition	6.906	11.690	5.470	1.610	7.080
	erosion	–6.171	–8.308	–4.527	–5.958	–10.485	–3.405

4.4. River Morphodynamics

The morphodynamic changes of the fluvial topography for the period from April 2013 to February 2014 are documented in Figure 10. It can clearly be seen from Figure 10a that the annual flood event in June 2013 (Q_{2013−06−25} = 140.0 m³·s⁻¹, HQ₁ = 136.6 m³ ·s⁻¹) resulted in a relocation of gravel bars (cf. detail A) as well as in substantial bank erosion especially at the steep bank within the meander bend (B). The dark blue and red colors indicate a vertical thickness of the eroded and deposited sediment of more than 0.50 m (mean/max thickness, A (deposition): 0.32/0.81 m, A (erosion): –0.33/–1.51 m) and the edge of the steep bank receded (B) by more than 1 m. The impact of the 2013 flood was mainly restricted to the main channel (cf. black MQ polygon marked in Figure 10).

Figure 10. Color-coded visualizations of the raw (unmasked) DoD documenting the changes caused by the annual flood event in June 2013; (a) Period April 2013–0ctober 2013; (b) Period October 2013–February 2014; (c) Period April 2013–February 2014.

Figure 10. Color-coded visualizations of the raw (unmasked) DoD documenting the changes caused by the annual flood event in June 2013; (a) Period April 2013–0ctober 2013; (b) Period October 2013–February 2014; (c) Period April 2013–February 2014.

Figure 11. Color-coded visualizations of the masked DoD in front of the hillshade of the more recent epoch. (a) Impact of 30-years flood in May 2014; (b) Total impact of the entire studied period, Photo documentation of morphodynamics; (c) bank erosion (i.e., detail B in Figure 10a); (d) gravel bank relocation (detail E in Figure 11a); (e) new meander shortcut channel (detail G); (f) deadwood pile (detail H).

Figure 11. Color-coded visualizations of the masked DoD in front of the hillshade of the more recent epoch. (a) Impact of 30-years flood in May 2014; (b) Total impact of the entire studied period, Photo documentation of morphodynamics; (c) bank erosion (i.e., detail B in Figure 10a); (d) gravel bank relocation (detail E in Figure 11a); (e) new meander shortcut channel (detail G); (f) deadwood pile (detail H).

Furthermore, the dynamics within the main channel from October 2013 to February 2014 should be noted. In that period of time only a single noteworthy discharge , but still less than HQ₁ (Q_{2014−11−24} = 80.0 m³·s⁻¹) occurred. Whereas occasional negative differences in pool areas (cf. detail C) can be attributed to the better depth performance of the February 2014 flight (low discharge, very clear water), detail D documents sediment transport triggered by discharges below HQ₁. The discharge threshold to start incipient sediment motion is still discussed [87]. The captured bathymetry data indicate that sub annual flood discharges are sufficient to trigger bed load transport within the study reach.

In contrast to the June 2013 event the extreme flood on 16 May 2014 (Q = 369 m³·s⁻¹, HQ₃₀ = 347 m³·s⁻¹) did not only effect the channel but the entire riparian area. Figure 11 shows the significant changes over the hillshade of the prior of the two compared epochs. Besides massive gravel relocation (E) and bank erosion (F) both erosion and deposition also occurred in the alluvial forest. Detail G shows a new shortcut channel bypassing the meander bend. The deposition areas (H) could belinked to accumulations of woody debris. Entire trunks, branches and twigs were transported from upstream by the flood wave and were caught between the large trees within the alluvial forest. In the course of a local inspection in November 2014, the fluvial changes caused by the 30-years flood event could be verified in-situ and some examples have been mapped with a low cost GNSS navigation device (Garmin Oregon 550t). The captured tracks are marked in Figure 11b summarizing the total changes caused by fluvial activity within the entire captured time range (I: newly formed shortcut channel; J+K: new gravel banks; L+M: areas with piles of deadwood; note the overlap of I and L with a massive pile of deadwood covering the engraved channel). Furthermore, the photographs in Figure 11c–f document the aforementioned features.

The last column in Table 4, finally, draws the conclusion that the 30-years event in 2014 caused a noteworthy erosion of the entire study reach (–4.082 m³), whereas the annual flood event in 2013 only exhibited a moderate erosion rate (–523 m³). Scouring is more pronounced in the riparian area (cf. new meander shortcut) whereas the changes in the main channel are in the state of equilibrium.

4.5. Flood Studies

Floods are the driving agents for morphological changes and thus the source for creation and renewal of instream habitats (e.g., [88]). In the present study DTMs based on topo-bathymetric LiDAR have been tested for flood modeling and the creation of flow-velocity and water depth maps as the basic input for habitat modeling (Figure 12).

The results show a complex distribution of the hydraulics (e.g., depth-averaged flow velocities) in the main channel caused by the mesounit scale variability of the river bed morphology. These high-resolution and high-accuracy data have to be seen as an important improvement for further numerical analysis and developments (e.g., detection of the initiation of sediment motion, numerical modeling of bank erosion) as well as a link to instream habitat studies (e.g., flood pulse investigations according to [89]) as integrative approaches are required by the aims of the European Water Framework Directive [5].

Figure 12. Results of two-dimensional depth-averaged hydrodynamic-numerical modeling at the Pielach River based on the annual flood discharge as hydrological input (136 m³·s⁻¹) for (a) water depth; (b) depth-averaged flow velocities; (c) detail including flow vectors).

Figure 12. Results of two-dimensional depth-averaged hydrodynamic-numerical modeling at the Pielach River based on the annual flood discharge as hydrological input (136 m³·s⁻¹) for (a) water depth; (b) depth-averaged flow velocities; (c) detail including flow vectors).

4.6. Habitats

To evaluate the impact of morphodynamic changes on instream habitats, integrative modeling approaches are required to predict the effects of e.g., habitat alterations on population of aquatic organisms (macroinvertebrates, fish). Here, so called habitat selection models (often referred to as “preference or habitat index” models or simply habitat models) are widely used [90]. These models are based on observing the frequency with which animals use various habitat types and the availability of the habitat types. The ratio of habitat use to habitat availability is then transformed into a measure of habitat selection [91]. Habitat selection can be demonstrated if e.g., fish are found in higher densities in particular suitable habitats (habitat suitability), or if fish are found at higher frequencies in particular units of the river relative to the frequencies of those units in the aquatic environment (habitat preference) [92]. For instream flow studies the selection and application of habitat investigations on various scales has to be seen a crucial component for the integrative analysis between possible changes of the physical environment (e.g., changes in flow) and the response of the biota to those changes [93,94,95,96,97]. Instream habitats, which are strongly hierarchical, can be analyzed (e.g., numerically modeled) at a variety of spatial scales [93,98], for which several scaling concepts are available (e.g., [93,95,99]).

For testing bathymetric LiDAR both the microunit (10⁰ m) and the mesounit (10¹ m) scale have been selected [93]. For characterization of the physical environment (e.g., flow velocity, water depth) a depth-averaged hydrodynamic-numerical model has been applied using DTMs derived from topo-bathymetric LiDAR for modeling. At the Pielach River the microunit scale has been investigated based on habitat suitability data of the target fish species nase (Chondrostoma nasus). Suitability for spawning and for juvenile fish have been obtained by scuba-diving for 1250 observation points of fish under low flow conditions [57]. The highest suitability of spawning has been detected in areas with shallow water depths (< 0.3 m) and high flow velocity (> 0.8 m·s⁻¹) [57]. Moreover, habitat suitabilities of juvenile nase exhibited the highest distribution in shallow water depths (<0.3 m) with low flow velocities (< 0.1 m·s⁻¹). Applying this information of habitat requirements to the numerically modeled depth-averaged flow velocities and water depths at the Pielach River, (micro) habitat maps could be derived based on the PHABSIM approach (multiplication of suitability indices, [100], presented in Figure 13b,c. The results for mean flow conditions (common during spawning and the juvenile stage) showed on the one hand suitable habitats for both life stages and on the other hand a clear separation between those two tested biotic requirements (Figure 13b,c). The differences have to be related to morphological heterogeneity and, thus, variability in flow velocity patterns and water depth. Morphological heterogeneity is also one of the main drivers for the second applied habitat modeling approach, the mesohabitat evaluation based on the MEM-concept [9]. Here, discrete hydraulic patterns are grouped concerning similarities in flow velocity, water depth and bottom shear stress (substrate). The MEM-approach allows a differentiation of six discrete mesohabitats ranging from the high energetic riffle to low energetic backwater habitats (Figure 13a). Especially the mesohabitat results underline the quality of the (nearly) undisturbed river environment of the Pielach River (e.g., self-forming processes due to bank erosion) as all habitat types are present for the investigated study site (under mean flow conditions).

Comparison of the mesohabitat distribution for the entire monitoring period (April 2013–October 2014) showed that the morphydynamic processes caused only minor changes in the habitat composition. Although the mean flow simulations showed minor fluctuations in the habitat distribution for runs, which exhibited a decrease due to morphodynamic changes for the one-years flood and an increase due to the 30-years event, the distribution of most of the mesohabitats remained in a dynamic equilibrium which is a quality indicator for the hydro-morphological status of the investigated reach.

Figure 13. Application of LiDAR bathymetry for instream habitat studies on the mesohabitat scale (a) as well as on the microunit scale documenting the habitat suitabilities of (b) spawning of nase (Chondrostoma nasus); and (c) juvenile nase; all Figures are derived for mean flow conditions = 7.16 m³·s⁻¹; Suitability Index (SI)data based on [57]).

Figure 13. Application of LiDAR bathymetry for instream habitat studies on the mesohabitat scale (a) as well as on the microunit scale documenting the habitat suitabilities of (b) spawning of nase (Chondrostoma nasus); and (c) juvenile nase; all Figures are derived for mean flow conditions = 7.16 m³·s⁻¹; Suitability Index (SI)data based on [57]).

Habitat modeling based on the topo-bathymetric LiDAR point cloud indicated that all six types of hydro-morphological units were present for the studied discharge range in all captured epochs (cf. Figure 13). Important low-energetic habitats like shallow water and backwater sites are available for small fish for so called flow-pulse events within the bankfull width. The same preferred habitat types for small fish, however, may function as refugee habitats in case of high flow events (flood-pulse). Here, the need for habitat shifting at the Pielach River occurs only within the spatial extent of hydro-mophological units with the additional opportunity to use downstream drift to reach these suitable or refugee habitat conditions. This has to be seen different in regulated rivers, or in rivers which have been only restored within a regulated profile. In such semi-natural or regulated rivers there are huge declines in available and suitable habitats if the wetted width has already reached the banks. The increase in fast run habitats (flow velocity > 1 ms⁻¹) with increasing flow is dominantly reducing shallow water- and backwater-areas to almost zero due to the regulated banks [101] or is causing excessive downstream drift of juvenile fish. Thus, a stable overall habitat distribution at the reach scale is an indicator for the quality of the investigated study area as a natural self-forming reach of the Pielach River. This especially holds true as the dynamic component of hydrology and river morphodynamics (erosion and deposition due to floods) are considered on a small scale. Although, the impact of the annual flood event and the thirty-years flood exhibited significant changes in morphology (downstream migration of gravel bars), however, no significant overall quantitative habitat changes on the mesounit scale occurred (cf. Figure 14).

Figure 14. Variability of area distribution of hydro-morphological units in the period April 2013–October 2014 analyzed for mean flow (7.16 m³·s⁻¹, lower part of the bars) and annual flood (136.6 m³·s⁻¹, entire bars) conditions

Figure 14. Variability of area distribution of hydro-morphological units in the period April 2013–October 2014 analyzed for mean flow (7.16 m³·s⁻¹, lower part of the bars) and annual flood (136.6 m³·s⁻¹, entire bars) conditions

5. Conclusions

In this article it was shown that monitoring of alpine rivers (riffle-pool type gravel rivers, mean discharge <10 m³·s⁻¹, stream order ≤ 4) by topo-bathymetric LiDAR is feasible. A 1500 m long reach of the Pielach River in Lower Austria was repeatedly surveyed with topo-bathymetric LiDAR.

The suggested method for modeling the water surface was compared to terrestrial reference measurements and showed a bias below 2 cm and a standard deviation of 3 cm. In comparison to previous studies evaluating the offset between air/water interface and surface reflections our modeling approach shows smaller deviations.

To overcome the problems of classification into land surface and vegetation points at steep banks we introduced an additional explicit terrain break line modeling step. For the wetted perimeter a new classification method based on measured and derived point features allowed separation of the laser echoes into river bottom, volume backscatter, and water surface. Together with the high and homogeneous point density (>20 points per square meter), this lead to a highly detailed description of the river bed topography. This high resolution enabled a detailed analysis of the river morphodynamics, both for the main channel and the floodplain. Terrain changes could be related to bank erosion and sediment transport induced by a 30-years flood (May 2014), an annual flood (June 2013), and even smaller events (December 2013).

In this study the t-test was applied to identify significant changes in each cell of the DoD independently. The extension of Wheaton et al. [46] to consider spatial coherence of height changes was not applied because of the (small) systematic errors found in the data. Especially in the riparian forest the terrain elevation is too high during leaf-on data acquisition. This systematic error is not necessarily constant but can exhibit a spatial pattern, e.g., due to small clearances. It will require further development of sensors, terrain reconstruction algorithms, or testing procedures to reliably determine differences between data acquisition performed at different phonological states. For now, the interpretation of elevation changes in the riparian forest should be performed for leaf-off acquisitions only.

Modeling of DTM errors has shown that, apart from sealed roads, the channel ground surface is determined most accurately followed by the grassland in the floodplain. The elevation of the terrain in the riparian forest is determined with the least accuracy. The reason is that the channel features a clear surface definition (gravel bed) and a high point density of the LiDAR bathymetry, whereas the surface definition for the grassland exhibits higher uncertainty. For the riparian forest additionally the point density drops as the ground is only partially visible through the vegetation. Thus, in contrast to the overbank areas the river channel accuracy depends less on the flight season. Monitoring of fluvial changes within the channel area based on topo-bathymetric LiDAR is therefore feasible throughout the entire year.

Mesohabitat modeling was performed using the MEM-approach. Based on flow velocities, water depths and bottom shear stress, a distinction into riffle, fast run, run, pool, backwater, and shallow water habitats was possible. It could be demonstrated that all six mesohabitat classes are available in each of the analyzed datasets indicating the resilience of the study reach against disturbances. Repeated survey and habitat modeling showed that the overall distribution of the hydro-morphological units did not change largely, although the individual elements notably shifted. Based on the results of the mesohabitat modeling, suitability indices for the target fish species nase were derived (microunit scale) and underline the suitability of the studied reach for both spawning and juvenile life stages of nase.

We finally conclude that topo-bathymetric LiDAR at high spatial and temporal resolution is a viable method for monitoring of instream habitats but also for studying potential floods, thus contributing to various European directives (Fauna-Flora-Habitat, Water Framework, Flood).