Reading the bed morphology of a mountain stream : a geomorphometric study on high-resolution topographic data

Abstract. High-resolution topographic data expand the potential of quantitative analysis of the earth surface, improving the interpretation of geomorphic processes. In particular, the morphologies of the channel beds of mountain streams, which are characterised by strong spatial variability, can be analysed much more effectively with this type of data. In this study, we analysed the aerial LiDAR topographic data of a headwater stream, the Rio Cordon (watershed area: 5 km 2 ), located in the Dolomites (north-eastern Italy). The morphology of the channel bed of Rio Cordon is characterised by alternating step pools, cascades, and rapids with steps. We analysed the streambed morphology by means of ad hoc developed morphometric indices, capable of highlighting morphological features at a high level of spatial resolution. To perform the analysis and the data interpolation, we carried out a channel-oriented coordinate transformation. In the new coordinate system, the calculation of morphometric indices in directions along and transverse to the flow direction is straightforward. Three geomorphometric indices were developed and applied as follows: a slope index computed on the whole width of the channel bed, directional variograms computed along the flow direction and perpendicular to it, and local anomalies, calculated as the difference between directional variograms at different spatial scales. Directional variograms in the flow direction and local anomalies have proven to be effective at recognising morphologic units, such as steps, pools and clusters of large boulders. At the spatial scale of channel reaches, these indices have demonstrated a satisfactory capability to outline patterns associated with boulder cascades and rapids with steps, whereas they did not clearly differentiate between morphologies with less marked morphological differences, such as step pools and cascades.


Introduction
Channel bed morphology, resulting from the interactions between hydrology, topography and sediment supply, is sensitive to climate changes, land use changes, and to disturbances (natural or anthropogenic) in the upstream river network.Accordingly, the study of channel morphology provides information on the influencing processes mentioned above.
Small high-gradient streams in mountainous areas have a complex and dynamic morphology, strongly influenced by the topographic ruggedness of the basin and the spatial variability of sediment supply.Geomorphic processes, such as debris flows and landslides are widespread and frequent in mountainous basins, and they often interact with the channels.These extrinsic episodic processes have a quite immediate effect on local bed morphology.
The characterisation of topography and morphology of mountainous channels, the assessment of changes caused by large floods and the objective comparison between different channels require reliable data and procedures.However, quantitative studies of mountain streambed morphology are challenging both with regard to (1) the acquisition of extensive high-resolution topographic data and (2) difficulties in classifying complex stream bed morphology.
1.The topographic survey of mountainous streams requires collecting large amounts of data (Grant et al., 1990;Wohl et al., 1997;Chin, 1999;Lenzi, 2001;Milzow et al., 2006).Surveys carried out using traditional techniques, such as theodolite and stadia, may lead to high-quality representations of channel-bed topography, but are expensive and time-consuming.This can be overcome by means of Airborne LiDAR technology, which makes it possible to quickly collect large amounts of high-resolution topographic data, allowing careful surveys of large areas at reasonable costs and enabling the implementation of high-resolution DTMs.
S. Trevisani et al.: Reading the bed morphology of a mountain stream 2. The morphology of channel beds is characterised by a hierarchical structure with features that can be recognised at both the spatial scale of morphological units (e.g., steps and pools) and of stream reaches (e.g., step pool sequences, cascade, and plane bed).The classification of stream reaches into unequivocally defined classes is difficult when the channel bed morphology is complex, with short channel reaches of different morphologies and fuzzy transitions between them.
A geomorphometric approach, i.e., the quantitative study of morphology (Pike, 2000;Hengl and Reuter, 2009), can provide valuable elements for interpreting the complex morphologies of mountainous streambeds and for handling the subjectivity intrinsic in the classification of bed morphologies by means of field observations.Airborne LiDAR data have been widely applied in many geomorphological fields, especially for studying landslides (Glenn et al., 2006;Ardizzone et al., 2007) and sediment dynamics (Scheidl et al., 2008;Brown et al., 2009).In studies on riverine environments, early applications of LiDAR data were mainly intended to provide improved DTMs for applying numerical flood models.Recent applications of aerial LiDAR data to the analysis of fluvial environments include the geomorphological mapping of gravel-bed rivers (Charlton et al., 2003), tidal channels (Mason et al., 2006) and river valley environments (Jones et al., 2007), the evaluation of river bank erosion (Thoma et al., 2005) and long-term changes caused by fluvial and debris-flow activity (Magirl et al., 2005), the assessment of the slopes in headwater channels (Vianello et al., 2009), the hillslope to valley transition (Tarolli and Dalla Fontana, 2009) and the delineation of the water surface (Höfle et al., 2009).
In this study, morphometric indices derived from terrain analysis (Wilson and Gallant, 2000) and geostatistics (Goovaerts, 1997) are used to characterise channel bed morphologies through the analysis of a high-resolution DTM (cell size 0.5 m) developed from aerial LiDAR data.We read the streambed morphology using ad-hoc developed morphometric indices capable of highlighting morphological features at a high level of spatial resolution.Findings arising from the analysis of the morphometric indices are checked against field evidences to evaluate the correspondence with the morphological units classified by means of field observations.

