Development of GIS-based Python scripts to calculate a water surface profile on a landscape for wetlands decision-making

The data, code and tools developed in this research are available via https://github.com/TrauthK/Wetlands/tree/master/Hydraulics/Spatial%20backwater%20Curves. The tools were developed in Python 2.7 with ArcGIS 10.3.

As the name suggests, wetlands are locations that are associated with wetness. This wetness causes a transformation of the soil because of the absence of free oxygen, which in turn supports plants that thrive in this specific anaerobic environment. However, this wetness does not need to occur continuously, and may only be required for a short period of time, as during the growing season for the vegetation of interest.

Wetlands provide multiple beneficial services for the environment and human society. Different kinds of plants and animals are supported by wetlands, and the dynamic relationship among these creatures is very complex, which makes wetlands one of the most productive ecosystems on Earth, providing crucial habitats for many plants and animals, including birds, amphibians and reptiles (Snodgrass et al. 2000; Gibbons et al. 2006). Conversely, human activities, such as construction, pollutant discharges, deforestation and farming, can significantly impact the quality of water resources. For these reasons, interest has grown in the concept of mitigation wetlands – undertaking an effort to create or recreate wetland conditions. Informed decision-making regarding mitigation wetlands requires a thorough understanding of the water balance over time.

For some time, researchers have been interested in the connectivity between wetlands and how wetlands connect to downstream locations (Amezaga et al. 2002; Roe et al. 2009; Lang et al. 2012). Many geographically isolated wetlands can connect with other water bodies during wet seasons. During storms, wetlands receive water from their own contributing areas. After wetlands fill up, surface water may start to flow to neighboring wetlands at lower elevations or to downstream streams and rivers. Quantification of the processes governing hydrologic connectivity between wetlands is thus important for the scientific and decision-making communities (Golden et al. 2014). These wetlands that are connected, even if only periodically, can form a habitat that helps to support plant and animal populations.

The flow in a channel with a regular cross section (rectangular or trapezoidal) has been studied for many years (Manning 1891). Modeling the flow over a landscape is a difficult undertaking and introduces additional complexity. As water drains into a depression from the surrounding contributing area, it will exist at a certain depth that will vary based on the size of the contributing area and the type of storm generating the runoff. As water flows between depressions, there exists a water surface profile that connects the depth of water in one depression with the depth of water in the connecting upstream or downstream depression. A water surface profile depicts the changing flow depth in the longitudinal direction of flow. On the hydraulically mild slopes of interest for wetlands mitigation assessments, the profiles of interest are backwater (M1) and drawdown (M2) curves.

Water surface profiles play an important role in engineering practice. Infrastructure, such as channels, sewers or culverts, should be designed with the consideration of gradually varied flow. Different methods to calculate the water surface profiles have been developed, including the Standard Step method, the Direct Step method and the Modified Euler's method (Chow 1959; Sturm 2010). These methods are based on the energy, momentum and continuity equations.

The Standard Step computational procedure is a trial and error process. Based on the given discharge and the boundary and channel conditions, the method assumes an initial water surface elevation for a target cross section. It then recalculates a new water surface by calculating the headloss as a result of friction between the two cross sections. Iterative calculations with different assumed initial water surface elevations continue until the difference in the total heads is within an acceptable range.

Chow (1959) indicated that profile calculations should be carried out from downstream to upstream if the flow is subcritical and from upstream to downstream if the flow is supercritical. Performing calculations in the wrong direction could cause the water profile to diverge from the correct one. The Standard Step method can be used to perform calculations in either direction depending upon the flow conditions.

HEC-GeoRAS (U.S. Army Corps of Engineers 2011) has been developed to perform hydraulic calculations utilizing landscape data input through ArcGIS to produce water surface profiles, among other outputs. This tool is able to derive cross sections at specified locations from a digital elevation model (DEM), once a user identifies the centerline of the flow path, the locations of cross-sections of interest and the locations of the channel banks. Manning's n values must be entered manually, as there is no capability to translate a land use/land cover layer to specific roughness values for a channel and the overbanks. While HEC-GeoRAS is strictly hydraulic in its capabilities, it is able to perform calculations based on the hydrologic output from HEC-GeoHMS (U.S. Army Corps of Engineers 2013). HEC-GeoHMS does incorporate various loss methods, such as the Green and Ampt infiltration method (Green & Ampt 1911), but the input parameters are held constant over time and do not represent the actual fluctuations that will occur. For example, the usage of the Green and Ampt Loss method requires the specification of a moisture deficit (the difference between the saturated moisture content and the initial moisture content) that is used for the entire period of calculation. In reality for a site, the soil moisture deficit will change over time with precipitation, infiltration and evapotranspiration (ET). ET itself is not incorporated into HEC-GeoHMS at all. Additionally, a set value for initial water loss must be entered, rather than being calculated as water that is infiltrated until the point of surface saturation.

