Tuolumne Meadows floodplain modeling using HEC-GeoRAS (Part 2)

After a long break, today we will see how to import the preprocessing data generated in part 1 of this tutorial into the HEC-RAS hydraulic model.

HEC-RAS is a program that models the hydraulics of water flow through natural rivers and other channels. It is one-dimensional, meaning that there is no direct modeling of the hydraulic effect of cross section shape changes, bends, and other two- and three-dimensional aspects of flow.

Launch the HEC-RAS program (Start > Programs > HEC > HEC-RAS > HEC-RAS 4.1.0). The latest version of HEC-RAS can be downloaded for free from the Hydraulic Engineering Center of the U.S. Army Corps of Engineers website: http://www.hec.usace.army.mil/software/hec-ras/hecras-download.html

Save the new project. Go to File > Save Project As… Save it as “TuolumneMeadows.prj”

To import the GIS data into HEC-RAS, first go to the geometric data editor (Edit > Geometric Data…)

Next, click on File > Import Geometry Data > GIS Format. Browse to “GIS2RAS.RASImport.sdf” file previously created in HEC-GeoRAS, and click OK. An import wizard will appear. In the first window (Intro), confirm SI (metric) units.

In the second window (River Reach Stream Lines), make sure all import stream lines boxes are checked, and click Next. Finally, in the last window (Cross Sections and IB Nodes), make sure all Import Data boxes are checked for cross-sections and click OK.

Click on Finished-Import Data. Data will be imported into the HEC-RAS geometry editor.

Save the geometry file as “TuolumneGeometry” (File > Save Geometry Data). Note that the geometry data has a .g01 extension.

To check the quality of the imported data, we can use the graphical cross-section editor (Tools > Graphical Cross-section Edit) in the geometric data editor.

We can use the editor to move bank stations, change the distribution of Manning’s n, add or delete ground points, etc.

The next step is to edit the structures data (bridges). Click on “Bridge/Culvert” editing button. In our case, we have two bridges at river stations RS=5360.882 and RS=4123.291.

Pacific Crest trail bridge (RS=4123.291)

We see that, in both cases, the deck elevation is not accurate. So we need to edit this information. Click on Deck/Roadway editor, and delete all information. Enter the data for the deck hight and low chords height as shown in the next picture:

Click on the “Copy US to DS” button to copy the information from upstream to downstream, and click OK. The result should look something like this:

Now, we will enter the pier information. The first bridge (RS=5360.882) has three piers and the second (RS=4123.291) has only one. Click on the “Bridge Design” button and enter the information as shown below:

Close the bridge/culvert editor. Save the geometry data and close the geometric editor.

The next step is entering Flow Data and Boundary Conditions. In HEC RAS every flow to be simulated is called a profile. In our case, we are going to define one profile that represents the peak discharge at the most upstream location of the Tuolumne River for a 100 year return period.

100-year flows were assessed from a frequency analysis of the nearest gauging station data (1127479 Tuolumne River at Grand Canyon of Tuolumne), with the maximum discharge being 200 cubic meters per second (7060 cubic feet per second). Enter 1 for number of profiles, and click on “Apply Data”.

We now define flow conditions upstream for steady flow (Edit > Steady Flow Data). To define downstream conditions, click on “Reach Boundary Conditions”, select “Downstream” for Toulumne River, click on Normal Depth, and enter 0.001.

The basic computational procedure of HEC-RAS for steady flow is based on the solution of the one-dimensional energy equation. Energy losses are evaluated by friction and contraction – expansion coeficients.

To learn more about boundary conditions in HEC-RAS, check chapter 3 of the Hydraulic Reference Manual.

Click OK. Save the flow data as “TuolumneMeadows.f01”. Now we are ready to run the model.

In the HEC-RAS main window, click on Run > Steady Flow Analysis… Click on File > New Plan. Name it “TuolumneMeadows.p01” and “100year” for the short ID. Select the Subcritical Flow Regime, and click on “Compute”.

We have to take into account that the maximum number of station-elevation points is 500. The number of points can be reduced with the XS Points Filter Tool (geometric data window > Tools > Cross section point filter…).