Study area
The study area is the main channel of the Rio Cordon, a small stream in the Dolomites (Eastern Italian Alps).The drainage basin of the Rio Cordon (Fig. 1) covers 5 km 2 ; the elevation ranges between 1763 and 2748 m a.s.l., with an average value of 2200 m a.s.l.; the average slope is 27  setting of the basin is rather complex: dolomite crops out in the upper part, whereas in the central and lower parts volcanic conglomerates, sandstones and calcareous-marly rocks crop out; moraines, scree deposits and landslide accumulations are also widespread.
The basin presents a mean annual rainfall of approximately 1100 mm; precipitation occurs mainly as snowfall from October to April.Runoff is dominated by snowmelt in May and June; floods mostly occur in the summer and early autumn.Vegetation cover consists mainly of mountain grassland (60%) and widespread shrubs (15%), while forest stands composed by spruce and larch are found only in the lower part of the watershed and occupy 6% of the total area.Unvegetated areas (bare rock and scree) are common in the upper parts of the watershed.No relevant artificial structures, such as check dams or channel lining are present along the channels.
In this study, we consider the lower part of the main channel reach, starting from a rocky gorge in the central part of the basin, and ending at the basin outlet (Fig. 1).The analysed channel reach has a total length of 1625 m and ranges in elevation from 2098 to 1818 m (ellipsoid height) (Fig. 2).The drainage area varies from 2.23 km 2 at the upstream end of the studied channel reach to 5 km 2 at its downstream end.Lenzi (2001) presented a classification of the channel bed morphology of the Rio Cordon and reported the geometric parameters of step-pool structures, measured using a total station.The classification of Lenzi (2001), also followed in Cavalli et al. (2008), identified three morphologies in the studied channel reach: a cascade reach at the outlet of the rocky gorge, two riffle-pool sequences, separated by a short "mixed" reach, and an alternation of step-pool and "mixed" reaches in the lower part of the stream.Although no relevant changes in the morphology have occurred in the Rio Cordon after the surveys reported by Lenzi (2001), the analysis of bed morphology was revised to achieve a more detailed classification of channel reaches.The revised classification envisages three main morphologies defined following the channel-reach morphology classifications proposed by Grant et al. (1990) and Montgomery and Buffington (1997): cascade, step pool and rapid with steps.Cascade reaches are characterized by longitudinally and laterally disorganized bed material typically consisting of cobbles and boulders, whereas step-pool reaches are characterized by longitudinal steps formed by large clasts organized into discrete channelspanning accumulations that separate pools containing finer material (Montgomery and Buffington, 1997).The cascade reach at the outlet of the rocky gorge, which is characterised by abundant large boulders, was classified as a boulder cascade to better emphasize the coarse grain size of this channel reach.One of the two channel reaches previously classified as riffle-pool was interpreted as a rapid with steps (Grant et al., 1990), characterised by an organisation of small and medium boulders into irregular ribs oriented perpendicular to the channel and exposed to low flow (Fig. 3), whereas the second displays sufficiently defined steps and was ascribed to two neighbouring step-pool and stepped cascade reaches.The term "mixed reach," which identifies channel reaches similar to step-pool sequences but with coarse particle bars deposited upstream of the boulder steps or adjacent to isolated big boulders (Billi et al., 1998), was abandoned, and the channel reaches previously classified as "mixed" were classified as cascades.The presence of large woody debris in the studied channel is very limited and does not influence channel morphology.The upper and lower ends of channel reaches were checked using GPS, leading to slight differences from the previous sequences of step-pool and mixed reaches.Figure 2 shows a longitudinal profile of the studied channel, with the classification of channel bed morphologies.

