Abstract
This paper describes the methods and results of an investigation into the use of LIDAR1 derived digital elevation models as terrain input to the MIKE212 model. The purpose of this study is an attempt to utilize new developments in remote sensing technologies relating to elevation model derivation and implement these new sources in a GIS3 interfaced hydraulic modeling environment. The feasibility of simulating canal breaches is explored in areas near residential development.
Ortho-rectified photographic images are referenced to eliminate LIDAR signals falling within building footprints and a 'bare surface' elevation model is interpolated. The buildings are subsequently re-extruded in the raster domain utilizing a GIS to provide realistic structural definitions. These data are converted into a MIKE21 finite difference mesh and the numeric modeling is performed. Snapshots at desired time steps and animations are generated in both the GIS and MIKE21 software. ArcView 3D Analyst is used to assist visualization. The results of these studies are supplied to local authorities for emergency response planning and to evaluate the degree of risk posed to local residents.
KEYWORDS: GIS, LIDAR, MIKE21, HYDRAULIC, ARCVIEW
Introduction
In 1999, The U.S. Bureau of Reclamation (USBR) started a process to procure and evaluate LIDAR data in a number of areas about the United States. Although these data were collected for numerous applications, some of the project areas focused on the feasibility of utilizing these data as elevation sources for low flow hydraulic analysis.
An initial project area was set up in the Provo, Utah area where a USBR canal meanders through a populated area. Recently, in Ogden, Utah, a canal with similar characteristics breached and caused considerable damage. It was desired to determine if it is possible to simulate low flow canal breach scenarios. To model this would require high-resolution digital elevation models. LIDAR technology appeared to be a potentially viable solution, as time constraints existed and internal resources to perform stereo photogrammetry were impacted.
The study area is a 20 mile reach following the Provo Reservoir Canal, and extending away from the canal in the down slope direction for a distance of one mile. Several locations of interest were chosen to perform canal breach simulations. The selection of these locations was based on residential development density and overall condition of the canal embankment at those spots. Typical flow in the canal is 350 ft3/sec. The maximum design capacity of the canal is about 500 ft3/sec.
View of one of the selected sub study areas for breach simulation
Since 1990, USBR has been involved with integrating Spatial Analysis
Systems (GIS's) with numeric flow models. This project was designed using
the GIS as a data manager, manipulator, and analyzer for both pre- and
post-processing. The DHI Water and Environment's two-dimensional (2D) MIKE21
hydrodynamic finite difference model was implemented as the flow simulator.
Collection of LIDAR Data
The laser system is mounted on an aircraft, along with a Global Positioning System (GPS) and an Inertial Measuring Unit (IMU). The GPS enables derivation of the laser's latitude, longitude, and ellipsoid height. The IMU provides information on the aircraft's roll, pitch and yaw. Using these devices, a computer can calculate the position of the laser as a function of time.
As the aircraft proceeds along the flight path, the laser oscillates back and forth perpendicular to the aircraft's direction, while rapidly sending and receiving laser pulses that reflect off of the earth's surface. Utilizing the information on the position and attitude of the sensor, the elapsed time between laser pulse and sensor retrieval, and the fact that the speed of light is constant and known, a large series of x, y, and z values are collected.
Typically, the laser is flown at altitudes of 5,000 to 8,000 ft above the ground surface. Theoretically this can produce a horizontal accuracy of +/-0.4 meters and a vertical accuracy of +/- 0.15 meters.
Illustration of airplane, GPS units, satellites, and laser
View
of LIDAR point sampling density over photography
Analysis of LIDAR Data
The monochromatic beam that leaves the laser has a diameter of about .3 inches, however, by the time it strikes a surface feature, the beam has dilated to a diameter of 18 to 20 inches. This beam dilation can cause multiple returns: some of the photons in the pulse may bounce off of one object, and the rest may bounce off of another. A building is a good example- half of the laser pulse may hit the top of the structure and return with a given lag time, and the other half may miss the structure completely and hit the ground, producing a longer lag time.
The existence of multiple returns requires editing the LIDAR data to
determine which returns are to be used. Since we are flow modeling, we
are interested in vegetation removal and are usually interested in the
last return, as this represents the return closest to the ground, if not
the ground. The data are subsequently processed using 'vegetation removal
software' (which is usually proprietary) in an attempt to remove vegetation
and structures. The data are then loaded on to a softcopy photogrammetry
workstation to edit the remains of the vegetation and buildings. This processing
was performed by Horizons Technology.
Example of raw LIDAR data re-sampled at 1 meter from TIN data structure
Zoom of raw LIDAR data- note roof pitches and vegetation
Consecutive zooms of LIDAR data after vegetation removal filters
Spatial
Analyst view of vegetation and buildings removed from LIDAR surface
Bringing Data into the GIS
The LIDAR data received from Horizons Technology4 came in the form of large ASCII files containing NAD83 coordinates and elevations in meters. Both the 'raw' unfiltered points and the 'bare earth' points (vegetation and structures removed) were obtained.
LIDAR data sets can become very large, typically millions of points per flight swath. This becomes problematic when using ARC/INFO on current Operating System capabilities. We found that it is better to use non-GIS data filtering techniques to extract the areas of interest and subsequently use the GIS to operate on smaller data sets. Filters were written in C program language to accomplish this.
Sometimes a study area contained multiple LIDAR swaths. It is advantageous to combine the points into one file, rather than to process the files individually and later merge them on the GIS. This eliminates problems with edge effects.
The building and vegetation removal algorithms essentially eliminated those points that were suspect. We found that these procedures are about 90% effective, and further editing is required using Ortho-rectified photography if additional refinement is needed. The result of this is a data set containing 'no LIDAR points' where the vegetation and buildings existed, and interpolation across these areas is performed.
The data were then loaded into ArcView as a Table, which was in turn used to create an Event Theme. The Spatial Analyst extension facilitated the creation of a TIN data structure, which was then converted to a GRID5 at the desired sampled cell size (two meters in our case).
At this point, we had two GRID structures: one 'raw' data containing buildings and vegetation, and a 'bare earth' data with the buildings and vegetation mostly removed. We have used two methods to replace the buildings back on the surface:
1. Use Ortho-rectified photography to create a polygon coverage of the structures. This is converted to a GRID, where a Map Algebra CON statement is used to extrude the buildings back on the 'bare earth' elevation model.
2. Use the raw LIDAR GRID and the 'bare surface' GRID to create a difference GRID by subtracting the two. Since buildings are higher and have solid returns, Map Algebra statements can be written against the difference GRID to extrude the buildings back on the 'bare earth' surface elevation model. Additional clean up can then be performed using the Ortho-photography.
While procedure (1) is philosophically more clean, procedure (2) is faster in practice for our purpose, which is to get a general idea of what will happen in these events. The elevation lattice is manipulated until it is considered satisfactory.
Ortho-rectified image of a modeled area
Arcinfo GRID representing difference between raw and cleaned up LIDAR data
ArcMap view of LIDAR hillshade with pseudo-transparent ortho-rectified image
Spatial
Analyst view of LIDAR surface after buildings are extruded
MIKE21 Finite Difference Elevation Arrays and Boundary Conditions
MIKE21 is a two dimensional, finite difference hydrodynamic flow model available from the DHI Water & Environment. We began using a two-dimensional flow model when we realized that one-dimensional models were not adequate to describe certain types of terrain, such as flood plain areas, overland topography, split flow scenarios, etc.
The Arc/Info elevation GRID was converted to an ASCII array with an information header using the GRIDASCII6 command. This file was converted to a MIKE21 array file (T2 file) using a C filter specifically designed for elevation array translations between Arc/Info GRID and MIKE21. The elevation array can then be visualized and manipulated within the MIKE21 T2View Software.
A flat input platform is placed on the elevation lattice where the flow is to enter the model. One end of the platform is open to the elevation lattice and the other end is where the variable flux boundary condition is placed (flows in MIKE21 are expressed as flux densities- m3/s/m width)7. The input platform facilitates a stable flow pattern into the model.
Input
platform for hydrograph inflow (right center); houses are red
Problems with Modeling Overland Flow
The modeling of water storage or conveying facilities that breach poses unique problems, one of which is flow over otherwise dry land. In a 'normal' one dimensional flow model, all of the finite difference elements are wetted with some water depth, and the finite difference equations are solved for this contiguous set of water points. This is not so for a two-dimensional overland flow model. At the beginning of a run, almost all of the cells are dry.
From a philosophical standpoint, this problem partially manifests itself due to the fact that when finite difference techniques are implemented, and thereby space is descretized, one is forced into the fact that either an element (or cell) is 'wet' (part of the solution domain) or it is not (a boundary). Analytically, another way to view this is that the boundary condition is part of the problem- it is not fixed. Compounding this problem is the fact that Saint Venant equations do not describe flow regimes at the leading edge of a wave as well as could be desired.
Facing these limitations, one can either set all the cells initially
with a small water depth, or implement some type of wetting procedure where
a cell is wetted as a function of neighboring cell fluid velocities and/or
water depths. Both approaches posed their own set of problems. Currently
with MIKE21 we are using a wetting approach, whereby a cell is wetted when
any of its neighboring cells attains a small, user set water level. The
idea here is that the Saint Venant equations 'apply' between the newly
wetted cell and its neighbors as soon as it is wetted, and a cell will
be wetted as soon as it realistically can be. DHI is currently researching
methods to address this area.
Other Modeling Considerations
When using the MIKE21 model, a dimensionless flow parameter comes into play called the Courant number, defined as C = c(D t/D x) where D t is the time step, D x is the grid spacing, and c is the celerity, defined by c = sqrt(g × h) where g is the gravitational constant and h is the average water depth. In cases where high current velocities may exceed the celerity, the maximum current speed, or Umax, may be more appropriate to use in the calculation. The Courant number describes how many grid points the celerity 'information' propagates within one time step. In numerical models using an explicit technique to advance the solution in time, the Courant number must be less than unity. Implicit methods allow for higher Courant numbers. In MIKE21 when modeling flow, the Courant number should be kept under 5 (in general). Because of the small cell size (D x = 2 meters), the time step (D t) must be set quite small to meet the Courant condition. The constraint on D t was usually around 0.1 second or less for a 2-meter grid spacing.
Two-dimensional models do not simulate vertical velocity profiles, and therefore inherently cannot fully describe supercritical flow, which is a 3 dimensional phenomenon and may really require solution of the full Navier-Stokes equations. In these situations, usually some technique is used to dissipate the energy of the flow in an attempt to restrain it to the sub-critical domain. MIKE21, as of yet, does not implement automated methods to dissipate energy when the Froude8 numbers approach unity (signifying supercritical flow), although in the next release these enhancements will be included. As such, manual-damping methods can be used in those situations where the flow becomes supercritical. Modification of the eddy-viscosity map in the suspect areas usually produces the desired correction.
Three-dimensional effects dominate the leading edge or front of a flood
wave, and these regions are characterized by substantial turbulence, which
implies considerable energy dissipation. We are not actually modeling these
effects, and it is uncertain how much error this introduces to what would
be the 'actual' propagation of the wave front. This is a considerable concern,
because the leading edge as a function of time is one of the main factors
used in emergency management planning.
Results and Conclusions
A GIS using the Esri suite of software running on UNIX/NT platforms was interfaced to DHI Water & Environment's MIKE21 two-dimensional hydraulic flow model. LIDAR data were obtained for a study area, manipulated in the GIS, and translated into the MIKE21 modeling environment. Different scenarios were evaluated, and the results were translated back into the GIS for visualization and analysis.
Various steady state flow rates were introduced onto the elevation lattice. The cell size was set to 2 meters. At a breach flow of 300 ft3/s (steady state), the water depths ranged up to 2 meters and the velocities up to 7 ft/s. The horizontal velocity profiles that were set up were within reasonable expectations. Model run times varied from 25 to 100 minutes. Future studies will include unsteady breach flows with varied duration.
It was clear from the results that a two-dimensional flow simulator was required to adequately describe the flow. MIKE21 seemed to function adequately once the Courant restraint was honored. A one-dimensional model would not have been able to model these scenarios properly. Translators between MIKE21 and Arc/Info supplied by DHI worked very well. Arc/Info GRID and Arc/View Spatial and 3D Analyst were the main GIS components used.
There were no real means to calibrate the simulations. Usually when
these types of events occur, flow and stage data are not recorded. Since
we are mandated to perform such analyses, we have tried to implement the
best available modeling and GIS technology available to us, tempered of
course by cost and resources. The modeling methods described in this paper
appear to give us a reasonable description of what may occur in these low
flow breach scenarios, and USBR will be applying this technology toward
future emergency management planning endeavors.
Maximum inundation - 300 ft3/s run over ortho-rectified image
MIKE21 T2View of depths and velocity vectors of flow around a house
ArcMap view of LIDAR lattice with maximum inundation area
A scenario showing maximum velocities attained for that run (50 minutes)
Spatial Analyst view (1) of maximum inundation on LIDAR surface at 300 ft3/s
Spatial Analyst view (2) of maximum inundation on LIDAR surface at 300 ft3s
Spatial Analyst view (3) of maximum inundation of LIDAR surface at 300 ft3/s
Animation of 300 ft3/s steady state flow onto LIDAR lattice
Animation
of 300 ft3/s steady state flow into LIDAR lattice (zoom)
Footnotes
1 Light Detection and Ranging
2 The DHI Water & Environment's two-dimensional flow model
3 Geographic Information System
4 Horizons Incorporated, Rapid City, South Dakota
5 We are referring to an Arc/Info GRID data structure
6 An Arc/Info command that converts a GRID to an ASCII array file.
7 This is the volumetric flow rate per meter perpendicular to flow
8 The Froude number is an indicator of the
type of flow present- subcritical, critical, or supercritical
Authors
Tom Heinzer
GIScience/Chemical Engineer
U.S. Bureau of Reclamation
theinzer@mp.usbr.gov
Phone: (916) 978-5273
Michael Sebhat
MPGIS Service Center Manager
GIScience/Electrical Engineer
U.S. Bureau of Reclamation
msebhat@mp.usbr.gov
Phone: (916) 978-5272
Bruce Feinberg, P.E.
Civil/GIScience Engineer
Denver Technical Services Center
U.S. Bureau of Reclamation
bfeinberg@do.usbr.gov
Phone: (303) 445-3241
Dale Kerper, M.S.C.E
Senior Engineer
DHI, Inc.
Encintas, CA
dkerper@home.com
