Floodplain Morphometry Extraction from a High Resolution Digital Elevation Model: a Simple Algorithm for Regional Analysis Studies

B. A. Dodov and E. Foufoula-Georgiou

AbstractIn this paper we introduce a fast algorithm for floodplain delineation and morphometric
analysis over large river basins. The underlying assumption of the delineation procedure is that the depth of inundation used to define the lateral extent of the floodplain is proportional to the depth of the channel at bankfull flow, a relationship which emerged from a detailed and objective regional analysis. Channel bankfull depth was obtained from a regional geomorphologic analysis and was found to follow approximately a power law with contributing area. The regional analysis of bankfull geometry was combined with traditional analysis of channel network topology from Digital Elevation Model (DEM) and high resolution hydrography data (NHD data) to extract the floodplain mask by means of a scaledependent kernel function (scale is used as synonymous to the contributing area draining to a point in the network). After the floodplain mask was extracted a medial axis decomposition procedure was used to compute the floodplain centerline (skeleton) and, respectively, the floodplain cross-sectional shape at any location in the river network. Analysis of floodplain geometry over a wide range of scales revealed important connections between fluvial and hydrologic processes and showed that the nature of these connections is scale-dependent. Index TermsDigital Elevation Model, Floodplain, Geomorphologic analysis, Hydrology.

I. INTRODUCTION Floodplain delineation procedures have been traditionally related to flood risk management and, as such, they require maximum precision in order to optimize the local insurance and reinsurance policies. The traditional floodplain delineation approach typically considers the floodplain inundated by a flood peak of particular magnitude corresponding to a given return period (i.e. the number of years within which such an event will be equaled or exceeded once, on the average) using one-dimensional river hydraulics models such as HEC-RAS (US Corp. of Engineers), MIKE11 (Danish Hydrologic Institute), SWMM-EXTRAN (US Environmental Protection Agency), etc. After the free water surface is obtained from a river hydraulics model, it is extrapolated (and/or interpolated) over the Digital Elevation Model (DEM) of the terrain and all submerged parts of the terrain surface are assumed to be part of the floodplain (e.g. see [9] for technical details). Since the simulation of stream hydraulics based on the full dynamic St. Venant equations is computationally expensive and a data-demanding process even in the one-dimensional case, hydraulic simulations (and respectively floodplain delineation) are usually performed with respect to separate river reaches or part of the network but not with respect to the whole stream network.

In contrast, the purpose of a regional geomorphologic analysis is to: (i) define the extent of the natural variability envelope of channel and floodplain morphologies (both cross-sectional and planform) as they vary with contributing area, and (ii) quantify the multidimensional statistical dependencies between these morphologies over a large range of scales within regions of similar climatic and geologic conditions. Immediately, it becomes clear that the implementation of the traditional floodplain delineation approach would be infeasible for such an analysis as one has to consider scales (contributing areas) from one to hundreds of thousands of square kilometers. In addition, since in a regional geomorphologic analysis the emphasis is on the statistical properties of the factors describing channel and floodplain morphologies (i.e. distributions, moments and covariance structure) the precision of the floodplain determination is a second order concern. Finally, in a regional geomorphologic analysis one is interested in floodplain morphology defined within the floodplains natural boundaries (which in most cases are well definable based on geomorphological grounds) rather than in the part of the floodplain submerged during an event of particular magnitude. The outcome of the above considerations is that for regional geomorphologic analysis one needs a floodplain delineation procedure which can be applied over large river networks (say of the order of 105 km2) and which incorporates some geomorphologic criteria that would allow the determination of the floodplain up to its natural boundaries. In this paper, we propose an automated floodplain delineation procedure based on the observation that the depth of inundation ( D fpin ) used to define the lateral extent of the floodplain is related to the depth of the channel at bankfull flow ( Dbf ). Through extensive geomorphologic analysis we show that this relationship is linear with a proportionality coefficient 0.6, i.e. D fpin = 0.6Dbf . In addition, we provide an automated procedure for floodplain centerline extraction and cross-section determination at any point within the network, which allows a continuous analysis of floodplain geometry along the river network. The paper is organized as follows. In section 2, we introduce the regional geomorphologic analysis necessary for the determination of floodplain inundation depth as a function of scale. In section 3, the floodplain delineation algorithm is introduced and some examples from basins in central USA are shown. A centerline extraction algorithm is provided, based on the so-called median axis transform procedure, by means of which a continuous analysis of floodplain properties can be performed along the channel network. The importance of the proposed analysis is demonstrated in section 4 by showing a comparison of floodplain and hydrologic data analysis from Oklahoma-Kansas and from the Appalachian Piedmont regions. Final comments and conclusions are given in section II. REGIONAL GEOMORPHOLOGIC ANALYSIS: CHANNEL BANKFULL AND FLOODPLAIN INUNDATION DEPTHS AS A FUNCTION OF SCALE The floodplain can be defined as a flat area adjoining a river channel constructed by the river in its present condition and climate [4]. As the floodplain is essentially formed by lateral migration and controlled by the water and sediment transport in a river, it is reasonable to assume that the channel and floodplain geometries (and particularly depths) are closely related. With an increase of contributing area and channel forming (bankfull) discharge, thin and narrow floodplains are replaced by wider and flatter floodplains of increasing thickness and inundation. Deeper inundation in turn provides more overhead