LiDAR data
LiDAR and photographic data were acquired from a helicopter using an ALTM 3100 OPTECH and a Rollei H20 digital camera.The survey was carried out in snow-free conditions in October 2006, flying at an average altitude of 1000 m above ground level.The flying speed was 80 knots, the scan angle 20 degrees and the pulse rate 71 KHz.The survey design point density was specified to be greater than 5 points/m 2 , recording up to four returns, including the first and last.The shallow depth of water in the analysed stream at the time of the LiDAR flight (not exceeding a few tens of centimetres) made it possible to survey the channel bed without using a bathymetric LiDAR sensor.However, some attenuation of the LiDAR signal has been observed in the pools of step-pool reaches.
LiDAR point measurements were filtered into returns from vegetation and bare ground.A bare ground dataset was generated by eliminating non-ground points by manual editing, through visual inspection of aerial photos and using 2-D plots (elevation versus distance along the channel, Fig. 4), which allows points located in vegetated areas to be identified.This approach, which is fairly time-consuming but still affordable thanks to the small extent of the analysed area (8144 m 2 ), was chosen because it permits a careful control on the filtering process.Manual filtering removes   misclassification errors (i.e., points falling on trees classified as ground), which would have appeared as strong local outliers from a geostatistical perspective.The filtering of less obvious errors, like those related to instrumentation, is handled during the interpolation sequence.The average density of ground points was 5.46 points/m 2 , with a standard deviation of 3.69.This value indicates that the manual filtering approach can preserve the high density of the raw dataset.Mean and median separations between ground points, determined with the nearest neighbour technique, are 0.19 and 0.17 m, respectively.

Domain delineation and coordinate transformation
This study required the definition of the spatial analysis domain, which corresponds to the channel bed, by delineating the lower edges of the banks.Elevation data falling on the channel banks, which could bias the interpolation and computation of morphological indices, were excluded from the analysis.
In order to digitise and georeference in a GIS the outline of the channel bed, different thematic layers were analysed interactively with the support of aerial orthophotos with a horizontal resolution of 0.15 m.The edges corresponding to the feet of the banks were detected using maps of shaded relief, slope, roughness (Cavalli et al., 2008) and openness index (Yokoyama et al., 2002).These maps were derived from a preliminary DTM with pixels of 0.5 m×0.5 m.The derived channel bed outline ranges in width from 1.5 to 10.5 m, with a mean value of 5 m (first and third quartiles of 4 m and 6 m).
Subsequent processing of LiDAR data was performed on a transformed coordinate system (Fig. 5).Spatial analysis in transformed coordinate systems has found many applications in the field of geostatistical geomodelling, in environmental geostatistics, and in river geomorphology (Deutsch, 2002;Merwade et al., 2005;Legletier and Kyriakidis, 2007a).The coordinate transformation was performed using spatial tools available in most GIS packages, such as spatial queries, spatial joins and shape discretisations.In the present study, we used ArcGIS 9.3 (ESRI, 2008).The procedure is composed of the following steps: -Derivation of the centreline: this is computed as the line equidistant from the two lines defining the boundary of the domain.This line must be smoothed to avoid having a radius of curvature less than the half-width of the channel; otherwise it would not be possible to perform the transformation (Goff et al., 2004;Legletier and Kyriakidis, 2007b).
-Centreline discretisation: the smoothed centreline is discretised as points with a close spacing (step of 1 cm), ensuring a high-accuracy transformation of the coordinate system.
-Computation of stream-centred coordinates: in the new system, the position of a given point in the river bed is defined by an abscissa, RX, and an ordinate, RY.RX indicates the distance (positive or negative depending on the side) from the nearest point discretising the centreline.RY corresponds to the position of the nearest point of discretisation on the centreline (Fig. 5).
The derivation of the centreline is quite straightforward in the Rio Cordon, which has a narrow channel with low sinuosity.In the case of high sinuosity and wider streams (e.g., meandering rivers), the centreline derivation could require more trials in the smoothing process, to get a curvature radius less than half of the channel width.
The transformation of the coordinate system from a Cartesian to a stream-centred one has advantages both for the interpolation (i.e., DTM generation) and for spatial analysis.DTM interpolation on the stream-centred coordinate system allows the various parameters and models (variogram anisotropies, search ellipse, trend calculation) to be adjusted in a system in which the main directions, along and transverse to the stream flow direction, have a sound physical basis.In the transformed coordinate system, the calculation of directional morphological indices in directions orthogonal and parallel to river flow becomes straightforward, and indices evaluating the asymmetry of bedform organisation with respect to the centreline can be easily calculated.Moreover, the stream-centred coordinate system emulates the viewpoint of the expert in the field, who surveys a river channel taking as the main directions the flow direction of water (channel profile) and the direction orthogonal to it (cross-sections).