Arc Hydro (Maidment & Morehouse 2002) provides users with a number of ArcGIS-based tools for landscape-based water resources investigations, which are similar to the spatial analysis techniques embedded in ArcGIS, including watershed delineation and generating flow paths. These capabilities include the ability to provide input data for HEC-GeoRAS. Also included within the category of hydrology and hydraulics is a Green and Ampt tool for the assessment of infiltration impacts. The tool requires the input of Green and Ampt parameters, which include the previously discussed moisture deficit and initial loss parameters that are treated as constants.

When modeling over a landscape, the capability to analyze a DEM and extract useful information for hydrologic and hydraulic purposes is quite mature. Various methods offer the ability to delineate a watershed boundary (Jenson & Domingue 1988) and generate the flow direction and flow path (Quinn et al. 1991). TauDEM (Tarboton 1997) is a set of tools capable of deriving landscape characteristics necessary for hydrologic analyses. Among its capabilities is a method to determine flow direction by proportioning flow between the two lowest downstream cells, recognizing that water flow is not limited to a single receiving cell. Being hydrologic in nature, the TauDEM suite does not incorporate infiltration or ET capabilities, nor does it incorporate the hydraulics necessary for the calculation of a water surface profile. In evaluating specific tools to utilize in a landscape analysis, the Height Above the Nearest Drainage (HAND) (Nobre et al. 2011) is a tool that can classify bare ground cells above a stream into different categories (e.g., 0–0.2, 0.2–0.5 and 0.5–1.0 m) above the stream bed or the water surface elevation). Such a classification could be used in developing inundation maps. However, HAND does not have the capability to work with the varying depths of flow associated with a backwater curve that may exist in slow-moving flood flows, nor is it integrated into either HEC-GeoRAS or Arc Hydro.

Wetland mitigation efforts may want to focus on natural depressions that, with minimal effort, could be converted into locations with sufficient water for sufficient periods of time to actually function as wetlands. When investigating the potential for essentially creating new wetlands, the questions that must be answered include how much water will exist at a site, where it will exist and for how long. Answering these questions requires a water balance evaluation that includes precipitation as well as infiltration and ET. Infiltration and ET are processes that are intertwined based on soil moisture content and change throughout the year and with precipitation. Modeling of moisture content must be conducted with a very small time step (e.g., 1 h) and becomes very computationally intensive over the long time periods over which potential mitigation wetland sites must be evaluated (a year or more). Both infiltration and ET are impacted by the distribution of water over a landscape over time, where the distribution is impacted by both hydrologic and hydraulic processes. Thus, a thorough water balance assessment requires an integrated process incorporating not just hydrology and hydraulics but also the land cover and soil processes of infiltration and ET. Land managers engaged in decision-making regarding mitigation wetlands may need to investigate multiple sites, and thus need a methodology that can be implemented relatively quickly and easily.

The objectives of this study are to investigate and hydraulically characterize the flow that could exist between two depressions. The characterization is in terms of a water surface profile that could ultimately be used to determine the aerial extent of water flowing over a landscape. It is not necessary for a location to be continuously covered with water for wetland conditions to develop and be maintained – for plants, it is only necessary for surface inundation during the growing season. As can be seen with HEC-RAS output, the water surface profile in the context of a cross section shows the aerial extent of water on the ground surface. A water surface profile, in combination with an assessment of infiltration, could indicate locations with potential surface saturation.

This study is part of a larger effort to assess the spatial and temporal occurrence of water on a landscape for decision-making with respect to the creation of mitigation wetlands. The overall effort requires the incorporation of precipitation and the hydrology and hydraulics of water flowing over a landscape with the subsurface processes of infiltration and ET. As a landscape-based assessment, a geographic information system (GIS) analysis is required on a scale that is commensurate with the scale of depressions observed on the landscape (<1 ha). Because of the ultimate usage of the modeling results, this determination of water surface profiles must be undertaken within an overall system that incorporates the flow of water over a surface with the possibility of infiltration and ET.