cover and greater channel depth. Intuitively, this line of thinking suggests a close connection between channel depth at bankfull and floodplain inundation depth. Although such an intuitive connection is in general accepted by geomorphologists (personal communication with Gary Parker and Chris Paola, University of Minnesota) to the best of our knowledge it has never been supported on physical grounds or by a geomorphologic analysis. Below, we provide such an analysis showing not only a close connection between bankfull channel and floodplain inundation depths, but also that the emergent relationship is linear to a good approximation over a range of scales. Before continuing with our analysis, we first need to establish a connection between channel bankfull depth at a given point in the network and the contributing area draining to this point. It has been shown that: (i) channel depth at bankfull follows a simple power-law relationship with the discharge corresponding to a given frequency of exceedance (e.g. see [1,7]) , and (ii) discharges of specific return period follow power-law relationships with contributing area [6,7]. A reasonable expectation is then to look for power-law relationships between contributing area (scale) and channel depth at bankfull. For the extraction of such a relationship we analyzed independent measurements of stage, maximum width, depth and discharge under different flow conditions (up to several hundred measurements per station) from 113 USGS gauging stations in Oklahoma and Kansas (see Fig. 1). Bankfull conditions were considered to occur when a break in slope was observed in both stage-discharge and width-discharge relationships (a transition to a very slow increase in stage and very rapid increase in width after bankfull discharge is reached). For small contributing areas, bankfull geometry derived from field observations was added to increase the reliability of our analysis [5,8]. In Fig. 2a we plot channel depths at bankfull as a function of contributing area. As can be seen, the relationship can be well approximated by a log-log linear relationship (solid line) with a break at approximately 700 km2. As we will see later on, the break at 700 km2 is related to a break in floodplain geometry occurring at the same scale, and also a break in flood statistics suggesting an intimate connection between fluvial and hydrologic processes that change their nature with scale (see also [3]). Having established the scale-dependence of bankfull channel depth, we proceed with the analysis of the connection between the depth of a channel at bankfull and the inundation depth of the floodplain adjacent to the channel. For this analysis, we apply the floodplain delineation procedure described in detail in the next section for four patches (shown in Fig. 1) of the DEM of Neosho River basin for 25 loguniform scale intervals between 1 and 15,000 km2 and log-uniform inundation intervals between 0.1Dbf and 3Dbf . In Fig. 3 we plot the relative floodplain area (RFA: area covered by the water surface corresponding to a particular inundation depth divided by the area corresponding to the maximum depth of inundation, i.e. 3Dbf) as a function of inundation depth for two scale intervals log-centered at approximately 9 and 9000 km2. A break is observed in each of these relationships denoting that below and above the break the rate of increase of the water surface for a given rate of increase of the inundation depth changes. The depths at which these breaks occur were interpreted as the inundation depths defining the lateral extent of the floodplain at a given range of scales. Curves such as that of Fig. 3 were examined for all 25 scale intervals and the inundation depths corresponding to the breaks were extracted. These depths were then plotted in Fig. 2a as a function of scale (solid circles) and compared with the bankfull depths. It is interesting to observe that a linear relationship between the two depths emerged with a