DTM derivation
The DTM was derived via ordinary block kriging interpolation, with a pixel size of 0. and the need to describe local (i.e., large-scale) morphological features.The high data density (5.46 points/m 2 ) allowed us to use a simplified interpolation approach, avoiding the use of more complex procedures, such as universal kriging or local kriging (Stroet et al., 2005).
Block kriging interpolation was carried out using an elliptical search window with a semiaxis of 1 m in the flow direction and of 2 m across it.The short semiaxis in the direction of flow was dictated by the need to attenuate the effect of trend, i.e., the downstream decrease of elevation.For the distances of interest, the variogram was modelled by a spherical variogram with a range of 4 m and a sill of 0.3 m 2 plus a nugget of 0.01 m 2 .We considered the nugget as representative of measurement error and not of microvariability (Cressie, 1993), corresponding to white noise with a standard deviation 0.1 m.The imposition of the nugget as error permits us to filter out the white noise (i.e., not spatially correlated) component of error measurement.Then, the interpolation via block kriging, filtering out the intra-block spatial variability, is capable of removing short-wavelength spatially-correlated error.
The quality of the derived DTM was evaluated quantitatively by means of standard deviation maps, cross-validation results and maps, and qualitatively during the field surveys visually comparing the shaded relief map with the real streambed morphology.