The overall research intent is to combine hydrology and hydraulics on a landscape in a GIS format. It is directed at decision-making for mitigation wetlands, and thus it is necessary to track parameters within a water budget (namely infiltration and ET) over time. This set of capabilities is not currently available in a single suite of software. Thus, it was necessary to embark on a research effort to integrate all these capabilities. Because of the growing utilization of Python scripts within the GIS platform, it was decided to create a Python-based platform that would be able to accommodate the individual scripts addressing overland flow and water surface profiles as well as infiltration and evaporation.

The overall project, of which this study is a part, requires the assessment of depressions on the landscape to identify locations that with minimal earthwork changes may have the potential to be converted to functioning wetlands. Such an analysis requires long-term continuous simulation (a year or more) of a water balance that incorporates hydrology and the hydraulics of water moving over a landscape but also the processes of infiltration and ET. Because existing tools have not been identified that meet all of these criteria, the overall plan became to develop a new Python-based tool. The study described in this paper regarding water surface profiles is only one part of the overall effort but must ultimately be able to be integrated into the larger model.

The Modified Euler's method is employed to calculate the water profile caused by backwater or drawdown effects between depressions on a mild slope. The Modified Euler's method was selected because it can directly calculate water surface profiles on the hydraulically mild slopes of interest here. In contrast to the iterative calculation process of the Standard Step method, the Modified Euler's method is more computationally efficient. This computational efficiency is essential given the fact that the overall model will perform a water balance over the period of at least a year, all the while calculating infiltration and ET on an hourly basis.

A work flow is proposed to simulate water surface profile development over a landscape (rather than in a defined channel) using a combination of readily available GIS tools and two Python scripts written for this study. Python is used because of its wide and growing usage in geospatial modeling. In this work flow, GIS tools extract detailed landscape information directly from a high-resolution DEM, which is very helpful, especially for the vast majority of depressions (potential wetland mitigation locations) where there are no field measurements. The two scripts written for this study are then used to derive hydraulic landscape information (i.e., cross-sectional area, wetted perimeter and hydraulic radius) and calculate the water surface profile for the landscape.

The consideration of locations for potential wetland mitigation sites may raise the issue of whether a DEM containing bathymetric information may be necessary for the analysis. Bathymetric DEMs are limited in availability and so could present a difficulty (U.S. Geological Survey 2019). Because the analysis is conducted on depressions, rather than on existing wetlands where there may be surface water, the scarcity of bathymetric DEMs is not an issue in this instance.

The computational process is shown in Figure 1. Once a user of the methodology identifies the general location of interest for analysis, a high-resolution DEM for the area is identified. Within the general area of interest, the user must identify two depressions for specific analysis. Tools readily available in ArcGIS are then used to delineate the contributing areas for each depression and identify the weir cross-section location (pour point) of the upper depression (which controls the outflow from the upstream depression). The starting point is normally recognized as the pour point of the upper depression which is the lowest point of the weir. Correspondingly, the lowest point of the downstream depression is considered as the ending point. The least-cost line (the steepest route) between these two points generated from the DEM is treated as the central line of the flow pathway.

The correct implementation of the Modified Euler's method using the scripts is demonstrated in two different ways. The first demonstration is conducted by using the scripts to extract the hydraulic input parameters for a trapezoidal channel from an artificial landscape and then calculating the backwater curve for an example provided by Sturm (2010). In this fashion, it is easy to ascertain that the Python script is carrying out the calculations properly. The second demonstration is performed for two different precipitation events on two connected depressions in the vicinity of Pershing State Park in Linn County, Missouri. This location is selected because of the historic presence of wetlands and current interest in wetland mitigation.

To verify the accuracy of the Python script, an artificial trapezoidal channel is created in ArcGIS using elevation data (Figure 2). It has the same channel geometry as an example provided by Sturm (2010). This trapezoidal channel has a channel slope of 0.001 m/m, a side slope ratio of 2:1 and a bottom width of 8 m. The Manning's n of the channel is 0.025. The Python script detects the DEM of the trapezoidal channel and extracts all geometric information directly from it. The script then calculates the hydraulic parameters (i.e., cross-sectional area, wetted perimeter and hydraulic radius), as they vary by water depth (Table 1).