Please note that this is a very simple example for educational purposes only. An actual hydraulic simulation requires the verification of results, identification of possible errors in the data and, if necessary, repetition of calculations. Assessment of the quality of the simulation comes with the knowledge of the study area and experience.

After the computation is finished, we will export the results to ArcGIS to delineate the inundation extent. To export the data to ArcGIS click on File > Export GIS Data… in the HEC-RAS main window. Leave the rest of the parameters as default.

Click on “Export Data” button. A “TuolumneMeadows.RASexport.sdf” file will be created in our working directory. Save the HEC-RAS project and exit. Now we will now return to HEC-GeoRAS to delineate the extent of the flooding.

Open ArcMap and load the “Tuolumne.mxd” file previously created. In the HEC-GeoRAS toolbar, click on “Import RAS SDF file” button to convert the SDF file into an XML file. Browse to “TuolumneMeadows.RASexport.sdf” and click OK.

Next, click on RAS Mapping > Layer Setup. In the Layer Setup window, select “New Analysis”  and name it as “Tuolumne Steady”. Browse to TuolumneMeadows.RASexport.xml for RAS GIS Export File. Select “Single” for Terrain Type, and browse to “tm_dem_flood” GRID. HEC-GeoRAS will create a geodatabase with the analysis name (Tuolumne Steady) in the output directory.

A new data frame named “Tuolumne Steady” will be added to ArcMap.

Click on RAS Mapping > Import RAS Data. This will create a bounding polygon by connecting the endpoints of the XS Cut Lines.

Click on RAS Mapping > Inundation Mapping > Water Surface Generation. Select PF1 as a profile and click OK. A TIN named “t PF 1” is created. This TIN define a zone that will connect the outer points of the bounding polygon.

And now the final step. Click on RAS Mapping > Inundation Mapping > Floodplain Delineation Using Rasters. Select PF1 (100-year flow), and click OK. At this stage the water surface TIN (t PF 1) is first converted to a GRID, and then the terrain surface (tm_dem_flood) is subtracted from the water surface grid. The area with positive results (meaning water surface is higher than the terrain) is flood area, and the area with negative results is dry.

The final result is a polygon (b PF1) which is the final flood inundation polygon, and a GRID (d PF1) that represents the final depth water map. We can see how the maximum water depth is 3.68 meters along the Tuolumne River main channel.

Ok, we are now done with delineating the 100-year flow floodplain in the Tuolumne Meadows area of Yosemite National Park.


Tuolumne Meadows floodplain modeling using HEC-GeoRAS (Part 1)

Tuolumne Meadows is a gentle, dome-studded sub-alpine meadowy section of the Tuolumne River, in the eastern section of Yosemite National Park (California). Yosemite National Park’s hydrologic resources include the headwaters of two magnificent rivers—the Merced and the Tuolumne. All of the creeks, streams, and lakes in Yosemite will eventually join with one of these two rivers.

Toulumne Meadows (Yosemite National Park)

The headwaters of the Tuolumne River arise on the slopes of Mount Dana near Tioga Pass and at the base of Lyell and Maclure glaciers, the two largest glaciers remaining on the west slope of the Sierra Nevada Mountains. The river and its tributaries drain the northern half of the park, passing through deep canyons like Matterhorn, Virginia, Lyell, and the Grand Canyon of the Tuolumne. The Tuolumne watershed is also home to two important California water sources: Lake Eleanor and Hetch Hetchy reservoirs.

We are going to use the HEC-GeoRAS extension of ArcGIS 10 for pre-processing of GIS data for flood inundation mapping. The only dataset required for HEC-GeoRAS is the terrain data (DEM). In our case, we are using a digital elevation model (DEM) of the Tuolumne Meadows area of Yosemite National Park derived from LIDAR that was flown in September, 2006. All above ground features such as vegetation, buildings, and constructed infrastructure has been removed through processing algorithms to leave a flat earth model of the ground surface.

This dataset can be obtained and used for non-commercial purposes from the U.S. National Park Service Integrated Resource Management Applications website: https://irma.nps.gov/

Start ArcMap 10. Save the ArcMap document (File > Save As…) as “Tuolumne.mxd” in your working folder.

