HEC-RAS 5.0 versus TUFLOW versus MIKE21 – HEC’s official response

In case you haven’t seen it, there has been a very popular LinkedIn discussion on the Hydraulic/Hydrologic Modeler’s Forum debating the advantages, disadvantages, and merits of HEC-RAS 5.0, TUFLOW, and DHI’s MIKE21.  There was a lot of great information (and some misinformation) and insight provided in that thread.  It is well worth the read.  Some of the misinformation was directed at the beta version of HEC-RAS Version 5.0.  Enough that HEC decided to publish a response to clear up any confusion.  

At play is a result of HEC-RAS’s 2D solution scheme using the full shallow water equation where for highly dynamic events where flows severely contract, using too small of a time step could lead to a divergence from the true solution, rather than a convergence.  This is not an issue when a computation interval is selected within the guidance presented by HEC in their user’s manual.  While HEC disagrees that this is a necessarily a “problem” with its software, as some on the discussion claim,  HEC has elected to make this a non-issue by “improving the portion of the full shallow water equation formulation, such that user’s will be able to use very small time steps without the results changing significantly.”   I’ve included HEC’s response, written by Gary W. Brunner, to the LinkedIn discussion below, but I highly recommend you read the LinkedIn discussion first by clicking here: 

I’m posting this not just to allow HEC to reach a larger audience with their rebuttal, but also because there is a LOT of great information in this document about how HEC-RAS 2D works and how we, as HEC-RAS 2D modelers should use it.  Please enjoy!  A downloadable pdf is available here.

The following text is copyrighted by the Hydrologic Engineering Center and Gary W. Brunner:

This document is the official HEC response to the discussion on the LinkedIn forum:

Hydraulic/Hydrologic Modeler’s Forum

Discussion Tile:

Tuflow vs Mike DHI products vs HEC-RAS 5

Gary W. Brunner, HEC

            I am the lead Author/Developer of HEC-RAS.  In general, I do not respond to these types of discussions, because I feel it is wrong for software developers to comment about other software products, since they are generally not experts in software products, other than their own.  I also believe it is impossible for any software developer to be completely unbiased in such discussions, as we all think our individual products are the best and we know the most about our own software.  However, there has been so much speculation and misinformation stated about HEC-RAS in this discussion, I feel compelled to respond.

            I read through the discussion (several times), and will try to clear up any misconceptions, and what I consider to be misleading statements about HEC-RAS and its new 2D capabilities.

1.  HEC-RAS solves the full 2D shallow water equations, including Coriolis effects and representation of horizontal turbulent dispersion of momentum using an Eddy viscosity approach.   There was some discussion that some 2D models actually solve the 1D form of the Shallow water equations on a 2D mesh.  HEC-RAS does not do that. It solves the full 2D shallow water equations. 

2.  HEC-RAS uses a semi-Implicit, Eulerian-Lagrangian Finite Volume scheme.  Momentum advection is performed through a Lagrangian tracking step.  Friction, Coriolis, and hydrostatic pressure terms are solved implicitly with second order accuracy in space.  Eddy diffusion is formulated explicitly.  A modified Crank-Nicholson scheme is used in time, where the order of the scheme varies between linear and second order in time depending on the “theta” time weighting coefficient (1.0 is the default which is first order).   The solution scheme is described in full detail in the RAS technical documentation.  All of the details of how HEC-RAS works will be available in Chapter 2 of our technical reference manual, as always.  HEC always tries to be completely open about what equations we use, how we solve the equations, and what assumptions are made.  I believe this is very important in order for a user to truly understand how the software works (if they want to!).

3.  HEC-RAS has been applied to the European Environmental Agency (EA) 2D tests data sets by myself and others.  I have obtained excellent results on all of the data sets (1 -7), but was unable to run Test 8, as we do not currently have the capability to perform a combined 2D overland flow directly connected to a 1D underground pipe network yet. 

4.  There has been a lot of discussion and speculation about HEC-RAS and its sensitivity to the selection of the time step.  All of this discussion has centered around the results for a single test case, and only at a single location within that test case.  This is the Environmental Agency (EA) Test Case No. 5, and the results at only one location for test 5, (Point Location 5).   However, HEC-RAS is more than capable of producing exceptional results for this data set, at all locations (including test point 5), if the user picks a time step based on the guidance provided in the HEC-RAS 2D User's manual.