Table 2 shows the results of the backwater calculations from Sturm (2010) and the Python script. Column 1 shows the distance from the starting point 0 to the upstream ending location 20 m away. At the beginning of the calculation, the x intervals are very small in order to obtain precise results where there may be significant curvature of the flow, as with an M2 water surface profile. Column 2 lists the water depths as given by Sturm (2010). Column 3 depicts the water depths calculated using a spreadsheet created for this purpose, while Column 4 displays the depths calculated using the Python script. A spreadsheet is used as an interim step in the coding to establish the calculational strategy. The very minor differences between these two processes indicate that the Modified Euler's method can be implemented using this Python script.

The procedure of calculating the water surface profile at an actual landscape location is demonstrated in this section.

A user of the methodology must identify a location of interest for potential mitigation wetlands. Pershing State Park in Linn County, Missouri currently contains wetlands. Because of the proximity to wetland locations and the potential for habitat connectivity with mitigation wetlands, this demonstration is undertaken in a nearby location. A 1-m resolution DEM is created based on the LiDAR data from the Missouri Spatial Data Information Service (MSDIS 2019). The scale of the demonstration location was specifically chosen to be representative of the scale of the areas over which this process would be implemented. Natural wetlands may be small, so this methodology is intended to be used at a very small scale (tens of hectares) to identify those locations that could possibly support small wetland mitigation sites. In addition, the simplicity of two depressions was specifically selected to mimic the assessments decision-makers may be conducting.

A user of the methodology must identify specific depressions for analysis. The study area includes two depressions connected by a landscape with a mild natural slope (Figure 3). In the DEM, lighter areas represent higher elevations, while darker areas represent lower elevations. During a rainfall event, these two depressions would start to accumulate surface runoff from their respective contributing areas. Once the runoff would fill up the upper depression, excess runoff would flow out of it and contribute to the lower depression. The lower depression, with a larger contributing area, may have a deeper water depth which could cause drawdown or backwater effects.

The contributing area for each subbasin can help to explain the flow accumulation. The first step is using the embedded spatial analysis tools in ArcGIS to generate boundaries for each depression (Figure 3(a)). In this example, depressions are identified from the DEM.

With sufficient rain, surface water will flow out of the upper depression and contribute to the lower depression. The starting point for flow (Figure 3(a)) is located at the lowest point along the boundary of the upper depression (i.e., the pour point of the upstream watershed). The ending point is the lowest point in the lower depression (i.e., the sink).

Surface runoff always flows along the steepest slope over the landscape. The distance tools in ArcGIS can create a cost map based on the slope of the landscape and then generate the least-cost path (Figure 3(b)) which can be recognized as the centerline of the flow. For this study site, the elevation of the starting point is 228.65 m and the elevation of the ending point is 227.52 m. The total length of the flow is 229 m, and the slope of the overall flow pathway is 0.0049 m/m.

In general, separate reaches are identified where the hydraulic characteristics of slope, land cover or geometry are consistent. For purposes of representing the different landscape characteristics in this demonstration, the flow path is divided into three reaches of equal length (Figure 3(c)). Each reach has different slope and cross-section parameters which are automatically extracted from the DEM using the Python script. First, the script creates perpendicular lines from the middle point of each reach (Figure 3(d)). The elevation data along each perpendicular line are extracted from the DEM. Second, the cross section for each reach is built from these elevation data. The width of each cross section is set at 100 m and can be modified by the user as necessary. In standard hydraulic convention, reach lengths are numbered starting from the downstream-most location.

Table 3 shows the elevation data of the starting and ending points for each reach, along with the length and slope. Figure 4 shows how the geometry of different cross sections can vary significantly, having a great impact on the hydraulic calculations. Compared with Reaches 1 and 3, the geometry of Reach 2 has steeper sides (impacting the cross-sectional area of flow), resulting in an increasing depth when water flows from Reach 3 to Reach 2 and a decreasing depth when the water flows from Reach 2 to Reach 1.

The water surface profile for each reach is calculated from downstream to upstream in a continuous process. Starting with a given water depth, the script calculates the water depth at the end of Reach 1 and uses that depth as the starting water depth for Reach 2. Here, the profiles are calculated for two different flow circumstances (i.e., from a more intense Storm 1 generating greater flow and from a less intense Storm 2 generating less flow).