Morphometric indices
Morphometric indices have been computed to characterise the morphology of the channel bed of the Rio Cordon.Except for the local channel slope, which is calculated based on elevation data, the other indices were calculated based on the residuals of elevation (Wooldridge and Hickin, 2002;Zimmermann et al., 2008;Trevisani et al., 2009), which were derived by removing the large-scale spatial variation from the interpolated DTM, i.e., removing the trend.The surface of residuals represents the high-frequency spatial variability component of the streambed morphology.The trend surface was derived by means of a local polynomial approach, fitting a planar model in a moving window covering the whole width of the stream for a length of 10 m in the direction of flow.
The morphological indices were calculated within moving windows; the calculated value was attributed to the centre of the window (Fig. 6).The indices were computed every 1 m along the stream centreline, using moving windows covering the whole width of the channel.
Both the DTM resolution and the moving window sizes were chosen taking into account the sizes of the topographic features being studied.The DTM cell size of 0.5 m is suitable to correctly represent the step pools' grain size, which is characterised in the Rio Cordon by a D 84 of about 1 m (Lenzi, 2001).Two different moving windows lengths (2 and 10 m) along the flow direction were used.The former (SW, small window) was chosen to detect single morphological units or small-scale features like steps, small pools and boulders; the latter (LW, large window) can describe the average spatial characteristics of a longer stream reach covering more morphological units (the mean step-pool spacing reported by Lenzi (2001) ranges from about 2.5 to 5 m), thus facilitating the recognition of morphological sequences.
The computation of the indices at regular intervals along the longitudinal channel profile, for windows covering the whole channel width, makes possible to revert the 2-D analysis to a 1-D issue, which can be tackled using time series analysis (Malamud and Turcotte, 1999); the distance along the stream centreline corresponding to time.An important advantage arises in the interpretation of the indices, which can be plotted as a profile along the channel, making it straightforward to visually compare with the topographic channel profile and morphological classification.A number of trials on different indices and computation parameters have been carried out to select the indicators to analyse the morphology of the Rio Cordon.The tested indices span from terrain analysis indices, such as directional slope and directional curvature, to spatial statistical indices, such as directional variograms calculated at different lags and parameters of variogram models fitted on experimental local variograms (Butler et al., 2001).The selection of geomorphometric indices has been based on a balance of complexity and informative content.
One basic morphometric index is the longitudinal slope, calculated by fitting a line, via least squares, to the elevation versus RY-coordinate plot.The regression was carried out on the subset of points falling within each moving window.This index, when calculated within SW (SW slope), can detect local features of channel bed morphology.When calculated on LW (LW slope), it gives overall information on the hydraulic energy of the investigated reaches.It is important to stress that the slope index is not calculated along the centreline of the channel, but considers the whole channel width, and is thus able to synthesise the average slope of the 2-D channel on a 1-D representation.The effectiveness of a 1-D representation considering the whole channel width for morphological analysis on step pool channels was demonstrated in a flume study by Zimmermann et al. (2008).The authors used average long profiles, constructed by averaging the bed elevation across the channel for different cross-sections, to objectively classify step pool units.A step toward a finer characterisation of channel-bed morphology is the calculation of directional variograms γ (h) (1) for the shortest lag (0.5 m) (Fig. 6).
(1) with γ (h) = γ (−h) As outlined in Eq. ( 1), the variogram gives, for each value of the vector h, an estimate of the half-mean squared differences between the N couples of points z(u α ) separated by the vector h.Thus, the variogram, for a specific lag, gives a measure of spatial variability functional to separation distance, i.e., a scale-dependent directional roughness index.The variogram was calculated along the RX and RY directions.The variogram calculated along RY (Fig. 6) can be viewed as a proxy for the hydraulic resistance (higher variogram values indicate higher flow resistance).The variogram calculated in the directions orthogonal to the river flow (RX) is intended to give information on the transverse organisation of the sediment.The calculation of the indices was performed via ad hoc routines developed inside the open-source statistical programming and graphics environment R (R Development Core Team, 2009).For variogram calculation, some of the functions of the Gstat R packages were used (Pebesma et al., 2004).
A great improvement in the readability of the directional variograms index can be achieved by calculating local anomalies (A) as the difference between the variograms calculated in small and large windows (a method analogous to the sharpening technique applied in image processing): In Eq. ( 2), the value of the variogram on LW is calculated from the original data, but similar results can be obtained (the differences are related to the different cardinalities of data encountered inside the moving windows) by calculating it, by means of moving averages, directly from the time series of the variogram calculated on SW, which represents the signal at the highest possible resolution.The interpretation of the local anomaly is straightforward: it determines whether the local variability on SW is higher or lower than the mean variability calculated on LW.This approach can be generalised to a multiresolution analysis by increasing the size of the windows.If S(t n ) is the discrete time series representing the analysed signal (here the variogram calculated on SW for a lag of 0.5 m), defined at times t n (t n = nδ; n=1, 2, 3. . .N; where δ is the time step, and N the number of time steps), the local anomaly at the level of resolution l, A l (t n ), can be derived as: where m(λ l ,S(t n )) indicates the moving average process of the signal S(t n ) centred at every time step t n within moving windows of size λ l , where λ l = λ 1 ,λ 2 ,λ 3 ,. . ., λ L indicate windows of increasing size up to the level of resolution L.
According to this formulation, the anomalies are calculated at different scales: by increasing the size of the windows it is possible to analyse the characteristics of the signal at lower frequencies.This approach gives similar results, at least for the datasets analysed in our case studies, to an analysis performed via wavelets (Lark and Webster, 1999;Percival and Walden, 2000), and, thanks to its simplicity, could be preferred to wavelets for an explorative analysis of data.For other tasks, such as image segmentation, signal compression and pattern recognition, wavelets analysis, after a careful choice of wavelet type and calculation parameters, is necessary.
In this case, to get a more smoothed signal, we used a modified version of Eq. ( 3), practically calculating the means at a given level l from the means of level l − 1 and not from the original signal: with m 0 (t n ) = S(t n )