Here is a plot of HEC-RAS 2D results for Test Case 5, location 5, using multiple time steps of: 12, 10, 6, 5, 4, and 2 seconds:

Here is the plot from the 2013 Environmental Agency report (This plot was extracted from the Environmental agency report from 2013, entitled: “Benchmarking the latest generation of 2D hydraulic modeling packages”  Report – SC120002).  I think everyone involved in this discussion should download this report and review it closely.

            As you can see, the results from the original report have a very wide spread in the results from a range of 2D models.  Two of the models show straight lines, these are from simplistic 2D models that do not solve the full shallow water equations.  However, even for all the models that do solve the full shallow water equations, the results very widely.  If you compare the HEC-RAS results to the results shown at this location you will see that the HEC-RAS results for a wide range of time steps are well within the spread of all the other model results.  If you take away the HEC-RAS results for the 2 second time step run, all of the HEC-RAS results are very close together and fall right on top of the most widely used 2D models (TUFLOW, MIKE21, ISIS 2D, and SOBEK 2D).   All of the RAS results fall within the results shown in the study (even the ones run for the 2s time step).   For this test case (EA Test Case 5), the HEC-RAS results at all of the other locations are even tighter/much more consistent than this location.  Location 5 is the very most downstream location in which model results are being reported in this study, and is the one location where HEC-RAS shows significant variation in the results if you pick extremely small time steps. 

So, if a user follows the HEC-RAS guidance (for the full shallow water equation option) and picks a time step that will produce a Courant number around 1.0 for the high contraction and expansion zones in this data set (it can be lower or higher than 1.0), they will get very good results.   The HEC-RAS time step selection guidance is to find the locations within your model with the highest velocities, then to compute a time step that would produce a Courant number of around 1.0 for that location.

For this test case, the grid resolution specified by the test is 50 meters.  The highest velocities in the model occur at several very tight contractions in which the velocities go up to 5.0 to 6.0 m/s at different locations.  Based on these velocities, the best time step to be used for HEC-RAS would be around 10 seconds.  Why would anyone pick a time step that produces a Courant number of 0.1 or lower for the high velocity portion of the model?  In any practical application of the model they would not, as it would cause their model to run 10 times slower than necessary.  So, in general a user would not select a 1 or 2 second time step for this model if they were following the guidance in the HEC-RAS User’s Manual.

5.  Also shown in the HEC-RAS plot are two runs made with a 100 meter grid.  One with a 10s time step and one with a 15s time step.  As you can see, the results, even at Location 5, are very good.  The HEC-RAS model with the original grid (50 m grid) took 1 min 40s to run Test 5.  However, the 100 m grid with 15s time step only took 19 seconds.  This is one of the major benefits of what we are doing with the subgrid terrain.  The HEC-RAS 2D cells are pre-processed into detailed elevation-volume curves based on all of the information in the underlying terrain.  The elevation profile for each cell is extracted as a detailed cross section based on the underlying terrain.  Next, detailed elevation vs area, wetted perimeter, and roughness curves are developed for each face of each cell.  This allows HEC-RAS 2D to compute very accurate volumes based on the full sub grid terrain.  Because the cell faces are detailed cross sections, HEC-RAS has a very accurate estimate of the terrain under each cell face.  None of the other models mentioned in this LinkedIn discussion do that (to my knowledge at least.  If I am incorrect I apologize to whatever software may have this capability).  The other software packages reduce the elevation data within each computational cell down to a single elevation for each cell and each face is also a single elevation flat line.  Those software packages that use triangles, reduce the terrain down to three elevations (One elevation at the end of each side of the triangle) and each Face is a straight line with two elevations.  This is a huge difference compared to what we are doing in HEC-RAS 2D.

What HEC-RAS does by preprocessing the cells and faces into detailed curves based on the underlying terrain allows the user to use larger cells and still retain the detail of the underlying terrain.  Larger cells means fewer cells, which means faster run times.  As I pointed out above, HEC-RAS was able to run the EA Test Case 5 with a 100 meter grid and get basically the same results as using the 50 meter grid.  But it was able to do this at a speed that was more than 5 times faster.   This is a major feature of HEC-RAS that user need to take advantage of in order to get the most efficiency out of the HEC-RAS software.