Enable ArcGIS’s Spatial Analyst and 3D Analyst extensions by clicking on Customize > Extensions…

Load the HEC-GeoRAS toolbar into ArcGIS (Customize > Toolbars > HEC-GeoRAS).

Create empty GIS layers using the RAS Geometry menu on the HEC-GeoRAS toolbar. Click on RAS Geometry > Create RAS Layers > All.

In the “Create All Layers” window, accept the default names, and click OK. HEC-GeoRAS creates a geodatabase in the same folder where the map document is saved with its same name. In this case, Toulumne.gdb (you may have to set pathnames in ApUtilities > Set Target Locations).

After creating RAS layers, these are added to the map document with a pre-assigned symbology.

Now we use the Editor toolbar (Editor > Start Editing) to digitize the main stream centerline in River layer (IMPORTANT: always digitize from upstream to downstream).

Digitize the main channel of the Tuolumne River.

Click on ID icon to assign River/Reach Codes and type in name the name of the river (in our case, we only have one river segment):

Click on RAS Geometry > Stream Centerline Attributes > Topology. This creates Topology, Lengths/Stations and assign Elevations to the Rivers/Reaches (creates “River3D” layer).

Again, we use the Editor toolbar (Editor > Start Editing) to digitize stream banks in Banks layer. Bank lines are used to distinguish the main channel from the overbank floodplain areas. Always start from the upstream end and, looking downstream, digitize the left bank first and then the right bank.

Click on RAS Geometry > Create RAS Layers > Flow Path Centerlines. Click “Yes” on the message box that asks if you want to use the stream centerline to create the flow path centerline.

Use the Editor toolbar (Editor > Start Editing) to digitize the left and right overbank center flowlines in Flowpaths layer. The left and right flowpaths must be digitized within the floodplain in the downstream direction.

Label the flowpaths by using the “Select Flowpaths and Assign LineType Attributes” button. click on one of the flow paths (left or right, looking downstream), and name each flow path (Left, Channel, Right).

Use the Editor toolbar (Editor > Start Editing) to digitize the cross sections locations for in XSCutlines layer (IMPORTANT: always digitize from left to right looking downstream).

Go to RAS Geometry > XS Cut Line Attributes to populate River/Reach Names, Stationing, Bank Stations, Downstream Reach Lengths and Elevations data for XSCutlines layer (this creates “XSCutlines3D” layer).

After creating cross-sections, the next step is to define bridges and other structures along the river.

Use the Editor toolbar (Editor > Start Editing) to digitize the bridges location for in Bridges layer. In our map, the Tuolumne River is crossed by the California State Route 120 (Tioga Pass Rd) and by the Pacific Crest Trail, as shown in the next picture:

Next, go to RAS Geometry > Bridge/Culverts to assign river/reach names and station numbers and elevations (this creates “Bridges3D” layer).

Ineffective flow area

Use the Editor toolbar (Editor > Start Editing) to digitize Ineffective flow areas. These are areas with water but no velocity within the floodplain. Areas behind bridge abutments can then be considered as ineffective flow areas.

Go to RAS Geometry > Ineffective Flow Areas > Position, to calculate the position of ineffective areas.

The final step before exporting the GIS data to HEC-RAS is assigning Manning’s n value to individual cross-sections.

We delete the LandUse layer (as we created it empty and has no data), and add the “Manning.shp” shapefile. This file store the Manning’s n values in its attribute table (N_VALUE field).

The data come from the Soil Resources Inventory (SRI) for Yosemite National Park area, which meets the geospatial requirements of the Soil Survey Geographic (SSURGO) database.

Click on RAS Geometry > Manning’s n Values > Extract n Values. Select “Manning” for Land Use, choose “N_Value” for Manning Field, “XSCutLines” for XS Cut Lines, name “ManningTable” for XS Manning Table, and click OK.

Finally, to create an import file into HEC-RAS, go to RAS Geometry > Layer Setup > Required Layers tab. The “Required Layers” tab should have River, XSCutLines and XSCutLines3D for Stream Centerline, XSCutmLines and XSCut Lines Profiles, respectively. In the “Optional Layers” and “Optional Tables” tab, make sure the layers and tables that are empty are set to Null, as show in the picture below:

Click on RAS Geometry > Export GIS Data. This process will create two files: “GIS2RAS.xml” and “GIS2RAS.RASImport.sdf”.

That’s it. We are done for now. In later posts I will explain how to import these data into the HEC-RAS hydraulic model.

Geostatistical modelling of Arsenic in groundwater

Arsenic in groundwater is largely the result of minerals dissolving from weathered rocks and soils. Arsenic is found naturally in rocks in the earth’s crust and it is recognized as a poison and a carcinogen substance. It occurs within organic compounds (combined with hydrogen and carbon), and within inorganic compounds (combined within sulphur, chlorine or oxygen). In water, arsenic has no smell or taste and can only be detected through a chemical analysis.

The concentration of arsenic in ground waters in the Duero river basin is relatively low, usually ranging from 0.01 to 0.02 milligrams per litre (mg/L). Sample data (arsenic_sample.shp) for this example comes from the central region of the Duero basin, where a recent investigation has shown increasing concentrations of Arsenic in various linked primarily to illegal agricultural practices.

By using ArcGIS Geostatistical Analyst, we are going to perform a prediction model for the data. In this case the geostatistical method that best suits the spatial variability of the data is the Kriging, which is an stochastic method.

A comparison of the different methods (deterministic and stochastic) can be found in the following document: http://www.esri.com/library/whitepapers/pdfs/geostat.pdf

Open ArcMap. In the menu bar click on Customize > Extensions… and activate Geostatistical Analyst extension. Click close.

In the Geostatistical Analyst toolbar, click on “Launch Geostatistical Analyst Wizard” button. This will prompt the following window:

Methods: We must select the method by which you want to analyze the data, in this case is Kriging.
Source dataset: the shapefile on which we are going to carry out the modelling. In this case, we select “arsenic_sample”.
Data Field: The field you want to perform the analysis. In this case, the concentration levels of Arsenic in mg per litre (“AS_”).

Click next. If you get a message saying that two or more sample point exist at the same location, select the option “Use Mean”. The following window will appear:

In Kriging type, select “Ordinary Kriging”. Select “Prediction” as the output type.
In Transformation, select “Log”, since is necessary for the data to apply a logarithmic transformation.
In Order of trend removal, select “First”. To find out which order of trend removal is necessary for a given dataset, check this link: http://webhelp.esri.com/arcgisdesktop/9.3/tutorials/geostat/Geostat_3_1.htm

Click Next.  The following window allows us to find out whether the data have directional anisotropy present or not. If the surface resemble a circle, there is no directional anisotropy. If there is anything like the figure, the data have directional anisotropy.

In the next window, we choose the following options:

Model: 1. Type = “Spherical”.
Anisotrophy = “True”.
Show Search Direction = “True”.

Now we must change the angle until the blue line that comes from the center of the Semivariogram Map match with the direction of the major axis of the ellipse, as shown in the figure:

Click Next. In the next window, we change the Maximum Neighbours to 12 and Minimum Neighbours to 6. As sampled locations get farther from the prediction location, their influence on the prediction location decreases. So, to speed up interpolation calculations, we can ignore those far-off sampled points.

Click Next. The Cross Validation window will appear. You can see a table showing a the model’s predicted data compared with the measured values, along with the Standard error.

On the lower-right side of the window you can see the prediction errors:

  • Root-Mean-Square: 7.356
  • Average Standard Error: 10,836.42
  • Mean Standardized: 0.002047
  • Root-Mean-Square Standardized: 0.0161

Click Finish. A report summary of the selected method will appear.

The final result is as shown below:

In the final result, we can see how there is a growing gradient in Arsenic concentration from East to West. Note the area colored in red, where the predicted concentration of Arsenic is greater than 15,3 mg/L, which in 300 times higher than the legal limit.

Anyway, we can try to refine the model by changing the model parameters. The lower the Root-Mean-Square (RMS), the better our prediction map is. To compare two models, we can right-click on the model in the table on contents and go to “Compare…”

Additionally, we can calculate the Standard Error Map. A Standard Error Map quantifies the uncertainty of our predictions. The larger the standard errors, the more uncertain are our predictions.