Field surveys
Field surveys were carried out with two main objectives: -to revise the previous classification of channel bed morphology (Billi et al., 1998, Lenzi, 2001, Cavalli et al., 2008) for achieving a more detailed identification of morphological sequences; -to check the correspondence between the results of digital terrain analysis and ground truth.
The revised classification of channel-bed morphology of the Rio Cordon has been presented in the description of the study area (Figs. 2 and 3).The field survey aimed at evaluating the capability of the morphometric indices to depict morphological features of the channel bed was conducted using mobile GIS software (ESRI ArcPad 7.1) installed on a rugged tablet PC with GPS.The thematic layers of the GIS were as follows: -Orthophoto at a high resolution (0.15 m); -Shaded relief maps of the DTM; -Residuals of the interpolated surface; -Contour lines.
Longitudinal profiles of the calculated indices (directional variograms and local anomalies) were used for a field comparison with morphological features of the channel.
The mobile GIS integrated with GPS made it possible to identify in the field the exact position along the centreline in the transformed coordinate system.The survey of the main channel of the Rio Cordon allowed us to compare values of the proposed morphometric indices with actual channelbed morphology.Furthermore, the representativeness of the LiDAR-derived DTM was visually checked during field surveys, paying particular attention to the areas where a deterioration of LiDAR data quality or density could be expected (i.e., vegetation and pools with relatively deep water).

Results
Figure 7 presents the slope calculated on SW (Fig. 7a) and LW (Fig. 7b) along the channel profile.The patterns on the graphs of the SW and LW slopes are similar, with the strongest oscillations of slopes located in the boulder cascade reach and the lowest at the rapid with step reach.Obviously, the spatial variability of slope and its range of variations are larger for SW, due to its capability to detect high detail features.Negative values are visible in the plot of SW slope.Negative slopes in step pool reaches correspond to a rise of the channel bed at the downstream end of pools; in the boulder cascade reach and the cascade reaches, these are mainly related to the presence of scattered large boulders.
The slope index is useful for a general characterisation of the channel profile and should be used preliminarily to the geostatistical indices, which focus more closely on channel bed roughness.
The directional variograms calculated along the RX and RY directions for a lag of 0.5 m are presented for SW (SWG) and LW (LWG) in Fig. 8.The overall patterns of the RX and RY variogram graphs are similar, although the RX variogram is generally more scattered, especially in the lower part of the channel.This could be related to the fact that the transverse organisation of the sediments, also given the limited width of the stream, is not clearly distinctive for the different morphologies.Accordingly, the discussion will largely focus on the analysis of variogram-based indices calculated along RY (Fig. 8c and  d).
-In the upper part of the profile (from 0 to 250 m), which mostly corresponds to the boulder cascade reach, the variability is very high.The highest values are observed in the upper 150 m, where the Rio Cordon flows between steep rocky slopes.Here large boulders accumulated by rockfalls cause a high roughness of the channelbed surface.
-From 250 to 480 m, the roughness of the step pool and cascade reaches is markedly lower than in the upstream boulder cascade reach.
-In correspondence to the rapid with steps reach and part of the cascades reaches from 480 to 630 m, the variability is very low.
-From 630 to 750 m, the variability increases again, with the presence of definite spikes (cascade and step-pool reaches).
-From 750 to 980 m, the graph appears less variable (in the same zone the RX variogram graph is still quite erratic); this happens in correspondence to a quite long cascade reach.
-In the remaining part of the stream profile, a succession of high and low variability zones is present, with clear spikes related to the step-pool sequences, and to the presence of large boulders.In this section, an alternating sequence of step-pool and cascade reaches is present, with a morphology strongly influenced by lateral tributaries and the deposits of debris-flows.
Figure 9 presents multiscale local anomalies (Eq.4) of RY together with the original signal (bottom, Fig. 9a), represented by the SW variogram along RY (SWG ry ), and the remaining smoothed signal (Fig. 9f), in a format analogous to the one used in wavelet analysis.The main patterns visible in the anomaly graphs persist along the first three levels of resolution, while in the fourth level the structure is quite simplified.The smoothed signal (Fig. 9f, moving window of 80 m), which can be interpreted as a simplification of the original signal (Fig. 9a), outlines the main patterns just described for SWG ry .In the present study because of our interest in characterising the morphology at the finest possible scale, we are especially interested in the anomalies at the highest level of resolution (Fig. 9b, anomaly A1).Even if the interpretation of the signal is similar to that of the variogram on RY computed on SW in the anomaly plot, the zones of different variability are more distinct.Peculiar patterns and values observed in the plot of anomalies at high resolution, in combination with slope and variogram indices, have been checked in the field, looking for corresponding features in channel bed morphology (Fig. 10).
Figure 10a depicts a detail of the upstream part of Rio Cordon, classified as a boulder cascade, which is characterised by a scattered pattern with a wide range of variation and the presence of marked positive spikes.In Fig. 10b, two moderate negative peaks separated by a strong positive peak represent, respectively, the smoothed shapes of two pools and the high variability of the step.Figure 10c  symmetric distribution of positive and negative anomalies in the reach classified as a rapid with steps.Here steps built by relatively small grain size material are very low and short: this explains the plane trend in the anomaly graph.The highest positive peaks in the lower part of the channel are indeed related to significant structural roughness features, such as boulder clusters and big steps (Fig. 10d and e).In Fig. 10d, the isolated strong positive spike represents the big step visible in the photo.In Fig. 10e, an isolated step pool appears as negative peaks divided by a strong positive peak.
The pattern and range of variation of the anomalies represent the different character of the analysed stream portion and can detect particular morphological units, like step pool sequences shown as negative peaks divided by a strong positive peak.In the Rio Cordon, the pools, highlighted by negative anomalies, can be considered as the more distinctive features between the different channel bed morphologies, whereas clear positive anomalies can be found in all morphologies, but rapid with steps.
The field observations summarised above show that the index is also quite robust in the presence of an articulate morphology.A previous study on the channel bed of the Rio Cordon (Cavalli et al., 2008) analysed the main morphological classes by means of morphological indices of roughness based on LiDAR-derived DTM.The analysis was conducted by pooling all channel reaches belonging to the same morphological class.The comparison between different morphological classes demonstrated that roughness indices of the channel bed surface have significantly different values for different morphologies.In the present study, a geostatistical approach has permitted an enhanced representation of the spatial variability of channel bed morphology, allowing the recognition of even single morphological units.