6.  Given that the other models discussed in this forum reduce each cell to a single elevation (RiverFlow 2D uses triangles), and each cell face to a single straight line, this in itself could be producing bad results.  How accurate is the estimation of the channel shape, wetted perimeter, and volume underneath the computed water surface if you are only using 5 or so cells across the main channel, with models that reduce each cell to a single elevation.  This could produce a bad estimate for wetted perimeter, cross section shape and volume of the channel/floodplain, and will definitely affect the results the model produces.  Additionally, with these types of models, if you change the cell size, you will get different estimates for the channel shape, wetted perimeter, and volume within the channel and floodplain.  How does one pick a consistent Manning’s n value, if the wetted perimeter and channel shape and volume change when you change the grid resolution?  This means that roughness coefficients picked for these models are tied to the grid resolution, and they will need to be changed when changing grid resolutions.  But no one seems to be talking about the significance of this issue.  HEC-RAS allows the grid resolution to change and still maintain accurate estimates of the channel and floodplain shape, wetted perimeter, and volume underneath the water surface.  This reduces the necessity to have different Manning’s roughness values for different grid resolutions.

Assuming the user understands this issue about these models, and why it is true, they would need to first get a good grid resolution before performing any final calibration of the model.  In areas where there is no data to calibrate, it will require that the user has significant experience in using these models, under conditions where they did have observed data, to gain experience in picking roughness coefficients for those types of models.   Experienced modelers will know and understand this issue, but many user’s do not know or understand this. 

7.  HEC-RAS allows the user to pick time steps that are higher than a Courant number of 1.0.  So our solution has a wide range of possible time steps that will work and produce consistent and good results.  Some of the other models cannot have a time step with a Courant Number higher than 1.0, as they will tend to go unstable.  This is also true of HEC-RAS if you pick way too large of a time step.  This is especially true of software products that use explicit solution schemes, which is what all of the software products that have been released with GPU versions of their software are doing.  So in that sense, HEC-RAS is more flexible in that it allows time steps for Courant numbers greater than 1.0.  But it is currently less flexible when picking a time step that produces Courant numbers much lower than 1.0 (C = 0.1 or lower) in the highest velocity zones of the model (tight contractions and expansions, as well as very sharp bends).

So, why is the current version of HEC-RAS 5.0 Beta sensitive to the time step for a few situations/data sets?

First of all, this issue only comes up for a very small range of data sets in which there are extremely tight constrictions and expansions, and the velocity changes are very extreme (i.e. going from a very low velocity to a very high velocity through the constriction, then back to a low velocity through the expansion).   This is not an issue of high velocity (HEC-RAS can handle high velocity data sets just fine).  Here is an example data set for a river system that is completely straight (i.e. has not contractions/expansions or curves).  The data set is a rectangular channel (200 ft wide) on a constant slope (0.001) that is 10,000 ft long.  The inflow hydrograph goes from 0.0 to 4000 cfs in 10 seconds (like an instantaneous dam or levee breach), then remains constant at 4000 cfs.  The grid size used 40 X 40 ft. grid cells.  So there are 5 cells across the channel.  The maximum velocities vary from a high of around 8 ft/s at the upstream end to around 4.5 ft/s at the downstream end.  Give the 40 ft cell size, and the max velocity of around 8.0 ft/s, a 5 second time step would result in a Courant number of around 1.0 at the upstream end during the max velocity.  All other locations and time will have courant numbers less than 1.0.  I ran this data set with the following time steps: 5, 2, 1, 0.5, 0.2 seconds.  The 0.2 second time step would results in maximum Courant numbers of 0.05 at the upstream end, and even lower everywhere else, as well as when the flood wave is at lower velocities.  The plot below shows the inflow hydrograph and all of the downstream hydrographs for the runs with different time steps.  As you can see, the downstream hydrographs are exactly the same.  This shows that the problem is not an issue of rapidly changing velocity with respect to time, or even distance, if the special velocity changes are not highly curved (as would be the case in severe contractions/expansions, or around sharp bends).