coefficient of proportionality p = 0.6, i.e, D fpin = 0.6Dbf . The above analysis puts on solid grounds the commonly used rule of thumb that the floodplain inundation depth is a percentage of the bankfull depth and provides an objective way of estimating the proportionality coefficient for the region of interest. This relationship was adopted in our floodplain extraction analysis explained in detail in the next section. III. FLOODPLAIN DELINEATION AND ANALYSIS A. Floodplain Delineation The idea of the floodplain delineation algorithm is to locally fill the DEM up to the depth of the inundation of the floodplain using a horizontal disk centered at pixels along the river network. The radius of the disk was chosen proportional to the median radius of channel planform curvature at any given scale. The relationship between channel planform curvature and contributing area was extracted from high resolution hydrography data of the Neosho River basin shown in Fig. 1 (USGS National Hydrography Dataset) and approximated by a power-law (see Fig. 2b). The median radius of curvature criterion was chosen because this parameter was considered the one which is most related to the width of the floodplain. Since the computed extend of the floodplain is much less sensitive to this parameter compared to the depth of inundation, only a trial and error procedure was applied for its determination (essentially, the only requirement in this case is that the diameter of the disc function has to be larger than the width of the floodplain at a given scale). All the pixels covered by the disk and with elevations less tha that of the disk were assumed to belong to the floodplain. More specifically, the floodplain delineation algorithm consists of the following: a) Input data: DEM grid, Area Accumulation Grid (AAG), Downstream Slope Grid (DSG). b) Choice of an appropriate number of scale ranges. Areas less than 1 km2 should not be considered since the resolution of the DEM does not allow proper treatment of floodplain cross-sectional shapes for such contributing areas. c) Determination of the pixels assigned to the floodplain: for each range , { define the set AAG , (the set of all pixels with area within the range ) for each Pixi DEM , { DEM define the set Pix , (d , Pixi ) (all DEM pixels within a distance d from Pixi ). DAG DEM define Pix , = Pix , I AAG, DAG if not empty( Pix , ) {
i i i i

DAG ElevPixi = E[ Elev Pix j ], Pix j Pixi ,

if ( Elev Pix DEM Pix ) then { Pixi FP } else { Pixi FP }

}
i

where DEM Pix is the elevation of Pixi extracted from Pixi DEM, and FP is the set of all pixels

belonging to the floodplain with contributing areas within the range . By trial and error, optimal values for the distance d was found to be 6 times the median radius of curvature of the channel centerline for a given range . In cases where subsets of floodplain pixels were discontinued, the pixels representing the river network were used to keep the continuity of floodplain network (see Fig. 4 for two examples of the floodplain binary grid FBG mapped on the shaded relief of two regions representing parts of Neosho and Blue River basins in Oklahoma and Kansas, as shown in Fig. 1). B. Floodplain skeleton extraction For the extraction of the floodplain skeleton the so-called medial axis transform was used. The procedure consists in the following: for every pixel from the floodplain calculate the unit vector pointing to the nearest pixel on the floodplain boundary. Calculate the maximum angle between a given pixels unit vector and the unit vectors of the surrounding pixels. The skeleton is formed by pixels with maximum angles larger than / 6 (see Fig. 5 for the skeletons of regions 1 and 2 shown in Fig. 1 and 4). More specifically the steps in this procedure are: a) Compute the direction grid DG from FBG (i.e. from every pixel on the floodplain to the floodplain boundary). b) Derive the skeleton grid SKG from the direction grid DG using the following algorithm: for each Pixi DG {
Ang Pixi = max[sin( Pixi ) sin( Pix j ) + cos( Pixi ) cos( Pix j )] DG Pix j Pixi ( X Pixi , YPixi )

if ( Ang Pix / 6) then { Pixi SKG } else { Pixi SKG }

}
DG where SKG is the set of all pixels belonging to the floodplain skeleton and Pix is the subset of pixels surrounding a given pixel Pixi of DG . After the floodplain skeleton was derived, it was vectorized using the public-domain WINTOPO software and transformed back to raster format in order to: (i) thin the skeleton, (ii) clean all branches and noise pixels not connected to the floodplain network.
i