Conclusions
The main outcomes of the geomorphometric analysis of channel bed morphology of the Rio Cordon are reported below.
-Working in a river-centred coordinate system was very useful for both the interpolation procedure and the geomorphometric study.according to a 1-D data representation, easily interpretable by means of various graphical, statistical and analytical tools.
-A slope index computed for the whole width of the channel bed using moving windows of different sizes makes it possible to depict the variability of slope along the channel profile, including the presence of negative slopes, mostly associated with pools in step-pool channel reaches.correspondence between values of local anomalies and particular features of channel morphologies, such as pools, large steps, clusters of boulders and pools.These morphological units, although also visible in topographic profiles from high-resolution data (Cavalli et al., 2008), need directional indices of topographic variability, such as directional variograms and local anomalies, to be adequately emphasised.
-Boulder cascades and rapids with step reaches can be easily recognised in graphs of longitudinal variograms and local anomalies, whereas other classes, such as step pools and cascades, show similar general patterns.Nevertheless, well-defined pools between steps which represent distinctive features between step pool and cascade morphologies, can be easily recognised as negative values in the graph of local anomalies (Fig. 10b and  e).The indices proposed in this paper are thus effective for recognising single morphological units (Fig. 10), whereas further procedures are needed for the automatic or semi-automatic classification of the morphology at the channel reach scale.
The Rio Cordon is a challenging case for the geomorphometric study of streambeds because of the presence of many short channel reaches with different characteristics (Fig. 2) and fuzzy transitions between the various morphologies.Moreover, the narrow channel width and DTM resolution limit the calculation of variogram to the shortest lag, leading to results similar to those achievable by means of directional derivatives.As a consequence, the potentials of the proposed geostatistical indices are not fully exploited.However, directional derivatives are less general and scalable than the variogram (or other spatial continuity indices such as the more robust madogram), which could be fully exploited in other morphological contexts (e.g., larger streambeds) and with higher resolution datasets acquired, for example, using terrestrial LiDAR.In many alpine headwater streams the morphological conditions of the channel are complex, like those of the Rio Cordon, and topographic datasets from aerial LiDAR, similar to that used in this study, are increasingly available.For these reasons, the Rio Cordon can be deemed as representative of achievements and problems that may occur in the geomorphometric analysis of channel bed morphology of small alpine streams by means of LiDAR data.