The issue arises for data sets that have extremely high changes in velocity over short distances, and when the path of the velocity is highly nonlinear with respect to space (velocity path is curving over short distances).  However, even in these types of conditions, if the user picks a reasonable time step for these zones, the HEC-RAS model will produce good consistent results.  This is exactly what the EA Test Case 5 shows.  This test case has several very tight constrictions/expansions.  It is also very steep upstream and then becomes very flat downstream.  The velocities in this model are extremely high, in some of the contractions the velocities are between 5 to 6 m/s (16 to 20 ft/s).  Additionally, the models all start out with the downstream area completely dry, then a hydrograph is introduced that goes from zero flow to the peak flow in 5 minutes.  As shown above, if the user picks a reasonable time step with HEC-RAS it can compute a solution that is similar to all of the other listed models in this test case, and is very close to the same results as TUFLOW, Mike21, ISIS 2D, and SOBEK 2D.

The reason HEC-RAS shows variation in the results, when selecting small time steps (far away from a Courant number of 1.0), is due to how we are calculating the convective acceleration terms in the momentum equation.  The current implementation of the 2D “Full Momentum” equations in HEC-RAS 5.0 Beta, uses an Eulerian-Lagrangian formulation (ELM). This is a non-conservative form of the momentum equation which is efficient and accurate for a very wide range of flow conditions.  However, when the convective acceleration terms in the momentum equation are dominant (rapid changes in velocity over short distances), this formulation is accurate within a more limited range of Courant values.

 

When the convective acceleration terms in the equations are dominant, the numerical approximations made in the convective acceleration terms become important. In HEC-RAS 2D the convective acceleration terms are discretized using a Lagrangian formulation which has two steps: a numerical integration on the velocity field, followed by an interpolation process.

If the time-step is too large, the numerical derivatives in the solver will not be appropriate to capture the non-linear behavior of the fluid, driven by high acceleration terms. If the time-step is too small, repeated interpolation to derive the convective acceleration will have a diffusive effect in the velocity field, similar to the smoothing effect of a moving average curve.

This means the user’s need to pay more attention when selecting a time step for the Full Momentum solution for extremely rapidly varying floods that have rapid changes in velocity in space.  Users need to pick a time step that results in a Courant number around 1.0 for the high velocity zones of the simulation.  Picking a time step that is too large (resulting in Courant numbers that are too high, i.e. > 2.0) may result in computational instabilities.  Additionally, for these extremely rapidly varying flow conditions, picking a time step that is way too small (resulting in Courant numbers that are less than 0.1), may result in some numerical diffusion of the computed velocity for certain types of conditions (as described above).

The HEC-RAS Diffusion Wave method does not use the ELM approach in solving the equations.  Therefore this method is not as sensitive to time steps that are too small, and it will ultimately converge on an answer that does not change as you reduce the time step.

Users should always test the consistency of their computational mesh and selected time step.  The consistency principle requires a reduction of both the space (grid) and time steps in order to guarantee convergence of a solution. If the grid is refined and the time-step is reduced simultaneously, the method will achieve convergence.  The user should always test different cell sizes (ΔX) for the computational mesh, and also different computational time steps (ΔT) for each computational mesh.  This will allow the user to see and understand how the cell size and computational time step will affect the results of your model.  The selection of ΔX and ΔT is a balance between achieving good numerical accuracy while minimizing computational time.

Note: people using the current Beta version of HEC-RAS 5.0 will have no trouble getting good results for any data set if they pick a time step based on the guidance provided in the HEC-RAS 2D user’s manual.  So, as far as I am concerned, this is not even an issue for the current Beta release.   However, HEC is currently working on improving this portion of the full shallow water equation formulation, such that user’s will be able to use very small time steps without the results changing significantly.

8.  I read Bill Syme’s contribution to this forum.  First off I would like to say “Thank you” to Bill for being professional and discussing modeling issues in an appropriate professional manner.  I agree with almost all of Bill’s discussion, as I feel that the modeler is the most important thing in producing a good hydraulic model, not the software.  However, I do have a few concerns with Bill’s discussion, as it might imply that HEC-RAS is performing computations in a less robust way than it actually is.  I do not think that Bill Syme was singling HEC-RAS out on these issues, but I do think many people reading his response may think that these comments applied to HEC-RAS.  So I want to set the record straight on a few of his discussion points, as it would pertain to HEC-RAS.  Bill’s comment under Number 2 of his document, does not mention HEC-RAS directly, but it talks about models that solve the 1D shallow water equations on a 2D grid.  HEC-RAS is not doing that!!! HEC-RAS is solving the full 2D Shallow Water Equations with the inclusions of turbulence modeling (Eddie viscosity) and the Coriolis affects.  Point 4 of Bills document talks about 1^st and 2^nd order accuracy.  As I mentioned previously, in general we are 2^nd order accurate in space and time, unless Theta is set to 1.0, then it is first order accurate in time.  Comment 5 also implies that HEC-RAS is 1^st order accurate, which is not correct.