Figure 5 shows the results from Storm 1 (with a larger discharge and a deeper downstream depth) as generated by the Python script. In this test, the discharge of the flow is 200 m³/s with an initial depth of 1.3 m at the bottom of the lower depression. The roughness coefficient for all three reaches is set to 0.025. The thicker solid line indicates the channel bottom, and the right vertical axis is elevation. The dotted line is the calculated uniform depth. The uniform depth changes between the different reaches because the landscape characteristics are different between the reaches. The thinner solid line represents the water surface profile, which always varies in terms of approaching uniform depth asymptotically. The 1.3 m depth in the downstream depression is greater than the calculated uniform depth of 0.90 m. Therefore, the water surface must decrease with an M1 curve in the process of approaching the equilibrium condition of uniform depth. Even as the uniform depth changes, as from Reach 1 to Reach 2, the water surface is shown to be continuous, as would actually occur on the landscape. The transition from Reach 2 to Reach 3 shows the continuous water surface, with the shape adjusting from the steeper sides of Reach 2 to the flatter sides of Reach 3. For most of Reach 3, the plot shows the water surface profile coinciding with the uniform depth line, indicating that the water surface profile is very close to the uniform depth in its asymptotic approach, as would be expected.

For Storm 2, the discharge is 5.0 m³/s with a given depth 0.3 m at the bottom of the lower depression. The roughness coefficient for Reach 2 is set to 0.05, and Reach 1 and 3 are unchanged at 0.025. In this scenario, the flow depth in Reach 1 decreases gradually following an M1 curve (Figure 6) because the water depth in the downstream wetland is larger than the uniform depth of 0.20 m in the region connecting the two wetlands. Unlike the plot from Storm 1, there is an M2 curve in the transition between Reach 1 and 2, as the water surface now needs to increase to approach uniform depth. This condition exists because a larger roughness coefficient is assumed for Reach 2 resulting in a deeper uniform depth of flow, even without changing any of the other parameters. Once in Reach 3, water depth again follows an M1 curve to decrease the depth of flow to approach the lower uniform depth.

This paper introduces an efficient methodology for using the Modified Euler's method to calculate water surface profiles over a natural landscape on the mild slopes connecting wetlands. Python scripts have been created to perform the calculations, including the extraction of landscape data directly from a DEM, eliminating the need for a field survey at every potential depression of interest. The efficiency of the methodology allows for widespread scoping calculations to assess locations for more detailed investigations for wetland protection or mitigation. Its simplicity and automation can help users with little hydraulics background to calculate profiles for decision-making.

Multiple tests have been performed to verify the accuracy of this methodology. The Python scripts created for this project automatically extract the necessary data and compute all necessary hydraulic parameters. The first test recreated a drawdown curve demonstration in a textbook. This was accomplished by generating an artificial landscape matching the trapezoidal channel of the example. The Python scripts extracted the necessary landscape parameters from the DEM, calculated the appropriate hydraulic parameters and developed the drawdown curve, exactly matching the published results. The second test was performed on a natural landscape with a discharge that would be associated with an intense storm. This set of conditions produced a sequence of backwater (M1) curves expected as landscape characteristics change, causing the water surface profile to adjust accordingly in moving from one reach to another. The third test was conducted on the same natural landscape (while changing one surface roughness value) but with a less intense storm. The lower flows associated with these more common, less intense storms, are of particular interest for continued wetland functioning. With the lower discharge, the script demonstrated the ability to calculate a water surface profile that incorporates both backwater (M1) and drawdown (M2) curves. For the development of the methodology, in order to focus on the surface processes, infiltration was assumed to be negligible.

This effort addresses the challenge of applying open channel principles to quantify the flow characteristics (i.e., a water surface profile) over a landscape in contrast to the current focus on flow in defined channels. It is able to accommodate the inherent variability in the landscape elevations and land covers. It is now possible to make assessments that could provide a first-step screening tool for wetland scientists to assist in the selection of locations that may warrant further investigation with respect to the preservation, restoration or mitigation of wetlands.

Now that the water surface profile methodology has been developed and demonstrated, it will be appropriate to extend the depiction of the water surface profile on the landscape to include the extent of the land surface that would be covered by water. Such water on the surface would be available for infiltration to support wetland function.

This work was supported by US Environmental Protection Agency, Region 7 Wetland Program Development Grant CD97748001.