Figure 1 .Fig. 1 .
Figure 1.Geographical location and shaded relief map of the Rio C 5 outlines the studied portion of the channel.6

Figure 2 . 4 Figure 3 . 6 Fig. 2 .
Figure 2. Longitudinal profile and plan view of the studied channel reach, showing the 2 classification of bed morphology.3

Figure 2 .
Figure 2. Longitudinal profile and plan view of the studied channel reach, showing the classification of bed morphology.

Figure 3 .
Figure 3. Picture of the rapid with steps reach.

Fig. 3 .
Fig. 3. Picture of the rapid with steps reach.

Figure 4 .Fig. 4 .
Figure 4. Elevation of LiDAR points in a sample area of the study domain, plotted along the 2 centreline of the channel (a).Non-ground points (highlighted in a and in b), which correspond 3 to trees (see orthophoto in b), are easily recognizable, and can be expunged from the dataset.4

Figure 7 .Figure 7 .
Figure 6.a) the SW (blue) and LW (red) moving windows are superposed 2 residuals.b) a sketch of the computation of directional variograms via 3 presented.4

Fig. 6 .
Fig. 6.(a) the SW (blue) and LW (red) moving windows are superposed on the surface of residuals.(b) a sketch of the computation of directional variograms via SW windows is presented.

25Figure 6 Figure 7 .
Figure 6.a) the SW (blue) and LW (red) moving windows are superposed on the surface of 2 residuals.b) a sketch of the computation of directional variograms via SW windows is 3 presented.4

Fig. 7 .
Fig. 7. Longitudinal slope [m/m] calculated on SW and LW, plotted along the centreline of the streambed.

Figure 8 .Fig. 8 .
Figure9presents multiscale local anomalies (Eq.4) of RY together with the original signal (bottom, Fig.9a), represented by the SW variogram along RY (SWG ry ), and the remaining smoothed signal (Fig.9f), in a format analogous to the one used in wavelet analysis.The main patterns visible in the anomaly graphs persist along the first three levels of resolution, while in the fourth level the structure is quite simplified.The smoothed signal (Fig.9f, moving window of 80 m), which can be interpreted as a simplification of the original signal (Fig.9a), outlines the main patterns just described for SWG ry .In the present study because of our interest in characterising the morphology at the finest possible scale, we are especially interested in the anomalies at the highest level of resolution (Fig.9b, anomaly A1).Even if the interpretation of the signal is similar to that of the variogram on RY computed on SW in the anomaly plot, the zones of different variability are more distinct.Peculiar patterns and values observed in the plot of anomalies at high resolution, in combination with slope and variogram indices, have been checked in the field, looking for corresponding features in channel bed morphology (Fig.10).Figure10adepicts a detail of the upstream part of Rio Cordon, classified as a boulder cascade, which is characterised by a scattered pattern with a wide range of variation and the presence of marked positive spikes.In Fig.10b, two moderate negative peaks separated by a strong positive peak represent, respectively, the smoothed shapes of two pools and the high variability of the step.Figure10cshows a pattern of anomalies characterised by small variations and by a quite

1Figure 9 .Fig. 9 .
Figure 9. Anomalies of the RY variogram for a lag of 0.5 m, for four levels of spatial 2 resolution.The results are presented in analogy to wavelets analysis.a) Original signal (i.e., 3 longitudinal variogram for a lag of 0.5 m).b) to e) Multiscale anomalies for decreasing levels 4 of resolution calculated according to eq. (4).The sizes of the windows are indicated in the 5 figure.f) Remaining smoothed signal.6

-Fig. 10 .
Fig. 10.Comparison of local anomalies of the RY variogram for a lag of 0.5 m with field observations of the channel bed morphology for some portions (highlighted by boxes along the profile).Boulder cascade (a), step-pool sequence (b); rapid with step (c); isolated big step (d); and isolated step pool (e).