9.  There was a comment in the forum that stated that “HEC-RAS was the least accurate and the slowest”.  I think this is an extremely unfair comment. 

First off, we have made significant speed improvements with each Beta version of HEC-RAS.  So it would depend on which version of HEC-RAS 5.0 2D was used.  Second, the HEC-RAS 2D scheme makes use of multiple CPU cores, so the more cores you have, in general the faster it will run.  So what machine was this run on?  How many Cores did it have?  What was the CPU speeds. 

There are other models using GPU solutions.  Models that use GPU solutions are explicit solvers and require time steps that equate to Courant numbers less than 1.0 for all of space and time during the simulation.  However, GPU’s often have 256 or 512 GPU cores on the graphics card.  For large models (100,000 cells and larger) this can produce much faster run times than models like HEC-RAS that are parallelized based on multiple CPU’s.  A GPU formulation will run faster for models with large numbers of cells, but will tend to require a much finer resolution in time and space.  A fair comparison of HEC-RAS 2D against a GPU code would have RAS running with a courser grid that effectively utilizes subgrid bathymetry, and produce the same results.  Under those circumstances HEC-RAS-RAS 2D may provide accurate results with significantly improved run times.

HEC is in the process of working on a GPU version of RAS 2D, but we are just getting started, so it will be sometime before we have a version that can be run on a GPU.

As to the comment about being “Least Accurate”.   This is based on only running the Environmental Agency Test cases, and only about one of the test case, at a single location within that test case (Test case 5, location 5).  As I stated previously, I have run 7 of the 8 test cases and have obtained just as good as results as any of the models applied to those test cases.  This is another example, of the fact that the modeler is more important than the software product.  If a person knows the software really well, which means knowing its strengths and its limitations, they can produce good results for almost any test case.  Unless the model does not have a required capability to perform that test case.  I will be publishing the official HEC results for the Environmental Agency Test cases shortly after the final release of HEC-RAS 5.0.  I want to wait until the final release so the results show exactly what the user would get if they run RAS 5.0 with the data sets in question.

10.  There was also a comment that “TUFLOW was the most efficient from a project perspective (14 times faster than HEC-RAS)”.  This is also an unfair comment.  Project efficiency is in general not about which model runs faster (unless model run times are exorbitant).  Efficiency is about the experience of the modeler using that specific software product, their knowledge of the river system being modelled, their ability to develop a model that accurately reflects the river/floodplain and the events being modelled, and their ability to understand the model results.  In general, HEC-RAS is very easy to use and efficient.  Developing a 2D model is even easier now with the tools we have added to develop a computational mesh that reflects the unique aspects of the river and floodplain, as well as all of the barriers to flow movement. 

What version of HEC-RAS was used for this time comparison?  How many CPU cores did it have? Which version of TUFLOW were you comparing it against?  Did you run it on the same machines?  Were you using the GPU version of TUFLOW, which would have the benefit of possible 512 GPU cores?  Even if you used the exact same machine, and TUFLOW classic, you failed to realize the benefit of HEC-RAS’s ability to utilize the subgrid terrain data in developing the hydraulic property tables of the cells and the faces, as stated previously.  This means that HEC-RAS will allow users to use far fewer cells than the other models, and still produce similar results.  So if you use HEC-RAS, and realize that it has this capability, you will be able to get similar results with much faster run times.  This was the whole purpose of HEC adding this capability into our 2D software.  By ignoring this significant capability within HEC-RAS, you are not really making a fare comparison of run times.  This was not done for this comparison, but I have done this on many data sets and have shown that HEC-RAS can do this with great success.

