Jun 17, 2025
A watershed-scale approach to 2D flood modeling
- Meredith Swallom1,
- Jason Dortch1,
- Mackenzie Choffel2,
- Desiree Cunningham1,
- Evelyn Bibbins1
- 1Kentucky Geological Survey;
- 2University of Kentucky Earth and Environmental Sciences
Protocol Citation: Meredith Swallom, Jason Dortch, Mackenzie Choffel, Desiree Cunningham, Evelyn Bibbins 2025. A watershed-scale approach to 2D flood modeling. protocols.io https://dx.doi.org/10.17504/protocols.io.261ge8y9wg47/v1
Manuscript citation:
Swallom, M.L., Dortch, J.M., Thigpen, J.R. 2025. Assessing the effects of anthropogenesis on Appalachian flood severity: An eastern Kentucky case study. Geomorphology 482. https://doi.org/10.1016/j.geomorph.2025.109818.
License: This is an open access protocol distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Protocol status: Working
We use this protocol and it's working
Created: June 16, 2025
Last Modified: June 17, 2025
Protocol Integer ID: 220354
Keywords: flood model, hazard, Appalachia, lidar, scale approach to 2d flood modeling, assembling 2d flood model, 2d flood modeling, 2d flood model, scale approach, scale
Funders Acknowledgements:
National Science Foundation
Grant ID: EAR 2242120
National Science Foundation EPSCoR CLIMBS
Grant ID: OIA 2344533
Abstract
This document is meant to serve as a guideline for assembling 2D flood models at large (watershed) scales.
Materials
All software needs are detailed in the main protocol text. To ensure your hardware is capable of modeling, consult with the TUFLOW manual or TUFLOW support.
Introduction, Software, and Data Preparation
Software Selection and Installation
Two-dimensional (2D) modeling is useful in studies in which the spatial distribution of maximum flood height may be informative, as was the case with the July 2022 flood in eastern Kentucky. TUFLOW Heavily Parallelised Compute (HPC) was selected for 2D modeling as it allows GPU acceleration and implementation of subgrid sampling (SGS), which decrease modeling time. For full details on the SGS method, see the TUFLOW manual.
To install the latest version of TUFLOW HPC, visit Downloads, then download and install all programs listed under "Requirements and Downloads" on Tutorial Introduction - Tuflow.
Building Blocks of a Flood Model
A flood is, in simplest terms, a result of the interactions between the landscape, the hydrologic properties of that landscape, and water imparted onto the landscape. The three minimum key elements to consider when assessing flood likelihood or severity are thus: topography of the area, how water interacts with that topography, including effects of any modifications to the topography and/or surface hydrology (roads, bridges, mines, etc.), and precipitation quantity and rate (Fig. 1).
Figure 1 – Schematic of flood model inputs. A flood model includes (a) all components
of the natural landscape, (b) parameters that vary because of anthropogenic alterations,
and (c) robust precipitation data.
A working model that incorporates each of these elements can be assembled with i) a lidar-derived digital elevation model (DEM), ii) a materials raster parameterizing the DEM surface, and iii) quality hourly precipitation data. Details of each of these components and best practices for their assembly are discussed below.
Data Collection
To ensure proper model calibration and attribution of hydrologic properties, some preliminary data collection is necessary. In this example case, a field campaign was launched immediately following the July 2022 storm to obtain high-water marks to validate the model (Section 7), as well as field observations that would inform selection of landscape properties (Manningʼs n and infiltration rates, Sections 3.4, 5.2). Discharge data and ground precipitation data were also collected.
High-water marks to spatially validate model
If the aim of a model is to replicate a flood event, as much data as possible about the true event must be collected. Spatially, model performance is considered good if observed flood heights can be replicated within inherent uncertainties (GPS, elevation data resolution). For comparison, the study area should be traversed as soon as possible after a large storm, and quality indicators of maximum flood height should be identified and recorded. Care should be taken to record only reliable indicators of flood height. For example, debris in mature trees is reliable, but debris in young trees (which may have bent during the storm and picked up channel debris) can yield an overestimate of flood height. Koenig et al. (2016) detail full best practices for recording high-quality flood height indicators. In this example, we discuss a scenario where 21 high-water marks were collected following the July 2022 flood affecting eastern Kentucky (Fig. 2).
Figure 2 – Aerial extent of the catchment area examined in this study near the headwaters of the North Fork of the Kentucky River. In the July 2022 event, more than 350 mm of total precipitation fell over four days (inset map). Flood high-water marks (red circles) were measured across the catchment in September 2022. Black outline indicates the extent of the model in this study, brown polygons are mountaintop removal sites identified in aerial imagery, and the stream network is detailed in blue.
Discharge data to temporally validate model
Due to the abundance of reliable stream gages (typically installed by the United States Geological Survey, USGS) in major catchments, discharge data during extreme events is readily available. The utility of discharge data is two-fold: it is used to compare modeled discharge to confirm a temporal match and to assist in deriving infiltration rates. Once data is acquired for the event, a main storm peak should be identified. In the eastern Kentucky case study, the main storm surge is reflected in discharge data at the Whitesburg USGS stream gage (Fig. 2) between hours 57 and 72 (Fig. 3). Each hour of discharge data within the identified main storm surge should be averaged to match the temporal resolution of hourly satellite precipitation grids (discussed below) and summed to obtain a total discharge in m3s-1. This figure should then be multiplied by 3600 to convert to m3hr-1. This total discharge should be recorded for later use in Equation 1.
Figure 3 – Discharge data (example from the USGS gage at Whitesburg, KY; Fig. 2) is
used to identify a main storm surge, as well as to validate the temporal performance of
an established model.
Precipitation data
TUFLOW software can apply precipitation data to a model in multiple ways (for full options see the User Manual), but this example leverages a gridded rainfall approach and thus requires that a hydrology database be set up with one grid assigned to each model timestep (hourly timesteps are used in the eastern Kentucky model for 144 grids total). Satellite precipitation data are available at hourly intervals (Du, 2011) but should be calibrated to ground precipitation before modeling to account for (occasionally substantial) discrepancies between radar-forecasted precipitation and that measured in rain gages (Boushaki et al., 2009). Ground precipitation data for Kentucky are available at 15-minute intervals through Kentucky Mesonet (Mahmood et al., 2019); networks in other US states offer similar data resolution. Ground precipitation data should be gathered for the entire storm one wishes to model, and a comparable storm surge to that identified in the discharge data (Fig. 3) should be identified. One should expect a lag between the precipitation and discharge peaks. Point locations for each weather station from which ground precipitation data are downloaded should be added to an ArcGIS Pro project. After ground precipitation data are acquired and locations of weather stations have been added, proceed with the steps below for downloading, processing, and calibrating satellite precipitation data (if modeling in the United States) to ground values.
1. Obtain GRIB2 files from the storm for calibration and save unaltered in a folder. Hourly precipitation satellite grids can be batch downloaded from: NSF NCAR EOL data archive -- NCEP/EMC 4KM Gridded Data (GRIB) Stage IV Data (Du, 2011).
2. Once the GRIB2 files are downloaded and saved in a folder, open a Jupyter notebook in ArcGIS Pro and add four code blocks. Text in bold will need to be modified for each project/dataset. It is key that your ArcGIS Pro version is up to date, or the scripts will run without error but will not output to the map contents.
BLOCK 1:
# importing packages
import os
import pandas as pd
import arcpy
import tqdm
BLOCK 2:
# pathway to the folder that holds the multidirectional datasets (raw downloads from your link)
ras_direct = r"yourfilepath"
# pathway to the location that the .nc folder is being stored. Make an NC folder in your directory.
map_path = r"yourNCfolder"
# pathway to the project geodatabase
gdb_path = r"yourgeodatabase"
BLOCK 3:
for filename in os.listdir(ras_direct):
path = ras_direct + "\\" + filename
grb2_ras = arcpy.Raster(path, is_multidimensional = True)
test = grb2_ras.mdinfo
if test == None:
continue
else:
arcpy.env.addOutputsToMap = False
path = ras_direct + "\\" + filename
grb2_ras = arcpy.Raster(path, is_multidimensional = True)
# Process: Make Multidimensional Raster Layer (Make Multidimensional Raster Layer) (md)
raster_out = f'{filename[:-5]}_MDLayer'
arcpy.md.MakeMultidimensionalRasterLayer(in_multidimensional_raster = grb2_ras,
out_multidimensional_raster_layer = raster_out)
# Process: Raster to NetCDF (Raster to NetCDF) (md)
nc_path = map_path + '\\' + f'{filename[:-5]}_NetCD2.nc'
arcpy.md.RasterToNetCDF(in_raster = raster_out, out_netCDF_file = nc_path)
BLOCK 4:
for filename in os.listdir(map_path):
arcpy.env.addOutputsToMap = False
# Process: Make NetCDF Raster Layer (Make NetCDF Raster Layer) (md)
ras_name = f'{filename[0:-10]}_raster'
nc_path = map_path + '\\' + filename
nc_OP = arcpy.md.MakeNetCDFRasterLayer(in_netCDF_file = nc_path, variable= f'{filename[0:-10]}_MDLayer',
x_dimension="x", y_dimension="y", out_raster_layer = ras_name)
# Process: Saving NetCDF float as a raster in gdb.
arcpy.env.addOutputsToMap = True
new_name = filename[:-10].replace(".", "_")
ras_OP = gdb_path + "\\" + new_name
arcpy.management.CopyRaster(nc_OP, ras_OP)
3. Once rasters have added to the project (if they aren't added, make sure ArcGIS Pro version is current), run the Sample tool concurrently on all rasters at each location where there is a weather station with recorded ground precipitation. Export results as a spreadsheet using the Table to Excel tool and align satellite sampled values by time with the ground values. Make reasonable shifts (~4 hours in this case) to the data to obtain a good (r2 > 0.8) correlation on a crossplot (Fig. 4).
The resulting relationship between satellite and ground precipitation for the example model is:
Satellite, y = 0.8533 (ground, x) + 0.1166, or
Ground = (satellite grid value – 0.1166)/0.8533
Figure 4 – Cross plot of satellite and rain gage precipitation data used to correct hourly
satellite precipitation grids (after Boushaki et al., 2009). Gage precipitation was
collected from Kentucky Mesonet weather stations (n = 7) within the North Fork of the
Kentucky River watershed. After establishing a correction relationship, each satellite
grid for the modeled timespan was corrected with a raster calculation to uniformly
increase the precipitation to represent the volumes observed to have fallen in the model
space during the storm.
4. Now, each raster in the geodatabase needs to undergo this operation to correct satellite to ground values. Proceed with BLOCK 5.
BLOCK 5:
import arcpy
from arcpy import env
from arcpy.sa import *
import pandas as pd
arcpy.CheckOutExtension("Spatial")
print ("Checked out \"Spatial\" Extension")
BLOCK 6:
arcpy.env.workspace = r"yourgeodatabase"
arcpy.env.scratchWorkspace = r"yourscratchgeodatabase"
BLOCK 7:
raster_list = arcpy.ListRasters("*")
print (raster_list)
BLOCK 8:
for raster in raster_list:
name = str(raster)
if len(name) > 24:
continue
else:
ras = Raster(raster)
filename = 'cali_' + str(raster)[16:-4]
outraster = (ras-0.1166)/0.8533
outraster.save(r"yourscratchgeodatabase" + "\\" + filename)
5. Add all rasters from the scratch geodatabase to a new map. Check the values of a few rasters against their uncalibrated counterparts to make sure the calculation in BLOCK 8 was successful. Project the rasters into your working coordinate system (EPSG 3090 in this case) and export them as ASCII files into your hydrology database (an empty folder accessible by TUFLOW).
6. Type “Project Raster” in the Geoprocessing search bar. Right click the result and select “Batch." Input all calibrated rasters into the tool, specify coordinate system, and select “cubic convolution” for resampling technique. Run tool.
7. Type “Raster to ASCII” in the Geoprocessing search bar, right click, and select “Batch.”
8. Add all projected and calibrated grids to the input rasters, update the “environment” tab with your coordinate system, and select an empty hydrology folder as the output file location. Make sure the file extension is .asc instead of the default .txt.
9. Run tool.
If your output folder above has .asc files in it after the tool runs, add one to a map in ArcGIS Pro and confirm the coordinate system is correct and values are reasonable. If this is the case, TUFLOW will be able to read them correctly as a rainfall grid.
To ensure that the .asc files will be read properly by TUFLOW, refer to the User Manual for specific structuring of folders, file names, and hydrology spreadsheets.
Field observations to inform other model inputs
While collecting high-water mark data, general reconnaissance should be performed to inform modeling choices. In this example study, the goal was to assess if anthropogenesis exacerbated the July 2022 flood, so observations were made and conversations were had with landowners about what factors may have been at play. These included: channel narrowing, mountaintop removal mining, bridges becoming blocked with debris, and roads and buildings altering natural runoff and infiltration dynamics on the landscape. Other models will vary in terms of relevant field observations. In addition to any specific observations needed for a model, attributes such as roughness of surface materials, degree of vegetation and debris, soil type, etc. should be gathered so that parameters such as Manningʼs n values (Arcement and Schnieder, 1989) and infiltration rates can be reasonably assigned. If possible, infiltration rates can be directly measured (Bouwer, 1986) and applied to material or soil classes in the model. Alternatively, an expedient approach applying a catchment-averaged infiltration rate may be applied when field measurement of many infiltration rates is not possible. The details of this approach are discussed below.
Model Assembly
Base Topography
A functioning 2D flood model requires a representation of the landscape and a characterization determining how precipitation will interact with it. Both can be derived using quality lidar data. To begin this process, a point cloud should be accessed for the study area and a DEM should be generated by indexing last and single returns from the point cloud. In this example, a 1.5 m DEM was used for the materials raster derivation discussed below. For modeling efficiency, a lower resolution DEM (10 m) was then exported to be read by TUFLOW. Smaller models can of course operate at full DEM resolution.
Landscape Characterization
Materials raster
Once a base topography surface is established, a raster of values representing different material classes on the model surface is used to classify the landscape and assign infiltration rates and Manningʼs n values. In eastern Kentucky, material classes are based primarily on slope thresholds and vegetation coverage. To derive a materials raster for a model area, a slope map of equal resolution to the highest resolution lidar DEM (1.5 m in this case) should first be made. The raw lidar point cloud data that yielded the DEM should also be used to generate a classified vegetation layer. First and single returns from the point cloud represent the top of vegetation, and this surface minus the DEM yields a raster of vegetation height. Any unrealistic values for vegetation height (less than 0, greater than 45 m) should be eliminated (made to equal 0) with the Conditional tool, and the vegetation layer can then be classified by vegetation height: bare (0 m), grass (>0 to <0.5 m), shrub (≥0.5 to <5 m), and tree (≥5 m) (FGDC, 2008). Identifiers of 1, 2, 3, or 4 should be assigned using the Reclassify tool to represent bare ground, shrubs, grass, and trees, respectively. Based on reconnaissance of a study area, a reasonable threshold for the boundary between slopes and the floodplain should be determined. In eastern Kentucky, this is roughly 12°. Make the lidar-derived slope map binary according to this cutoff using the Conditional tool. In this instance, slopes greater than 12° were assigned a value of 20 and slopes less than 12° were assigned a value of 10. Finally, add the binary slope and classified vegetation layers together so that each possible sum serves as a material ID and represents a distinct material class (Table 1).
| A | B | C | |
| Material ID | Class | Manning's n Value | |
| 11 | Bare plain | 0.035 | |
| 12 | Grassy plain | 0.08 | |
| 13 | Shrubbed plain | 0.1 | |
| 14 | Treed plain | 0.2 | |
| 21 | Bare slope | 0.033 | |
| 22 | Grassy slope | 0.075 | |
| 23 | Shrubbed slope | 0.08 | |
| 24 | Treed slope | 0.18 | |
| 31 | Stream channel | 0.09 | |
| 32 | Buildings | 0.02 | |
| 33 | Roads | 0.02 | |
| 34 | Mines | 0.035 | |
| 35 | Ridges | 0.2 |
Table 1 – Material IDs and Manning’s n values for each class of the landscape surface.
Besides natural slope and vegetation, other materials known to affect surface rugosity per any field observations should be added to the materials raster at this point. For the eastern Kentucky example, stream channels, buildings and roads, and mine sites were incorporated into the materials raster with the following workflow, which may be adapted to suit other model areas and surface features: building footprints (ORNL and FEMA Geospatial Response Office, 2022) were downloaded, assigned a material ID, and converted from a polygon to a raster. Large mountaintop removal sites were identified using aerial imagery (Skytruth, 2024), mapped as GIS polygons, assigned a material ID of 34, and converted to a raster. A stream order raster was generated using a standard ArcGIS hydrology tools workflow on the high-resolution lidar DEM. The raster was converted to lines, with stream order retained as an attribute. To convert stream lines to polygons representing full channel widths, a series of buffer widths was tested until a best fit per stream order was established (Table 2). State and local road center lines were downloaded from KyGovMaps Open Data Portal and a buffering process similar to the one used for streams was applied (Fig. 5). Buffers used for each road category are shown in Table 3, and an example of the match to aerial imagery is shown in Fig. 6. Streams and roads were each assigned a material ID and converted to a raster. A series of conditionals in ArcGIS was used to add buildings, roads, streams, and mines to the base material raster for a final raster with 13 possible material IDs (Table 1). The materials raster was paired with a materials spreadsheet, where material IDs were linked to their respective Manningʼs n values.
| A | B | C | |
| Stream Order | Buffer (m) | Resulting channel width (m) | |
| 1 | 1 | 2 | |
| 2 | 2 | 4 | |
| 3 | 3 | 6 | |
| 4 | 5 | 10 | |
| 5 | 8 | 16 | |
| 6 | 13 | 26 | |
| 7 | 21 | 42 |
Table 2– Buffers and final channel widths applied to hydrologically-derived stream orders to reasonably match channel widths across the study area.
Figure 5 – Road center lines buffered variable amounts according to route category (Table 3). Each buffered road was assigned a material ID, converted to a raster, and added to the materials raster.
| A | B | |
| Route Type | Buffer Radius (m) | |
| private | 1.5 | |
| county | 2.25 | |
| state | 3 | |
| federal, city | 3.75 | |
| parkway | 4.5 | |
| interstate | 11 |
Table 3 – Buffer widths used to generate polygons from road center lines, then eventually a raster to add to the materials raster.
After burning in all extra classes (roads and buildings, streams, etc.), the final materials raster should look like this:
Figure 6 – Materials raster used to assign Manning’s n and infiltration values to the landscape. The raster is derived from slope and classified vegetation layers, except where streams, buildings, roads, and mine sites are present. These classes are each attributed separately.
The Raster to ASCII tool can be used to export this file into a folder to be read by TUFLOW alongside the model-resolution DEM.
Selection of Manning's n values
In this example, Manningʼs n values were determined with guidance from values established in literature (Jung et al., 2011) and were altered only when warranted by field observations (Table 1). Specifically, bare plains were assigned a value of 0.035 to represent a standard value for bare soil (Jung et al., 2011). Bare slopes were assigned a slightly lower value of 0.033 to account for the post-storm field observation of less debris on the slopes than in the floodplain. Grassy plains were given the value 0.08 to correspond to the modeled value for grasslands (Jung et al., 2011), and grassy slopes were slightly lowered to 0.075 due to less debris. Owing to the enhanced debris trapping observed post-storm when shrubs were present on floodplains, the Manningʼs n value was increased to 0.1 for this material class and lowered slightly to 0.08 for the equivalent slope class. Manningʼs n in forested areas is commonly on the order of 0.18 (Arcement and Schneider, 1989; Jung et al., 2011), which was maintained for forested (treed) slopes. For treed plains and ridges, Manningʼs n was increased to 0.2, as large trees trapped much anthropogenic and small vegetation debris. Though the value for water is typically ~0.03, channels were often observed to be completely clogged with debris and large vegetation, so a Manningʼs n of 0.09 was assigned to the stream network. The typical concrete/settlement value of 0.02 (Jung et al., 2011), was assigned to roads and buildings, and mines were given the bare soil value of 0.035.
Derivation of infiltration rates
Infiltration rates can be assigned to the same material classes discussed above. This can be done effectively using a catchment-averaged infiltration rate. This catchment averaged rate requires total precipitation from the calibrated satellite data (water in), total discharge (water out) determined above, the duration of the main storm surge determined from the peaks in discharge and precipitation (Fig. 3), and catchment area. Calibrated grids for the hours of the main storm surge should be extracted by mask using a polygon of the watershed of interest (the catchment area of this watershed polygon should also be determined). The Weighted Sum tool can be used with a weight of “1” and all rasters as inputs to obtain a raster of total precipitation for the main storm surge identified in hydrology data (Fig. 3). The Resample tool should be used on this sum to match the resolution of the DEM used for modeling (10 m in this case). The sum output raster should be divided by the resolution squared (10 × 10) to get an average precipitation total per pixel (m of precipitation per m of area). To convert this to a total precipitation volume, multiply the mean precipitation per pixel by the catchment area. This represents the total water input into the system. Using Equation 1 below and the values for total discharge, total precipitation, catchment area, and storm peak duration, calculate an infiltration rate.
Equation 1
To use this rate to assign rates to individual material classes, the materials raster (Fig. 7) can be separated by material class using a histogram and the total number of pixels within each class can be determined (Table 4).
Figure 7 – Process for deriving catchment-averaged infiltration. (a) Catchment area feeding USGS stream gage 03277300 at Whitesburg, Kentucky, underlain by total satellite precipitation grid for the storm peak (hours 57–72 beginning at 1 pm on July 25, 2022). (b) Total peak precipitation (hours 57–72) extracted within the watershed and resampled to match the 3 m DEM resolution. (c) Mean precipitation for the catchment on each 3 m cell. This mean multiplied by the catchment area yielded total precipitation volume, or water in, to consider in infiltration rate calculations.
Starting estimates for rates of each material class should be based upon soil type and degree of vegetation, then a weighted infiltration can be calculated by multiplying the number of pixels for each material class by the estimated infiltration. This weighted infiltration for each class represents the contribution of each class to the catchment-averaged infiltration rate. The sum of the weighed infiltrations divided by the total number of pixels in the materials raster should equal the catchment averaged rate determined above. Make minor and reasonable adjustments to the estimated rates of individual classes until the weighted infiltration divided by total pixels equals the catchment averaged rate. Once this is achieved, the materials raster can be exported as a “soils” raster, and the soils file required by TUFLOW can be populated.
| A | B | C | D | E | F | G | H | |
| Material ID | Class | Soil Type, Avg. Infiltration (mm hr-1) | # Pixels | % Area | Infiltration Rate Used for Model (mm hr-1) | Justification for Infiltration Modification | Weighted Infiltration | |
| 11 | Bare terrace | Allegheny Loam, 5.08 | 1198 | 0.20 | 3 | Decreased due to preceding storm, no vegetation | 3594 | |
| 12 | Grassy terrace | Allegheny Loam, 5.08 | 6344 | 1.05 | 3.5 | Decreased due to preceding storm, | 22204 | |
| 13 | Shrub terrace | Allegheny Loam, 5.08 | 10862 | 1.80 | 4 | Decreased due to preceding storm | 43448 | |
| 14 | Treed terrace | Allegheny Loam, 5.08 | 12031 | 1.99 | 5 | ~Average for soil, treated as baseline scenario | 60155 | |
| 21 | Bare slope | Shelocta-Highsplint-Gilpin complex, 7.62 | 10137 | 1.68 | 3 | Decreased due to preceding storm, no vegetation | 30411 | |
| 22 | Grassy slope | Shelocta-Highsplint-Gilpin complex, 7.62 | 20136 | 3.34 | 4 | Decreased due to preceding storm, little vegetation | 80544 | |
| 23 | Shrub slope | Shelocta-Highsplint-Gilpin complex, 7.62 | 80325 | 13.31 | 7.2 | Decreased due to preceding storm, no trees | 578340 | |
| 24 | Treed slope | Shelocta-Highsplint-Gilpin complex, 7.62 | 341146 | 56.52 | 8.6 | ~Average for soil, treated as baseline scenario | 2933856 | |
| 31 | Stream channel | Combs silt loam, 7.62 | 13290 | 2.20 | 0 | Assumed to be saturated | 0 | |
| 32 | Buildings | n/a | 2129 | 0.35 | 0 | Impermeable (Lederle et al., 2020) | 0 | |
| 33 | Roads | n/a | 10191 | 1.69 | 0 | Impermeable (Lederle et al., 2020) | 0 | |
| 34 | Mines | n/a | 69182 | 11.46 | 2.5 | Negley et al., 2006 | 172955 | |
| 35 | Ridges | Shelocta-Highsplint-Gilpin complex, 7.62 | 26639 | 4.41 | 10 | Increased from baseline scenario due to low gradient | 266390 | |
| 603610 | 4191897 |
Table 4 – Derivation of infiltration rates for each material class of the model. Infiltration rates were adjusted from baseline soil estimates until the sum of the weighted infiltrations divided by the total pixels in the model equaled a catchment-averaged rate of 6.94 mm hr-1.
Shapefiles to Layer atop the Model Landscape (Bridges, Culverts, Dams, etc.)
Clogged bridges may substantially alter flow through a catchment and should be considered in a model. Similarly, culverts may allow greater flow through embankments than suggested by elevation data, and should thus be incorporated as well. Both can be added as separate shapefiles read by TUFLOW. Bridges are modeled in 2D as having four layers: space below deck (L1), the deck (L2), rail (L3), and clearance above the rail (L4) (Fig. 8a). Although exact bridge dimensions vary, for simplicity, a deck thickness of 1 m and a rail height of 1 m were assigned to all bridges in the example model in accordance with field observations.
Figure 8 – Schematic cross-section view of bridge scenarios considered in model. (a)
Normal scenario, in which the channel is clear of debris and water flows unobstructed
beneath the bridge deck. (b) Flood scenario, in which debris enters narrowed channels
and blocks bridges. Bridges less than 5 m high were assigned blockages of 100% for
layers L1 (clearance), L2 (deck), and L3 (rail) based on field observations.
Height measurements for eastern Kentucky bridges are not included in state bridge geospatial data (KYGovMaps Open Data Portal, 2025). Fortunately, heights can be estimated manually from visible ramparts in the lidar DEM adjacent to each bridge in the feature class. Heights for eastern Kentucky bridges averaged 3.2 ± 1.8 m (Fig. 9a); based upon this value and field observations of washout and blockages, all bridges less than 5 m high were assumed to be fully blocked by debris. Thus, a modeled bridge with 5 m height plus the 1 m bridge deck and 1 m high rail yielded a 7 m zone of 100% blockage, which represents a worst-case scenario in which all “short” bridges (n = 98) across the area would be fully blocked (Fig. 8b). Taller bridges (>5 m high) were blocked under their decks by only a percentage of their height (5 m of below deck space, plus the 1 m deck and 1 m rail). For bridges over 7 m high (n = 22), no rail blockage (Fig. 8, L3) was applied. Specific degrees of blockage and distributions of heights will vary from these results in different model areas with different bridge types, but the same expedient approach for height estimation can be applied. One large culvert was identified in aerial imagery in the eastern Kentucky study area. Due to the options available for applying culverts to the model in TUFLOW, the culvertʼs parabolic shape was modeled as circular by matching the cross-sectional area (Fig. 9b). It was assumed to be only 10% blocked because it was immediately downstream of a low-lying bridge, which shielded it from most debris. A Manning's n value of 0.045 was assigned to the culvert (half that of the underlying channel, n = 0.09) because, although the channel bottom is exposed, the pipe above it is metal.
Figure 9 – (a) distribution of bridge heights measured on DEM across the study area.
Bridges less than 5 m tall were treated as dams in the pilot model (Fig. 8b). (b) Google
Earth imagery of the one large culvert in the study area, which was modeled as circular
and blocked 10%.
Model Validation
Models can be validated temporally and spatially. To assess model performance temporally, time series data should be loaded from the model by setting up plot output points at gage locations. Peak discharge extracted from the model should be compared to the true peak discharge from any gages in the study area. For the eastern Kentucky model, peak discharge was 282 m3s-1 at hour 72.5. If compared to the Whitesburg gage peak of 278 m3s-1 at hour 70, temporal accuracy is 98.6%.
For spatial validation, any high-water marks (as discussed above) should be brought into GIS space, and the uncertainties associated with both the GPS and the vertical uncertainty of the lidar data should be known. Model outputs include HMax, a map of maximum water height. A buffer with radius equal to GPS uncertainty can be applied to each high-water mark point location. The minimum and maximum HMax values output by the model can be extracted from within this buffered zone at each point. The vertical uncertainty inherent in the elevation data accounts for one sigma vertical uncertainty. In the example model, DEM uncertainty is 30 cm. A model may be considered accurate if the measured high-water mark falls within one sigma of the extracted HMax model value. In the eastern Kentucky model, 19 of the 21 points matched high-water marks within these criteria for a spatial accuracy of 90.5% (Fig. 10b). The two modeled points that did not match were investigated and found to be within one 10 m cell of being a match. To revise the model in the event of a temporal or spatial mismatch, individual Manning's n values or infiltration rates may be adjusted by reasonable degrees. Specifically, increases in Manning's n values across all material classes may delay peak discharge, while decreases will accelerate and increase it (Swallom et al., 2025).
Figure 10 – Results and validation of the model. (a) HMax map output displaying time-independent spatial patterns of flood height. Time series data extracted from the model at plot output points and the observed data from a USGS gage location (black dashed curve) show that the model yields a temporal discharge response that compares favorably with observed discharge peak (accuracy = 98.6%). (b) Zoomed in view of the HMax surface atop topography and buildings near the USGS gage and an example of the process for spatially validating the default model. Yellow circle indicates a high-water mark of 2.41 m. (c) 1:1 plot between modeled flood heights from the default model and flood heights recorded in the field following the July 2022 storm. The green line is a reference line with slope = 1, while the solid black line is the linear regression line between observed high water marks and modeled flood heights. The black dashed lines denote uncertainties of one sigma; the reference line falls within this range, indicating that the model performs well to a maximum flood height of 17.7 m. Beyond this height, the model will underestimate Hmax by more than one sigma.
Protocol references
Arcement, G.J., Jr., and Schneider, V.R., 1989, Guide for selecting Manningʼs roughness coefficients for natural channels and flood plains: U.S. Geological Survey, Water Supply Paper 2339.
Boushaki, F.I., Hsu, K.-L., Sorooshian, S., Park, G.-H., Mahani, S., and Shi, W., 2009, Bias adjustment of satellite precipitation estimation using ground-based measurement: A case study evaluation over the southwestern United States: Journal of Hydrometeorology, v. 10, no. 5, p. 1231-1242.
Bouwer, H., 1986, Intake rate: Cylinder infiltrometer, in Klute, A., ed., Methods of soil analysis: Part 1 [2nd ed.]: Madison, WI, Soil Science Society of America, p. 825–844, https://doi.org/10.2136/sssabookser5.1.2ed.c32.
Du, J., 2011, NCEP/EMC 4KM gridded data (GRIB) stage IV Data: UCAR, NCAR Earth Observing Laboratory, https://doi.org/10.5065/D6PG1QDD [accessed 03/15/2024].
Federal Geographic Data Committee Vegetation Subcommittee (FGDC), 2008, National Vegetation Classification Standard, version 2: Reston, VA, Federal Geographic Data Committee, FGDC-STD-005-2008, 119 p.
Jung, I.-K., Park, J.-Y., Park, G.-A., Lee, M.-S., and Kim, S.-J., 2011, A grid-based rainfall-runoff model for flood simulation including paddy fields: Paddy and Water Environment, v. 9, p. 275-290, https://doi.org/10.1007/s10333-010-0232-4.
Koenig, T.A., Bruce, J.L., OʼConnor, J., McGee, B.D., Holmes, R.R., Jr., Hollins, R., Forbes, B.T., Kohn, M.S., Schellekens, M., Martin, Z.W., and Peppler, M.C., 2016, Identifying and preserving high-water mark data: U.S. Geological Survey, Techniques and Methods 3-A24, https://doi.org/10.3133/tm3A24.
KYGovMaps Open Data Portal, 2025, Commonwealth of Kentucky, KyGovMaps Open Data Portal [accessed 01/22/2025].
Lederle, R., Shepard, T., and De La Vega Meza, V., 2020, Comparison of methods for measuring infiltration rate of pervious concrete: Construction and Building Materials, v. 244, article 118339, https://doi.org/10.1016/j.conbuildmat.2020.118339.
Mahmood, R., Schargorodski, M., Foster, S., and Quillligan, A., 2019, A technical overview of the Kentucky Mesonet: Journal of Atmospheric and Oceanic Technology, v. 36, no. 9, p. 1753-1771, https://doi.org/10.1175/JTECH-D-18-0198.1
Negley, T.L., and Eshleman, K.N., 2006, Comparison of stormflow responses of surface-mined and forested
watersheds in the Appalachian Mountains, USA: Hydrological Processes, v. 20, no. 16, p. 3467-3483,
Oak Ridge National Laboratory (ORNL) and FEMA Geospatial Response Office, 2022, USA structures: USA_Structures_View (FeatureServer) [accessed 2024].
SkyTruth. 2024. https://skytruth.org.
Swallom, M.L., Dortch, J.M., and Thigpen, J.R., 2025, Assessing the effects of anthropogenesis on Appalachian flood severity: An eastern Kentucky case study: Geomorphology, v. 482, article 109818, https://doi.org/10.1016/j.geomorph.2025.109818.
Acknowledgements
Parts of this work were supported by the National Science Foundation (EAR 2242120), which enabled us to expediently collect and leverage perishable data following the July 2022 flood in eastern Kentucky. Computational infrastructure in the Surface Processes and Modeling Laboratory that was used for this work was also supported, in part, by the National Science Foundation EPSCoR CLIMBS Project (OIA 2344533).