C. Floodplain centerline extraction In order to properly analyze the variation of floodplain width as a function of contributing area the branches within the floodplain have to be removed (e.g. pixels with contributing area 3 km2 cannot have 2 km wide floodplain). To do this, we take every pixel in a given range and compute its distance to the nearest pixel of range larger than (higher order pixel). If this distance is less than the distance from the higher order pixel to the floodplain boundary, the pixel of range is removed from the skeleton. The algorithm for this procedure is as follows: a) Compute the distance grid DISTG from SKG and FBG (i.e. the distance from every pixel on the skeleton to the floodplain boundary). b) Clean all the branches within the floodplain using the following algorithm:

for each Pixi SKG , { Dist ( Pixi , Pix j ) = min[ Dist( Pixi , Pixk )] find Pix j SKG , > such that Pixk SKG , > SKG if [ Dist( Pixi , Pix j ) > DISTGPix )] then { Pixi } else { Pixi SKG } }
j

To every pixel on the centerline, a contributing area, distance and direction to the floodplain boundary is assigned. IV. ON THE IMPORTANCE OF FLOODPLAIN GEOMETRY ANALYSIS OVER A LARGE RANGE OF SCALES: CONNECTION TO REGIONAL FLOOD FREQUENCY ANALYSIS Since every pixel of the floodplain centerline can be assigned a contributing area, a distance and a direction to the floodplain boundary, we can compute the floodplain transverse slope from (DEM) profiles and perform statistical analysis of floodplain properties for any range of scales. In Fig. 6 we plot the median floodplain half-width and median transverse slope as a function of scale (median because the distributions of the floodplain widths and transverse slopes are skewed) together with the coefficient of variation CV (standard deviation divided by the mean) of the maximum annual peak discharges for 72 out of the 113 stations shown in Fig. 1. Compared with Fig. 2a these plots clearly suggest a common break in channel/floodplain geometries at the scale of approximately 700 km2, which obviously affects the relative variability of flood peaks in a different way below and above the scale of the break. Although the physical explanation of the existence of the break goes beyond the scope of this paper, it is emphasized that the continuous analysis of floodplain geometry with scale reveals an important interplay between hydrologic and geomorphologic processes which changes character at the scale of approximately 700 km2 (for further details see [3]). To further support the importance of the floodplain geometry analysis on flood statistics let us consider another region with different geologic and climatic conditions, namely, the Appalachian region in the Eastern US (see Fig. 7). In their work Smith [10] and Gupta et al. [6] showed that the CV of maximum annual floods in this region increases up to approximately 100 km2 and then decreases. The result of our floodplain geometry analysis in Fig. 8, showing a break at 100 km2, once again supports the close connection between floodplain geometry and streamflow variability. It is noted that the scale of the break is different in the Oklahoma-Kansas and Appalachian regions as the physical mechanisms responsible for the feedbacks between fluvial and hydrologic processes differ significantly between these two regions. V. SUMMARY AND CONCLUSIONS This work proposes a simple and fast algorithm for floodplain delineation and analysis based on geomorphologic considerations, namely, on a documented from regional analysis linear relationship between the floodplain inundation depth and the depth of the channel at bankfull flow over a range of scales. A floodplain delineation algorithm is implemented by locally filling the DEM up to the depth of floodplain inundation using horizontal disks centered at pixels within the river network. The radius of the disk was chosen proportional to the median radius of channel planform curvature at any given scale. A