Additionally, since HEC-RAS can use unstructured grid cells, users can easily align cells along main channel banks (at the high ground that separates channel from overbank flow); along levees; major roads; and other barriers to flow.  This allows for larger cells, while still accurately defining the critical elevations of these flow barriers with the cell faces.  The other models will have to use much smaller cell sizes to pick up the details of these features.  So how can one ensure that the main channel flow carrying capacity is being accurately respected?  The HEC-RAS subgrid terrain approach, unstructured grids, and breakline tools make it easy for user’s to do this within the full 2D domain.

11.  There was a comment that stated: “The intent of the benchmarking testing was to compare each model against recorded data to determine their accuracy using identical inputs (i.e. Cell size, cell topography, etc…)” This statement is not correct.  First off, none of the EA data sets have observed data to compare against, except Test Case 6A, which is an instantaneous dambreak test from a flume experiment.  However, the observed data from this test case is very noisy (jumps around a lot), and none of the models were able to produce very accurate results compared to this observed data.  Some models did better than others of course.

Second, here is a direct quote from page iv of the EA 2013 Benchmarking Test Results document:

               “The overall objectives of this research are to provide:

            1.  an evidence base to ensure that 2D flood inundation modelling packages used for flood risk management by the Environment Agency and its consultants are capable of adequately predicting the variables on which flood risk management decisions are based

            2.  a data set against which such packages can be evaluated by their developers”

            Please take special note to Item 2 above.  It was always the original intent of the Environmental Agency to see if a model could actually perform these tests, but based on their most appropriate application of the test problem and the software being used.  As you know, it is very easy for someone who is not truly an expert in a particular piece of software, to unintentionally misapply that piece of software.  This is why the Environmental Agency wants the model developers to apply these test data sets and report the results in a public form.  This does not mean that user’s should not also run these tests, or any other tests that they feel are important in their model selection.  However, if you are not getting the same results as HEC on these specific test cases, then you may not be using HEC-RAS to its fullest extent/capability.

12.  There was a comment that stated: “a triangular mesh is better suited to describe a 2D flow behavior because you can compute 3 different flow directions for each cell instead of just 2 as happens for quadrangular meshes”.  This is not correct.  As far as I know, flow can cross a cell face at any angle for all of these models in this discussion.  This is definitely true for HEC-RAS.  HEC-RAS can use triangles, squares, rectangles, and cells up to eight sides.  So the mesh can be structured (square grids) or completely unstructured with cells of varying shapes and sizes.  Water can cross the face of any cell at any angle (360 degrees).  I have performed many test experiments where I have aligned the grid in the direction of flow, then took the same model and purposefully aligned all faces at 30, 45, and 60 degrees to the direction of flow.  HEC-RAS produces almost the exact same results, with very little variation.  This is something we spent a lot of time on to get it right.  Without this being true, no unstructured mesh software would ever be able to solve flow movement for any real world channel and floodplain.

Some final thoughts:

1.       All of these 2D models produce different results for all of the 8 EA test cases.  Some are closer together than others, but none are the same.  Yet they are all solving the 2D shallow water questions (except of few of the models, which are clearly pointed out in the study). 

 

2.      There are many things that will cause different results to occur with these different models, such as:

a.      Mesh size, shape

b.      Time step

c.       How each model depicts the actual terrain within their mesh (Cells and faces).  As I mentioned above.

d.      Friction modeling

e.      Temporal and convective acceleration modeling

f.        Turbulence modeling

g.      Etc..

3.      It is absolutely true that the modeler is more important than the software selected.  A person who is a good hydraulic modeler needs to have the following qualities:

a.      A good understanding of river hydraulics both in theory and practice (real world water movement).

b.      A good understanding of hydrology, as all of these hydraulic models require accurate inflow boundary conditions to have any chance of predicting a good result.

c.       A reasonable understanding of numerical analysis as it applies to solving equations such as the shallow water equations.  This means understanding how grid resolution and time step selection will influence model results.  As well as a reasonable understanding of how the other terms in the momentum equation are being computed, and there significance to your specific study.

d.      A good understanding of the system they are modeling and how it will react to the events they are modelling.

e.      Knowledge of special hydraulic features, such as: levees, bridges, culverts, dams, weirs, spillways, pump stations, etc… and how to best model these features with the software you are using.

4.        None of the test case discussed in this forum have any observed data.  So no one can truly say which model is producing better results.  Many things will affect the results that the models produce, and there are different plusses and minuses of each of the software packages.