medial axis transform algorithm was used to extract the floodplain centerline and consequently floodplain cross-sectional geometry continuously along the river network. The floodplain delineation and analysis algorithms were tested at several basins in Central (Oklahoma, Kansas) and Eastern (Appalachian Region) United States. Combined with a regional flood frequency analysis, the analysis of floodplain geometry continuously along the river network revealed an important interplay between hydrologic and geomorphologic processes which changes character with scale. VI. ACKNOWLEDGMENT This research was jointly supported by the STC program of the National Science Foundation via the National Center for Earth-surface Dynamics (NCED) under agreement EAR-0120914, and by the Minnesota Supercomputing Institute (MSI). We are thankful to Gary Parker for his insightful comments and encouragement during the course of this study. REFERENCES [1] Dodov, B., and E. Foufoula-Georgiou, Generalized hydraulic geometry: Derivation based on a multiscaling formalism, Water Resour. Res., 40 (6), W06302, doi:10.1029/2003 WR002082, 2004. [2] Dodov, Boyko and E. Foufoula-Georgiou, Generalized hydraulic geometry: Insights based on fluvial instability analysis and a physical model, Water Resour. Res., 40(12), W12201, 10.1029/2004 WR003196, 2004. [3] Dodov, B., and E. Foufoula-Georgiou, Fluvial processes and streamflow variability: interplay in the scale-frequency continuum and implications for scaling, to appear, Water Resour. Res., 2005. [4] Dunne, T. and L.B. Leopold, Water in Environmental Planning. W.H., Freeman and Company, New York, NY, 1978. [5] Dutnell, R.C., Development of bankfull discharge and channel geometry relationships for natural channel design in Oklahoma using a fluvial geomorphic approach, M.Sc. Thesis, University of Oklahoma, 2000. [6] Gupta, V. K., O. J. Mesa, and D. Dawdy, Multiscaling theory of floods: Regional quantile analysis. Water Resour. Res., 30(12), pp. 3405-3421, 1994. [7] Leopold, L. B. and T. Maddock, The hydraulic geometry of stream channels and some physiographic implications, U.S Geol. Survey Prof. Paper, 252, 57p., 1953. [8] National Water Management Center, Regional Hydraulic Geometry Curves, NRCS, 2004 (URL: wmc.ar.nrcs.usda.gov/technical/ HHSWR/Geomorphic/osage.html). [9] Noman, N. S., E. J. Nelson, and A. K. Zundel, A Review of Automated Flood Plain Delineation from Digital Terrain Models, ASCE Journal of Water Resources Planning and Management, Vol. 127, No. 6, pp. 394-402, 2001. [10] Smith, J. A., Representation of basin scale in flood peak distributions, Water Resour. Res., 28(11), pp. 993-2999, 1992.

Fig. 1. Locations of the 113 stations , the shaded Neosho and Blue River basins, and the patches used for the determination of D fpin .

Fig. 2. a) Bankfull depth, Dbf versus contributing area for the 113 stations of Fig. 1 (computed from stage-discharge curves or field surveys). On the same figure, the flood inundation depth D fpin computed for the 4 patches in Fig. 1 is also shown. In is noted that the ratio D fpin / Dbf is almost constant and equal to 0.6; b) Median radius of curvature vs. scale (based on the analysis of high resolution hydrography data of the Neosho River basin).

Fig. 3. Determination of D fpin for two scale intervals centered at 9 and 9000 km2.

Floodplain mask

a) Region 1: Neosho River basin

Floodplain mask

b) Region 2: Blue River basin. Fig. 4. Delineated floodplain (mask) mapped on the shaded relief of regions 1 and 2 (respectively parts of Neosho and Blue River basins shown in Fig. 1)

Fig. 5. Floodplain skeletons extracted from the floodplain masks of Regions 1 and 2.

Fig. 6. Top plots: Floodplain geometry of Neosho River basin versus scale. Bottom plot: coefficient of variation (CV) of maximum annual peak discharges versus scale for 72 stations in Oklahoma and Kansas.

Fig. 7. Locations of the 104 stations of Smith, 1992 [4] and the region chosen for floodplain geometry analysis.

Fig. 8. Top plot: Floodplain half-width of the shaded region in Fig. 6; Bottom plot: CV of maximum annual peak discharges versus scale for the 104 stations of Smith, 1992 [4].