a
Computational Technique and 

Performance of Transient Inundation 

Model for Rivers—2 Dimensional 

(TRIM2RD): A Depth-Averaged 

Two-Dimensional Flow Model 

By Janice M. Fulford 

U.S. GEOLOGICAL SURVEY 

OPEN–FILE REPORT 03-371 

Tacoma, Washington 

2003 

i
U.S. DEPARTMENT OF THE INTERIOR 

GALE A. NORTON, Secretary 

U.S. GEOLOGICAL SURVEY 

Charles G. Groat, Director 

Any use of trade, firm, or product names is for descriptive purposes only and 

does not imply endorsement by the U.S. Government 

For additional information write to: Copies of this report can be purchased from:

 

Director, Washington Water Science Center

U.S. Geological Survey 

1201 Pacific Avenue – Suite 600 

Tacoma, Washington 98402 

http://wa.water.usgs.gov

U.S. Geological Survey 

Information Services 

Building 810 

Box 25286, Federal Center 

Denver, CO 80225-0286 

ii
Table of Contents

Abstract............................................................................................................................... 1 

Introduction......................................................................................................................... 1 

Purpose and Scope .......................................................................................................... 2 

Acknowledgments........................................................................................................... 2 

Governing Equations .......................................................................................................... 3 

Numerical Solution ..................................................................................................... 5 

Numerical Method .......................................................................................................... 5 

Wetting and Drying....................................................................................................... 11 

Boundary Conditions.................................................................................................... 11 

Roughness Parameters.................................................................................................. 12 

Model Performance........................................................................................................... 12 

Uniform Depth Scenarios ............................................................................................. 13 

Scenario A: Roughness of 0.035 With Uniform Depth of 2 Meters......................... 14 

Scenario B: Roughness of 0.075 With Uniform Depth of 1Meter ........................... 16 

Scenario C: Roughness of 0.012 With Uniform Depth of 3 Meters......................... 19 

Possible Sources of Simulation Error ....................................................................... 20 

Laboratory Dam-Break Simulation............................................................................... 21 

Reach Downstream of Dam...................................................................................... 21 

Pool Behind Dam and Reach Downstream of Dam.................................................. 24 

Large-Scale River Simulations..................................................................................... 28 

Chattahoochee River Simulation .............................................................................. 28 

Snoqualmie River Simulation................................................................................... 33 

Summary and Conclusions ............................................................................................... 40 

References......................................................................................................................... 41 

Appendix A. File Specifications for Model Input and Output Files................................. 43 

 

Table of Contents 4
Figures 

Figure 1. Diagram showing relation between vertical datum and variables 

representing depth and water-surface elevation.................................................4 

Figure 2. Diagram showing computational mesh used in spatial discretization ...............6 

Figure 3. Diagram showing computational module for TrimR2D, showing how 

terms from non-adjacent cells are included in the computational 

node at (i,j).......................................................................................................10 

Figure 4. Graph showing schematic of a simple mesh and boundary conditions...........12 

Figure 5. Graph showing simulated water-surface profiles for scenario A with a 

roughness of 0.035 and a uniform depth of 2 meters.......................................15 

Figure 6. Graph showing simulated flow depth for scenario A with a roughness 

of 0.035 and a uniform depth of 2.0 meters.....................................................16 

Figure 7. Graph showing simulated water-surface profiles for scenario B with 

a roughness of 0.075 and a uniform depth of 1 meter .....................................17 

Figure 8. Graph showing simulated flow depth for scenario B with a roughness 

of 0.075 and a uniform depth of 1.0 meter ......................................................18 

Figure 9. Graph showing simulated water-surface profiles for scenario C with 

a roughness of 0.012 and a uniform depth of 3 meters....................................19 

Figure 10. Graph showing simulated flow depth for scenario C with a roughness 

of 0.012 and a uniform depth of 3.0 meters.....................................................20 

Figure 11. Hydrographs for site located 21.335 meters upstream of the outflow 

boundary ..........................................................................................................22 

Figure 12. Graph showing simulated water-surface profiles of reach downstream 

of dam at 1, 15, 30, and 55 seconds elapsed from start of simulation .............23 

Figure 13. Hydrographs for site located 36.58 meters upstream of flume end .................25 

Figure 14. Graph showing measured and simulated water-depth profiles for pool 

behind dam and reach downstream of dam at various elapsed times..............26 

Figure 15. Graph showing simulated water-surface profile for initial conditions 

for the modeled reach of the Chattahoochee River..........................................30 

Figure 16. Hydrographs for March 22-23, 1976, at Georgia Highway 120 for 

measured and computed water-surface elevations for 2- and 60-second 

time-step sizes..................................................................................................31 

Figure 17. Graph showing simulated water-surface profiles of Chattahoochee River 

at 30, 510, 990, 1,500, 1,740, and 2,190 minutes on March 22-23, 1976 ......32 

Figure 18. Hydrographs for the November 24, 1986, flood event for Snoqualmie 

River stream-gaging stations............................................................................35 

 Figures v
Figure 19. Map showing simulated peak water-surface elevations on the Snoqualmie 

River for the November 24, 1986, flood event for a Manning's 

roughness value of 0.12 ...................................................................................36 

Figure 20. Hydrographs of the December 3, 1975, flood event for Snoqualmie River 

stream-gaging stations .....................................................................................38 

Figure 21. Map showing simulated peak water-surface elevations on the Snoqualmie 

River for the December 3, 1975, flood event for a Manning's roughness 

value of 0.12.....................................................................................................39 

 Figures vi
Tables 

Table 1. Uniform depths for three scenarios representing a wide range of stream 

channel materials .............................................................................................14 

Table 2. Mean and standard deviation for percent error of simulated flow depths 

as a function of time.........................................................................................28 

Table 3. Mean and standard deviation for percent error of simulated flow depths 

as a function of flume location.........................................................................28 

Table 4. Comparison of measured and simulated peak water-surface elevations for 

calibration of the Snoqualmie River model, using the flood event of 

November 24, 1986, and a Manning’s roughness value of 0.12......................37 

Table 5. Comparison of measured and simulated peak water-surface elevations for 

verification of the Snoqualmie River model, using the flood event of 

December 3, 1975, and a Manning’s roughness value of 0.12........................40 

vii Tables 
Conversion Factors and Abbreviations 

Multiply By

Length To obtain 

centimeter (cm) 0.3937 inch (in.) 

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi) 

kilometer (km) 0.5400 mile, nautical (nmi) 

Flow rate 

square meter per second (m2/s) 10.76 Square foot per second (ft2/d) 

meter per second (m/s) 3.281 foot per second (ft/s) 

Conversion Factors and Abbreviations viii
COMPUTATIONAL TECHNIQUE AND 

PERFORMANCE OF TRANSIENT 

INUNDATION MODEL FOR RIVERS — 

2 DIMENSIONAL (TRIMR2D): A DEPTH–

AVERAGED TWO–DIMENSIONAL FLOW 

MODEL 

By Janice M. Fulford 

ABSTRACT 

A numerical computer model, Transient Inundation Model for Rivers – 2 

Dimensional (TrimR2D), that solves the two-dimensional depth-averaged flow equations 

is documented and discussed. The model uses a semi-implicit, semi-Lagrangian finitedifference method. It is a variant of the Trim model and has been used successfully in

estuarine environments such as San Francisco Bay. The abilities of the model are 

documented for three scenarios: uniform depth flows, laboratory dam-break flows, and 

large-scale riverine flows. The model can start computations from a “dry” bed and 

converge to accurate solutions. Inflows are expressed as source terms, which limits the 

use of the model to sufficiently long reaches where the flow reaches equilibrium with the 

channel. The data sets used by the investigation demonstrate that the model accurately 

propagates flood waves through long river reaches and simulates dam breaks with abrupt 

water-surface changes. 

INTRODUCTION 

River flooding is a significant natural hazard. The 1993 Mississippi River flood 

caused losses worth 15 to 20 billion dollars and 48 deaths (Perry, 2000). Flood-response 

planners need flow measurements, estimates of inundated areas, and prediction of river 

flows throughout a river valley to plan remediation actions and public-safety responses. 

Usually, measurements of river discharge and stage and predictions of flood stage from

runoff calculations are limited to a few locations in the river valley. Stage and flow 

predictions at a few additional locations typically are done with one-dimensional storagerouting models. These models do not adequately account for inertia or backwater effects. 

The depth and extent of inundation at other locations in the river valley must be 

approximated from the available data, either measured or predicted. The paucity of flood 

information can hamper mitigation and public-safety decisions. 

 Introduction 1
Methods are available that can increase the amount of flood information available for 

decision-making. Numerical computer models that solve the shallow-water equations can 

be used to simulate flow and provide information at locations where no measurements are 

available. Two-dimensional depth-averaged models can be used to predict flood 

inundation throughout the river valley from forecast flows at point locations (Bates and 

DeRoo, 2000). 

The U. S. Geological Survey (USGS) is investigating a numerical flow model, 

Transient Inundation Model for Rivers - 2 Dimensional (TrimR2D), for providing nearreal-time flood information on the Internet. The model will be used to calculate the extent 

of flood inundation and flow features such as stage and velocity distributions from 

forecast flood hydrographs from the National Weather Service. 

TrimR2D is a two-dimensional numerical flow model that is extensively based on 

the Trim model (Casulli, 1990; Cheng and others, 1993) written by V. Casulli of the 

University of Trento, Italy. The Trim model initially was developed for estuarine 

applications and used to simulate flows in San Francisco Bay (Cheng and others, 1993). 

The outstanding numerical stability characteristics of the Trim2D model made its 

adaptation to riverine modelling appealing. Trim2D was adapted to river systems by 

modifying the input and output routines and by modifying boundary conditions for the 

needs of river modeling. The model has been used to model flows in steep river systems 

(Ostenaa and others, 1999). 

Some work using the model has been published. Test cases confirming the lack of 

spurious modes and the second-order convergence rate for uniform grid refinement were 

performed by Walters and Denlinger (1999). However, no confirmation of model results 

with measured data was published for the river simulation described by Walters and 

Denlinger (1999). Examples comparing results from the TrimR2D parent code with 

measurements have been published for oceanographic conditions (Cheng and others, 

1993). 

Purpose and Scope 

This report describes the computational technique and demonstrates the performance 

in river systems of the TrimR2D computational flow model. The description of the 

computational technique includes a presentation of the governing equations and the 

numerical solution used to solve the governing equations. Model performance for river 

applications is demonstrated and discussed for three types of scenarios: (1) uniform-depth 

flows, (2) laboratory dam-break flows, and (3) large-scale riverine flows. The format of 

model input and output files is documented in the appendix. 

Acknowledgments 

The assistance of Dr. Roy Walters (U.S. Geological Survey, retired) is greatly 

appreciated. His advice and interest in the study has made this report possible. 

 Introduction 2
GOVERNING EQUATIONS 

The governing equations for TrimR2D are commonly known as the shallow-water 

equations. The shallow-water equations describe two-dimensional, depth-averaged, 

unsteady flow and are based on the conservation of mass and momentum. The equations 

assume that water is incompressible and that pressure distribution is hydrostatic. Density 

is constant throughout the water column and vertical velocities are considered small 

relative to the horizontal velocities. In differential form, the governing equations of

continuity (eq. 1) and momentum (eqs. 2 and 3) conservation that are solved by TrimR2D 

have the form 

[ ] ( ) [ ] ( ) = 0 ∂

∂ +

+ ∂

∂ +

+ ∂

∂

y

h v

x

h u

t

ζ ζ ζ , (1) 

( ) ( ) 1 2

h

w b

x x

u u u u v g A u fv

t x y x h

ζ τ τ

ρ ζ

∂ ∂ ∂ ∂

+ + = − + − + ∇ + ∂ ∂ ∂ ∂ + , (2) 

and 

( ) ( ) 1 2

h

w b

y y

v v v u v g A v

t x y y h fu ζ τ τ

ρ ζ

∂ ∂ ∂ ∂

+ + = − + − + ∇ − ∂∂∂ ∂ + , (3) 

where t is time; u(x,y,t) and v(x,y,t) are the depth-averaged velocity components in the xand y-axis directions in the horizontal plane, respectively; ζ(x,y,t) is the water-surface 

elevation measured from the vertical datum; h(x,y) is the depth measured from the 

vertical datum; g is the acceleration of gravity; f is the Coriolis parameter; b

x τ and b

y τ are 

the bottom-stress terms; w

x τ and y

ω

τ are the wind-stress terms in the x- and y-axis 

directions; ∇ = i ∂ ∂x + j∂ ∂y G G is a vector operator in the horizontal x-y plane; and Ah

 is 

the horizontal eddy viscosity. The vertical datum is arbitrary and is selected to minimize 

numerical truncation errors. Usually, for oceanographic applications, mean sea level is 

the vertical datum. For river applications, the approximate outflow elevation is the 

vertical datum. Figure 1 illustrates the relation between the vertical datum and variables

that represent depth and water-surface elevation. The water depth, H, is equal to (h+ζ). 

The depth, h(x,y), measured from the vertical datum, is positive in the down direction and 

negative in the up direction. The depth (or distance to the land surface) is elevation of the 

topography referenced to the vertical datum. The water-surface elevation, ζ(x,y,t), 

measured from the vertical datum, is positive in the up direction and negative in the down 

direction. 

 Governing Equations 3
Figure 1. Relation between vertical datum and variables representing depth and water-surface 

elevation 

The bottom-stress terms are formulated using a Manning-Chezy type expression as 

( ) u

h

b

x τ γ

ρ ζ = +

1 (4)

and 

( ) v

h

b

x τ γ

ρ ζ = +

1 , (5) 

where 

γ = ( )

2 2

2

g u v

C h ζ

+

+ (6)

is the bed friction factor. 

C is the Chezy coefficient. Manning's coefficient in metric units is related to Chezy 

coefficient by

( ) 1

h 6

n

C

+ζ = . (7) 

 Governing Equations 4
NUMERICAL SOLUTION 

The shallow-water equations are not amenable to direct solution and require that the 

continuous derivatives in the equations be approximated by finite-difference techniques. 

Finite-difference techniques solve the governing equations for a finite number of 

locations in space and time. The techniques require the subdivision of the application 

domain into a mesh (or grid) with a finite number of node points. 

The numerical solution of the governing equations uses a semi-implicit finitedifference method that approximates the advective terms of the momentum equations 

with a Lagrangian formulation. Simulation of the topography wetting and drying is 

permitted by the solution algorithm. The following is a brief overview of the numerical 

method, wetting and drying computations, boundary conditions, and roughness 

parameters. Details on the derivation and stability of the numerical technique are 

documented in Casulli (1990) and Casulli and Cattani (1994). 

Numerical Method 

The equations are spatially discretized as a staggered rectangular mesh that consists 

of square computational cells with length ∆x and width ∆y (∆x=∆y). Figure 2 illustrates 

the computational mesh used in the spatial discretization. Water-surface elevation, ζ , is 

defined at cell centers with an index of (i,j), and velocities u,v are defined at the middle of 

cell sides with an index of (i+1/2,j+1/2). The x-axis velocities, 1

2 i , u + j and 12 i , j , are 

defined at the right and left sides of the cell, and the y-axis velocities, 

u −

1

2 i j , v + and 12 i j , v − , 

are defined at the upper and lower sides of the cell. An arithmetic average of the 

surrounding values is used for variables not defined at a spatial location. 

Depth or land topography, h, is defined at the velocity locations. However, land 

topography is input to TrimR2D at node and half-node mesh locations. The resulting 

topography mesh has a resolution that is twice that of the computational mesh. The 

higher resolution of the topography mesh eliminates the need for users to construct a 

staggered topography mesh and eliminates a potential source of user error. 

 Governing Equations 5
Figure 2. Computational mesh used in spatial discretization 

The temporal discretization is a mixture of explicit and implicit methods. The 

gradients of water-surface elevation (

x

g ∂

∂ζ , y

g ∂

∂ζ ) in the momentum equations (eq. 2

and 3) and the velocity divergence (

x

u

∂

∂ , y

v

∂

∂

) in the continuity equation (eq. 1) are 

discretized semi-implicitly to remove the stability constraint for gravity waves. The 

velocity part of the friction terms (γu , γv ) in the momentum equations is discretized 

implicitly. The γ parts of the friction terms in the momentum equations and the water 

depth terms (h+z) in the continuity equation are discretized explicitly. 

 Governing Equations 6
The advective terms (

y

v v

x

v u

t

v

∂

∂ + ∂∂ + ∂

∂ , y

u v

x

u u

t

u

∂

∂ + ∂∂ + ∂∂ ) in the momentum 

equations are expressed as a substantial (or material) derivative, D/Dt, resulting in 

( ) ( ) 1 w b 2

x x h

Du

g u fv

Dt x h A ζ

ρ ζ τ τ

∂ = − + − + +

∂ + ∇ (8) 

and 

( ) ( ) 1 w b 2

y y h

Dv

g v fu Dt y h A ζ

ρ ζ τ τ

∂ = − + − + − ∂ + ∇ . (9) 

Equations (8) and (9) are discretized in time with a forward difference operator. 

The semi-implicit discretization of equations (1), (8), and (9) have the form 

( 1 1 ( )) ( 1 ( )) 2 2 2

1 1

, , , 1, , 1, , , ma , x , max , j

k k k k k k k k

i j i j i j i j i j i j i j i j i j i

t h u h x ζ ζ θ ζ ζ ζ ζ + +

+ + + − −

∆ ⎡ ⎤ = − + − +

∆ ⎢ ⎥ ⎣ ⎦ 1

2,

1 u +

−

( 1 1 ( )) ( 1 ( )) 2 2 2

1 1

, , max , , 1 , , , 1 , max , k k k k k ki j i j i j i j i j i j i j i j h v h

y

θ ζ ζ ζ ζ + +

+ + + − −

⎡ ⎤ − + − + ⎢ ⎥ ∆ ⎣ ⎦ 1

2 ,

t

v −

∆ 

( 1 1 ( )) ( 1 ( )) 2 2 2 , 1, , 1, , , , (1 ) max , max , k k k k k k

i j i j i j i j i j i j i j i j h u h x θ ζ ζ ζ ζ + + + − −

⎡ ⎤ − − + − + ⎢ ⎥ ∆ ⎣ ⎦ 1

2,

t u −∆

( 1 1 ( )) ( 1 ( )) 2 2 2

1 1

, , , , 1 , , , 1 (1 ) max , max , k k k k k ki j i j i j i j i j i j i j i j h v h

y

θ ζ ζ ζ ζ + +

+ + + − −

⎡ ⎤ − − + − + ⎢ ⎥ ∆ ⎣ ⎦ 1

2 ,

t

v −

∆ (10) 

( ) ( ) 1 1 1 1 2 2 2 2 1

2

1 1 1 1

1, , , , , ,

,

1

k

k k k k k k w

i j i j i j i j i j i j x

i j

t u Iu g t u fv

x h θ ζ ζ γ τ

ρ ζ

+ + + +

+ + + + +

+

⎛ ⎞ ∆ ⎛ ⎞ ⎜ ⎟ = − − − ∆ − ⎜ ⎟ +

∆ + ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

 

( ) 1, , (1 ) k k

i j i j g x θ ζ ζ − − + − ∆

∆t (11) 

( ) ( ) 1 1 1 1 2 2 2 2 1

2

1 1 1 1

, 1 , , , , ,

,

1

k

k k k k k k w

i j i j i j i j i j i j y

i j

t v Iv g t v fv

y h θ ζ ζ γ τ

ρ ζ

+ + + +

+ + + + +

+

⎛ ⎞ ∆ ⎛ ⎞ ⎜ ⎟ = − − − ∆ − ⎜ ⎟ − ∆ + ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

( ) 1, , (1 ) k k

i j i j g

y

θ ζ ζ − − + − ∆

∆t . (12) 

In equation 10, terms similar tomax ( , , 1, ) k k ζ ζ i j i+ j

 denote that the maximum value between 

the parentheses is selected and used in the equation. The variable θ is a temporal 

difference operator. I is a nonlinear finite-difference operator that only contains terms at 

the previous time step. The eddy viscosity terms in the momentum equations have been 

neglected for the sake of simplicity, but are included in the actual formulation used by

TrimR2D with an explicit discretization and hence a weak stability constraint. 

 Governing Equations 7
The momentum equations (11) and (12) are substituted into the continuity equation (10). 

The resulting equation is solved for water-surface elevation at the new time step and has 

the form 

( ) ( ) 1 1 1 1 2 2 2 2

1 1

2 2

2 , , , , 1 2 1 1 1 1

, , 1, 2

, , 1 1

k k

i j i j i j i j k k k

i j i j i j i j i j k k

i j i j

h h t

g

x t t

ζ ζ

ζ θ ζ ζ ζ ζ

γ γ

+ + − − + + +

+ −

+ −

, 1,

k+ k+

⎡ + + ⎤ ∆

+ − ⎢ + − ⎥

∆ + ⎢ ∆ + ∆ ⎥ ⎣ ⎦

( ) ( ) 1 1 1 1 2 2 2 2

1 1

2 2

2 , , , , 2 1 1

, , 1 , , 1 2

, , 1 1

i j i j i j i j k k k k

i j i j i j i j k k

i j i j

h h t

g

y t t

ζ ζ

θ ζ ζ ζ

γ γ

+ + − − + + + +

+ −

+ −

1 1

k k

ζ

⎡ + + ⎤ ∆

+ − ⎢ + − ⎥

∆ + ⎢ ∆ + ∆ ⎥ ⎣ ⎦

( ) ( ) 1 1 1 2 2 2

1

2

, , , ,

,

1 1

1

k k k

i j i j i j i j k

i j

t I u h

x t

ζ ζ θ

γ + + +

+

⎡ ⎛ ⎞⎤

∆ ⎢ ⎜ ⎟⎥ = − + ⎢ + ⎜ ⎟ + ⎥ ∆ ⎢ + ∆ ⎥ ⎣ ⎝ ⎠⎦

1 1 ( 1 ) 2 2 2 1

2

, , ,

,

1 1

1

k k

i j i j i j k

i j

t I u h

x t

ζ θ

γ − − −

−

⎡ ⎛ ⎞⎤ ∆ ⎢ ⎜ ⎟⎥ + + + +

∆ ⎢ ⎜ ⎟⎥ + ∆ ⎢⎣ ⎝ ⎠⎥⎦

( )

( ) ( ) ( )( ) ( ) 1 1 1 1 2 2 2 2

1 1

2 2

2 , , , ,

1, , 1, ,

, ,

(1 )

1 1

k k

i j i j i j i j k k k k

i j i j i j i j k k

i j i j

h h

t

g

x t t

ζ ζ

θ θ ζ ζ ζ ζ

γ γ

+ − + −

+ −

+ −

⎧ ⎫ + + ⎛ ⎞ ∆ ⎪⎪ ⎪⎪ − − ⎜ ⎟ ⎨ − − − ⎬ ⎝ ⎠ ∆ ⎪ + ∆ + ∆ ⎪ ⎪⎩ ⎭⎪

1 1 ( 1 ) 2 2 2 1

2

, , ,

,

1 1

1

k k

i j i j i j k

i j

t I u h

y t

ζ θ

γ + + +

+

⎡ ⎛ ⎞⎤ ∆ ⎢ ⎜ ⎟⎥ − + + +

∆ ⎢ ⎜ ⎟⎥ + ∆ ⎢⎣ ⎝ ⎠⎥⎦

1 1 ( 1 ) 2 2 2 1

2

, , ,

,

1 1

1

k k

i j i j i j k

i j

t I u h

y t

ζ θ

γ − − −

−

⎡ ⎛ ⎞⎤ ∆ ⎢ ⎜ ⎟⎥ + + + +

∆ ⎢ ⎜ ⎟⎥ + ∆ ⎢⎣ ⎝ ⎠⎥⎦

( )

( ) ( ) ( )( ) ( ) 1 1 1 1 2 2 2 2

1 1

2 2

2 , , , ,

, 1 , , 1 ,

, ,

(1 )

1 1

k k

i j i j i j i j k k k k

i j i j i j i j k k

i j i j

h h

t

g

y t t

ζ ζ

θ θ ζ ζ ζ ζ

γ γ

+ − + −

+ −

+ −

⎧ ⎫ + + ⎛ ⎞ ∆ ⎪ ⎪ − − ⎜ ⎟ ⎨ ⎬ − − −

⎝ ⎠ ∆ ⎪ ⎪ + ∆ + ∆

⎩ ⎭

,(13) 

where k

ζ is the average water elevation over a mid-side node. The Coriolis and wind 

terms have been neglected in equation (13) to simplify the presentation of the discretized 

equations. These terms are included in the formulation used by the model, however they 

are not needed for river applications. Equation (13) is a linear penta-diagonal system of

equations. The resulting matrix is positive definite and has a unique solution when 

( ) 0 , , 2

1 2 ± 1 + i± j ≥ k

i j ζ h AND ( ) 0 2

1 21 , ± + i, j ± ≥ k

i j ζ h . 

 Governing Equations 8
The matrix is solved using a conjugate-gradient method that is suitable for vector 

computations. 

In order to obtain an unconditionally stable and relatively accurate form for I, an 

Eulerian-Lagrangian approximation is used. The substantial derivative of a general 

variable, w, is 

t

w

v

t

w

u

t

w

Dt

Dw

∂

∂

+

∂

∂

+

∂

∂ = , (14) 

where D/Dt represents the time rate of change along a streamline defined by 

u

dt

dx = , v

dt

dy = . (15) 

As previously shown in equations 11 and 12, the substantial derivative can be discretized 

in time as 

1

, ,

k k

i j i j I Dw

Dt t

w w + −

≈ ∆ . (16) 

The term is the value of w at time t I ,

k

wi j k located at (i-a,j-b) which is convected in the 

grid in an elapsed time ∆t. Thus, 

1

, ,

k k

Dw i j i a j b

Dt t

w w +

− − −

≈ ∆ , (17) 

where 

x

t

a u

∆

∆ = , y

t b v

∆

∆ = . (18) 

Because a and b are not usually integers, bilinear interpolation is used to approximate the 

current value of w at tk. Letting a = n+p and b=m+q where n and m are integers and 

0 ≤ p < 1, and 0 ≤ q <1, ,

k

i j Iw is approximated by 

( ) ( ) ( ) , , , , 1 1, 1, 1 1 1 k k k k k k

i j i a j b i n j m i n j m i n j m i n j m Iw w− − p q w− − qw− − − p q w− − − qw− − − − = = − ⎡ ⎤ − + + ⎡ − + ⎣ ⎦ ⎣ 1

⎤

⎦. 

(19)

The computational module for TrimR2D is diagrammed in figure 3. It illustrates how

equation 19 includes terms from cells that are not necessarily adjacent to the 

computational node at (i,j). 

 Governing Equations 9
Figure 3. Computational module for TrimR2D, showing how terms from non-adjacent cells are 

included in the computational node at (i,j) 

Because u and v are not usually constants, the correct value of a and b is found from 

numerically integrating the ordinary differential equations in equation (17) from t0 to t1

by the Euler method (Conte and DeBoor, 1980). The time step, ∆t, is divided for the 

integration method into N equal time parts, so that ∆t' = ∆t/Ν . The time subdivision is 

restricted to be 

1 1

2 2 i,j i,j , ,

t' min ,

max max k k

i j i j

x y

u v + +

⎡ ⎤

⎢ ⎥ ∆ ∆ ∆ ≤

⎣ ⎦

 . (20) 

 Governing Equations 10
This restriction prevents the approximated streamlines from crossing solid boundaries. 

This condition is not required for the stability of the method, but improves the accuracy. 

The overall mesh shape for any application of TrimR2D is a rectangle. However, 

computation does not necessarily occur for every cell in the mesh. The model ignores the 

computation of cells with topography higher than a user-defined elevation. This reduces 

the solution matrix size and permits some flexibility in representing areas that are 

permanently dry. 

Wetting and Drying 

The discretized equation (13) inherently permits wetting and drying of cells without 

using complicated bookkeeping procedures. If the water elevation in a wet adjacent cell is 

higher than the topography at the "half" node, flow will enter the cell from the adjacent 

cell. To prevent the introduction of non-physical, negative values of total flow depth, H is 

limited to be max (0, h+ζ). When H is zero, the respective bed friction factor, γ, is 

assumed large. 

Boundary Conditions 

Two types of boundary conditions are available in the model, an external watersurface elevation boundary and an internal source. Both boundary conditions are time 

varying and exploit the properties of the discretized equation (13). The inclusion of these 

boundary conditions allows the solution matrix to maintain a penta-diagonal form with no 

off-diagonal terms. 

The external water-surface elevation boundary sets a spatially uniform water 

elevation on the perimeter of the model mesh. To apply the external water-surface

boundary, the model adds boundary cells to the perimeter of the user-defined mesh. In the 

boundary cells the water-surface elevation can be specified by either a time series or a 

harmonic function. Flow from the boundary cells occurs when the topography and water 

elevation of an abutting cell are lower than the water-surface elevation in the boundary 

cells. Flow out of the model occurs when the water-surface elevation in the boundary

cells is less than the water-surface elevation in an abutting wet cell. 

The internal boundary condition, point sources, is the only method in the model 

available to add flow. Flow is introduced to the model mesh by the addition of a source 

term to the continuity equation 

[ ] ( ) [ ] ( ) q

y

h v

x

h u

t = ∂

∂ +

+ ∂

∂ +

+ ∂

∂ζ ζ ζ , (21) 

where q is the discharge. The source term is expressed as on the right side of the 

discretized equation for water-surface elevation (eq. 13). The location of each source is 

specified by individual cell locations in the mesh. Flow for a cell containing a source 

(source cell) is specified by a time series of discharge values. Flow out of the source cell 

follows the rules for wetting and drying. Thus, when the water elevation in a source cell 

is higher than the water elevation and topography in an adjacent cell, flow out of the 

source cell occurs. The method for adding flow is analogous to a spring flowing 

1

,

k +

i j q

 Governing Equations 11
vertically into a water body. The method does not account for the momentum effects that 

occur when water flows into a water body with a non-zero horizontal velocity. Hence, the 

simulation of river inflows with source terms will result in simulation errors near the 

location of the inflow. 

The effect of the boundary conditions on the model is analogous to an island in an 

ocean. The external water-surface elevation boundary is the level of the ocean that 

surrounds the island. The sources are the springs on the island and they supply the river 

flow. A schematic of a simple mesh and the boundary conditions is shown in figure 4. 

Figure 4. Schematic of a simple mesh and boundary conditions 

Boundary cells are added by the program 

Roughness Parameters 

Roughness parameters for TrimR2D are Manning's n. Values of Manning's n can be 

varied as a function of flow depth, but cannot be varied as a function of location in the 

mesh. The model is not able to simulate flow patterns that result from spatially varying 

the roughness associated with differences in land use, soil, or vegetation types. 

MODEL PERFORMANCE 

This section demonstrates and discusses TrimR2D performance for typical river or 

open-channel flow type applications. No attempt is made to rigorously demonstrate 

spatial and temporal convergence of the solution by repeatedly increasing the resolution

of the applications. It is assumed that the numerical scheme used and described in the 

previous sections is convergent (Cassulli and Cattani, 1994). Model simulations for three 

types of conditions are presented: (1) uniform-depth flows, where gravitational forces are 

in equilibrium with frictional forces; (2) laboratory-measured dam-break flows that are 

highly advective; and (3) large-scale (greater than a kilometer in length) river flows. 

Simulation results are compared with well-known channel resistance equations and with 

measured data as appropriate. The applications herein demonstrate the model’s capability 

to simulate timing and magnitude of the water-surface elevation and peak water-surface 

 Governing Equations 12
elevations for a limited number of locations within each application. Because of a lack of 

measured field data on flooded extent, no simulations demonstrate the performance of 

TrimR2D in determining flooded extent. 

Uniform Depth Scenarios 

Uniform depth is a concept familiar to hydraulic engineers. It occurs when a constant 

flow in a channel of uniform shape and slope is in equilibrium with gravitational forces. 

When equilibrium is reached, the flow depth is constant and the water-surface slope 

equals the slope of the channel. Simulating uniform depth demonstrates how the 

TrimR2D approximation of bed friction in riverine systems affects the model solution of 

the governing equations. 

The uniform flow concept is expressed in the well-known engineering equation, 

Manning’s equation (Chow, 1959) 

2 1 3 2 k Q AR S n = , (22) 

where Q is the total discharge, n is the Manning’s roughness coefficient, A is the crosssectional area, R is the hydraulic radius, S is the channel slope (energy slope, equivalent 

to channel slope for uniform depth simulations), and k is a conversion factor that is 1 for 

meter-second units and 1.49 for feet-second units. The hydraulic radius is the crosssectional area divided by the wetted perimeter of the channel cross section. It is 

approximately equal to the flow depth when a channel is wide relative to the flow depth. 

The equilibrium flow depth is called the uniform or “normal” depth and is different for 

each flow rate for a given channel. The Manning’s equation, though somewhat empirical, 

has a long record of closely approximating uniform-flow depth in river systems. 

Originally determined by fitting measured data, the Manning's equation can be derived 

from the Navier-Stokes equations. 

The uniform depth simulations were for a 1,000-meter-long straight channel with a 

4-meter-wide rectangular cross section and a bed slope of 0.005. The topographic mesh 

(representing depth, h, in the governing equations) has half-meter cell resolution that 

result in 1-meter-sized computational cells. The computational mesh is 6 by 1,001 cells. 

Cell elevations on three sides of the computational mesh prevent flow through those sides 

for all the uniform depth simulations and restrict outflow to occur out one end of the 

channel. 

Uniform depths were simulated for three scenarios (table 1) that have different 

discharges and roughnesses representing a wide range of channel materials. Boundary 

conditions, downstream water-surface elevation, and upstream source discharge are 

constant for all of the uniform-depth simulations. Discharges were computed using 

equation (22). However, the hydraulic radius used in the discharge equation was 

 Model Performance 13
computed by neglecting the contribution of the channel side walls to the wetted 

perimeter. This method of computation accounts for the omission by TrimR2D of friction 

losses due to vertical side walls. 

At least two different downstream water-surface elevation boundary conditions were 

simulated for each scenario. For scenarios with subcritical flow conditions, backwater 

and drawdown water-surface profile curves were simulated. These water-surface profile 

curves are also referred to as M1 and M2 curves, respectively, by many hydraulic 

engineering texts (Chow, 1959). Backwater curves (M1) occur for subcritical flow when 

the downstream flow depth is larger than normal depth. Drawdown curves (M2) occur for 

subcritical flow when the downstream flow depth is smaller than normal depth. For both 

types of curves, the flow will converge toward normal depth in the upstream direction. 

Table 1. Uniform depths for three scenarios representing a wide range of stream 

channel materials 

Scenario Manning's

roughness 

Normal depth 

(meters) Froude No. 

Discharge 

(cubic meters 

per second) 

A 0.035 2.0 0.72 25.7

B 0.075 1.0 0.30 3.77

C 0.012 3.0 2.26 147

Simulations were for run times that usually allowed the computations to reach 

steady-state (equilibrium) flow conditions. All simulations began with a “dry” bed 

condition. Most two- and one-dimensional flow models cannot be started with all 

computational nodes having a zero flow depth. Usually, the user is required to create 

initial conditions for two-dimensional models by starting with a level water surface and 

for one-dimensional models by starting with a known water surface. TrimR2D does not 

require wet computational nodes to begin computations. However, because the 

simulations were started from a dry bed, a small time step was necessary. When the flow 

propagates over dry cells, the model can wet only the adjacent cell during a time step. For 

this case, the explicit terms in the discretized equation (eq. 13) dominate the first 

iterations of the solution and require time-step sizes similar to those used by fully explicit 

models. 

Scenario A: Roughness of 0.035 With Uniform Depth of 2 Meters 

The roughness for scenario A is typical of natural sand-and-gravel rivers. Watersurface profile simulations for outflow water depths of 2 m, 2.5 m, and 1.9 m are shown 

in figure 5. The outflow water depths of 2.5 m and 1.9 m represent a backwater (M1)

curve and a drawdown (M2) curve (Chow, 1959), respectively. Scenario A simulations 

are for 10,000 0.1-second time steps and reach flow equilibrium before the simulation

ends. Each simulated profile is for subcritical flow (F<1.0). 

 Model Performance 14
The computed uniform depth was 2.02 m for the three simulations and is 1 percent 

larger than the 2.00 m value computed by the Manning’s equation. The depth is in 

equilibrium for the backwater curve (M1 curve) at about 357 m from the outflow end and 

for the drawdown curve (M2 curve) at about 147 m from the outflow end. 

Figure 5. Simulated water-surface profiles for scenario A with a roughness of 0.035 and a 

uniform depth of 2 meters 

Davidian (1984) presents an empirical equation for estimating the distance needed 

for a downstream water depth that is 125 percent of uniform depth to converge to 

uniform depth for subcritical-flow conditions (M1 curve): 

o 2

n

LS =0.86-0.64F

Y , (23) 

where L is the required total reach length, So is the bed slope, Yn is the uniform depth, and 

F is the Froude number. Agreement between Davidian’s equation and the simulation 

results are good. For a water depth of 2.5 m at the outflow boundary (backwater 

simulation), the length of reach needed to converge within 3 percent of uniform depth is 

 Model Performance 15
210 m for equation (23), compared to 220 m for the simulation. The plot of flow depth 

versus channel distance in figure 6 illustrates the convergence of the simulation to 

uniform depth. The simulated backwater curve (M1 curve) slightly undershoots the depth 

before converging. The drawdown curve smoothly converges to the uniform depth. 

Figure 6. Simulated flow depth for scenario A with a roughness of 0.035 and a uniform depth of 

2.0 meters 

Scenario B: Roughness of 0.075 With Uniform Depth of 1Meter 

The roughness for scenario B is representative of streams with dense vegetation or 

with large boulders and steep water-surface gradients. Water-surface profiles for outflow 

water depths of 1.5 m and 0.75 m are shown in figure 7. The outflow water depths of 1.5 

m and 0.75 m represent a backwater (M1) curve and a drawdown (M2) curve, 

respectively. Subcritical flow occurs for each simulated profile (F<1). As in the previous 

scenario, scenario B simulations are for 10,000 0.1-second time steps and reach flow 

equilibrium before simulation ends. 

 Model Performance 16
Figure 7. Simulated water-surface profiles for scenario B with a roughness of 0.075 and a 

uniform depth of 1 meter 

The computed uniform depth for the two simulations is 1.06 m and is 6 percent 

larger than the 1.00-meter depth computed by the Manning’s equation. The depth is in 

equilibrium at about 870 m and 865 m from the outflow end for the backwater and 

drawdown curves, respectively (fig. 8). 

 Model Performance 17
Figure 8. Simulated flow depth for scenario B with a roughness of 0.075 and a uniform depth of 

1.0 meter 

Davidian (1984), in addition to the previous equation, presents an empirical equation 

for estimating the distance needed for a downstream water depth that is 75 percent of

uniform depth to converge to uniform depth for subcritical-flow conditions (M2 curve): 

2 0.57 0.79 o

n

LS F

Y = − . (24) 

For a water depth of 0.75 m at the outflow boundary, the length of reach needed to 

approach within 3 percent of uniform depth is 100 m for equation (24), compared to 110 

m for the simulation. However, the simulation does not converge to a flow depth that is 

within 3 percent of the uniform depth of 1.0 m. If the criterion is changed to be 3 percent 

of the depth to which the model converges, 1.06 m, the length of reach needed for 

convergence is about 594 m. The plot of flow depth with channel distance in figure 8

illustrates the convergence of the simulation to uniform depth. As in the previous 

simulation, the backwater curve undershoots the uniform depth before converging. For 

 Model Performance 18
the drawdown curve the simulation reaches a depth that is 98 percent of the uniform flow 

depth about 150 m from the outflow end. For distances greater than 450 m from the 

outflow end, the depth is greater than normal depth. 

Scenario C: Roughness of 0.012 With Uniform Depth of 3 Meters 

The third scenario is for a roughness of 0.012 that is typical for smooth concretelined channels. Supercritical flow occurs for each simulated profile. Water-surface

profiles for outflow water depths of 6.0 m and -0.27 m are shown in figure 9 for 

equilibrium conditions. Because the flow for this scenario is more advective than the two 

previous scenarios, a 0.04-second time step is necessary. Steady conditions were reached 

for the outflow depth that was less than the uniform depth within 10,000 time steps. For 

the outflow depth that is subcritical and greater than uniform depth, a hydraulic jump (an 

abrupt jump in water surface) initially occurs in the channel. During the simulation, the 

jump migrated to the outflow end. For this case, more time steps were required for the 

jump in the water surface to travel out of the channel reach. Steady conditions were 

reached within 20,000 time steps. 

Figure 9. Simulated water-surface profiles for scenario C with a roughness of 0.012 and a 

uniform depth of 3 meters 

 Model Performance 19
The simulated depth near the outflow boundary is 3.25 m and is 8 percent larger than 

the value computed by the Manning’s equation. Simulated depths are shown in figure 10. 

The computed depths near the outflow boundary may not be in equilibrium with gravity 

and friction. It is possible that a lower equilibrium flow depth near the outflow would be 

simulated if the channel were longer. 

Figure 10. Simulated flow depth for scenario C with a roughness of 0.012 and a uniform depth of 

3.0 meters 

Possible Sources of Simulation Error 

TrimR2D uses a constant flow depth within a cell. It is unlikely that this causes a 

significant error because the channel slope is constant. Error in the simulated normal 

depths did not decrease with increasing flow depth or correlate with roughness values. It 

may be that a higher mesh resolution could have reduced the error in simulated flow 

depth. Alternatively, the error in the simulated flow depths may be due to the numerical 

scheme used. No attempt was made to increase the model’s mesh resolution because of 

time constraints. 

 Model Performance 20
Laboratory Dam-Break Simulation 

Simulating a transient flow, such as a dam break, gives insight into how 

appropriately the advective terms are computed. Dam breaks typically generate steep, 

rapidly changing hydrographs that can challenge the computational stability of many 

numerical schemes. This section presents results from simulating a dam break. 

The dam break simulated is one of many laboratory experiments (Schmidgall and 

Strange, 1961) conducted by the U.S. Army Corps of Engineers (COE) during the 1950s 

to obtain stage and discharge data from dam breaks. The experiments were performed in 

a 121.92-meter-long laboratory flume with a 1.219-meter-wide rectangular section and a 

slope of 0.005. The dam was located at the midpoint of the flume length. The 

experiments simulated the sudden failure of a dam retaining a pool of water. 

The dam-break flow was for a constant inflow of 0.08767 m3

/s and an initial flow 

depth of 0.1707 m downstream of the dam. The base flow was passed beneath the dam. 

Roughness elements, consisting of 0.4445-centimeter (cm) aluminum angles, were 

installed on the flume floor perpendicular to the flow at 0.1524-m intervals. One leg of 

the aluminum angle was tacked to the flume floor, resulting in a roughness element 

height of 1.905 cm. Manning’s roughness coefficients estimated from uniform flow 

experiments were 0.04 at 0.21340-meter flow depth and 0.12 at 0.04572-meter flow 

depth. Total duration of the hydrograph following the dam failure was 240 seconds, with 

a peak discharge of 0.1500 m3

/s. Measured data are available for the experiment at 

several locations in the channel. Simulation results for two portions of the flume are 

presented: (1) a 38.1-meter reach downstream of the dam and (2) a reach that included 

both the impounded pool behind the dam and the reach downstream of the dam. 

Reach downstream of dam 

This simulation was for a portion of the laboratory flume downstream of the dam. 

The reach is located in the laboratory flume between 7.62 m and 45.72 m downstream of 

the dam. The computational mesh for the 38.1-meter-long reach has a square cell size of 

0.3048 m. The mesh is 6 cells wide by 125 cells long. The exit elevation of the flume in 

the computational mesh is 0.0 m. As for the uniform depth simulations, cell topography 

elevations on three sides of the computational mesh prevented flow through the mesh 

sides for all flows. 

The dam-break data for this reach are not very difficult to simulate accurately with 

typical one-dimensional, unsteady flow models. Fulford (1998) found that onedimensional models using a four-point implicit numerical method matched the dam-break 

hydrograph well. In the study, the roughness coefficients were adjusted for “dead” water 

behind the aluminum angles but were not altered to improve the fit of the onedimensional models to the data. 

The measured data for downstream water-surface elevation and upstream source 

discharge were used as boundary conditions. Discharge was divided equally among the 

four center nodes at the head of the flume. A converged solution computed by the 

TrimR2D model for a constant flow of 0.08767 m3

/s and a depth of 0.1707 m was used as 

the initial conditions for the simulation. A time-step size of 0.1 second was used. The 

 Model Performance 21
roughness value was determined by adjusting the roughness value until the simulated

flow depth approximately matched the flow depth measured for the initial discharge. The 

calibrated roughness was 0.055 and results in a 0.0035-meter underestimation of the flow 

depth at the inflow end of the reach downstream of the dam. 

Timing of the initial rise of the simulated peak (fig. 11) is good when compared with 

data measured at a location 21.34 m upstream of the outflow boundary (or 36.58 m from 

flume end). Simulated water-surface elevations (fig. 11) are generally lower than the 

measured values. During the recession of the peak, between 75 and 180 elapsed seconds, 

the simulated elevations match the measured elevations well. Simulated water-surface 

profiles at 1, 15, 30, and 55 elapsed seconds are plotted in figure 12. The plots of 

simulated water-surface profiles show that the rise in the water-surface elevation at 

21.34 m upstream of the outflow boundary is due mainly to the change in the downstream 

water elevation. 

Figure 11. Hydrographs for site located 21.335 meters upstream of the outflow boundary 

Modeled reach is downstream of dam 

 Model Performance 22
Figure 12. Simulated water-surface profiles of reach downstream of dam at 1, 15, 30, and 55 

seconds elapsed from start of simulation 

 Model Performance 23
Some of the error in the simulation may result from the use of point sources to inject 

the discharge into the model. Sources add flow volume to the model in the continuity 

equation and do not directly affect any of the terms in the momentum equations. Thus, 

any momentum (or the effects of velocity) from incoming flows is neglected by 

TrimR2D. The abrupt front of the wave from the dam break is likely smoothed by the

omission of the momentum forces associated with the inflow velocity. The omission of

inflow momentum effects may be causing the simulated peak to arrive later than the 

measured peak at the measurement site located 21.335 m upstream of the outflow 

boundary (fig. 11). Much of the fit in the hydrograph (fig. 11) at this location is the result 

of the propagation of the water-surface elevation boundary condition upstream from the 

outflow end. 

Pool behind dam and reach downstream of dam 

This simulation is for the portion of the laboratory flume that includes the 

impounded pool upstream of the dam and the reach downstream of the dam. This reach is 

especially challenging because of the abrupt change in water-surface elevation located at 

the face of the dam that occurs when the dam is removed. Many flow models, such as 

RMA2 (U.S. Army Corps of Engineers, 1996), Branch (Schaffranek and others, 1981), 

and FESWMS (Froehlich, 1989), are unable to compute for discontinuous or nearly 

discontinuous water-surface elevations or for changes in flow regimes. Typically, flow

models will have computational errors that will prevent the model from running over the 

specified period. 

Initially, the simulation was made for the portion of the flume beginning 49.28 m 

above the dam and ending 45.72 m below the dam. The initial-condition file used was a 

composite of two restart files computed with TrimR2D: a file with an outflow depth of 

0.1707 m and a file with an outflow depth that created an appropriate water-surface 

profile upstream of the dam. The resulting initial-conditon file represented the abrupt

jump in water-surface elevation at the dam face. Boundary conditions were a constant 

inflow of 0.08767 m3

/s and a time-varying stage measured in the laboratory at a location 

15.24 m from the flume outflow end. A solution was computed for all time steps. 

However, results from this simulation did not match measured values. 

The simulation therefore was repeated for the entire 121.92-meter flume length. The 

computational mesh used for the flume is a 6-by-410-cell mesh composed of 0.3048 msquare cells. A small change in elevation at the end of the flume was included in the 

topographic mesh to represent the height of an angle iron. Because the flow free falls out 

of the flume end in this simulation, it is likely that the flume edge condition affected the 

height of the backwater curve. The simulated flume has a roughness of 0.05. The selected 

roughness value was determined by matching the flow depth at a location 15.24 m 

upstream of the flume outflow end. 

Initial conditions were created by placing a weir at the location of the dam in the 

topographic mesh and by letting the simulation reach steady conditions. The weir height 

was determined by adjusting the height until the pool depths approximately matched 

measured values upstream of the dam. The constant inflow discharge of 0.08767 m3

/s was 

divided equally among the four center nodes at the head of the modeled reach, outside of 

the backwater zone produced by the dam. 

 Model Performance 24
The simulation of the dam break used the initial conditions with a topographic mesh 

that omitted the weir. Boundary conditions were the same as for the initial-condition 

simulation, constant inflow and free fall outflow. A time step of 0.02 second was used for 

the 240-second simulation. 

The simulated and measured water-surface elevation hydrographs for a location 

36.58 m upstream from the flume end are shown in figure 13. Timing of the initial rise of 

the simulated peak compares well with the measured data. The simulated water-surface 

elevations are biased higher than the measured values. This possibly is because the 

roughness calibration procedure focused only on the measured location at 15.24 m 

upstream of the flume end. 

Figure 13. Hydrographs for site located 36.58 meters upstream of flume end 

Modeled reach includes pool behind dam and reach downstream of dam 

 Model Performance 25
Figure 14. Measured and simulated water-depth profiles for pool behind dam and reach 

downstream of dam at various elapsed times. Dam is located 45.72 meters from outflow end 

 Model Performance 26
Figure 14. Measured and simulated water-depth profiles for pool behind dam and 

reach downstream of dam at various elapsed times. Dam is located 45.72 meters from 

outflow end—Continued

Simulated and measured flow depth profiles for the laboratory flume are shown in 

figure 14 for various simulation times. The computed depths match the measured depths 

very well for the simulation. Tables 2 and 3 present the percent errors between the 

simulated depths and the measured depths as a function of time and location, 

respectively. Generally, errors are less than +/-5 percent for the simulation. Based on 

location in the flume, average percent errors are largest near the dam location and 

smallest near the inflow end. Based on simulation time, average percent errors initially 

 Model Performance 27
declined to about 0.76 percent and then increased to about 4 percent over the duration of 

the simulation. 

Table 2. Mean and standard deviation for percent error of simulated flow 

depths as a function of time 

Time (seconds) 

0.5 1. 5. 10 20 30 50 70 100 120 180

Mean 1.44 2.71 0.76 0.76 0.90 2.06 2.27 2.99 3.01 4.23 3.90

Standard deviation 8.92 10.3 3.92 3.35 2.37 2.31 1.68 2.27 2.25 2.21 1.77

Table 3. Mean and standard deviation for percent error of simulated flow 

depths as a function of flume location 

Node location 

50 120 175 200 

Dam 201 204 210 220 250 280 330

Mean -2.94 2.37 2.97 9.07 6.06 3.19 1.50 1.00 -0.17 1.82 3.32

Standard 

deviation 1.21 1.42 2.22 12.6 5.70 3.66 3.19 3.40 4.04 2.73 0.79

Large-Scale River Simulations 

This section presents results for two large-scale river simulations—a reach in the 

Chattahoochee River, Georgia, and a reach in the Snoqualmie River, Washington. Both

simulations demonstrate the accuracy of the TrimR2D model applied to large river 

reaches. 

Chattahoochee River Simulation 

The Chattahoochee River simulation is for a 15-kilometer reach of the 

Chattahoochee River in northeastern Georgia. Data used in the simulation were collected 

by the USGS (Faye and Cherry, 1980) for the period of March 22-23, 1976. The data 

include flow, water-surface elevation, and cross sections. Cross sections are available at 

approximately half-mile intervals along the length of the river. The datum for the crosssection elevations is the National Geodetic Vertical Datum of 1929. Cross-section 

locations are referenced to river distance. Cross sections are not referenced to a common 

horizontal (northing and easting) datum. 

The simulated reach is from Littles Ferry Bridge, downstream of Buford Dam, to 

Georgia Highway 141. Flow data for the Chattahooche River at Littles Ferry Bridge and 

a tributary, Suwannee Creek, and stage data for Georgia Highway 141 were used as 

boundary conditions. Flow and stage data at Georgia Highway 120, located about 8.5 km 

downstream from Little Ferry Bridge, were used for comparison with the simulated 

values. Flow in the study reach is affected by the release of water that accompanies the 

 Model Performance 28
generation of hydroelectric power at Buford Dam. Flow released from the dam increased 

from 16.99 m3

/s on the morning of March 22 to 113.27 m3/s. This flow was maintained 

for about 20 hours before returning to the low flow. Subsequently, a peak flow of 226.53 

m3

/s, typical of a hydropower dam operations, was released on March 23.

The topographic mesh was interpolated from 19 cross sections using a square cell 

size of 4.8768 m. The cross sections were not referenced to a common horizontal datum 

in the report by Faye and Cherry (1980). This resulted in the mesh being a straightenedout representation of the actual reach and is similar to that used by one-dimensional flow 

models. Maps of the area could have been used to adjust the cross-section locations in the 

horizontal plane. However, it was not done because the horizontal adjustment would have 

required considerable effort and would have produced a large mesh of mainly non-active 

cells because the flow event was confined to the "low water" channel. The resulting 

computational mesh was 22 cells wide by 2,999 cells long. As in the previous 

simulations, flow out of the model domain occurred only at the downstream end. 

Initial roughness values were estimated from those presented by Faye and Cherry. 

Because roughness values vary spatially in the reach and TrimR2D does not permit 

spatial variation of roughness values, averaged values were estimated from data for cross 

sections located at the site of discharge measurements. Roughness values were adjusted 

to 0.035 for all flow depths. Various time-step sizes were tried, but a 60-second time step 

was the largest value that worked for this application. 

Initial conditions were computed using steady-state boundary conditions: 16.99 m3

/s 

inflow at Littles Ferry Bridge, 4.6723 m3

/s inflow at Suwannee Creek, and a downstream 

water-surface elevation of 268.34 m. The initial water-surface profile computed from the 

steady-state boundary conditions for the modeled reach is shown in figure 15. The watersurface elevation simulated at the Georgia highway 120 crossing was 0.234 m higher than 

the measured value. 

 Model Performance 29
Figure 15. Simulated water-surface profile for initial conditions for the modeled reach of the 

Chattahoochee River 

The arrival of the simulated peaks at the Georgia Highway 120 crossing matches the 

measured peak arrival well (fig. 16). The simulated water-surface elevation hydrograph 

for a 2-second time step matches the measured values well on the rises of the hydrograph. 

The computed hydrograph for a 60-second time step matches the measured hydrograph 

reasonably well. This hydrograph overestimates the peaks by approximately 0.05 m for 

the first peak and 0.1 m for the second peak. For both time-step sizes, the hydrographs for 

simulated water-surface elevations do not match as well on the recession of the peaks or 

at the lowest water-surface elevations. The simulated values for the lowest water-surface 

elevations were approximately 0.15 m larger than the measured values for both time-step 

sizes. 

It is possible that the relatively poorer match at the lower flow depths is due in part 

to how the topographic mesh was constructed. Errors in mean flow lengths for water 

parcels between ground features in adjacent cross sections are expected to be largest for 

the lower flow depths and least for the higher flow depths. In spite of this gross 

approximation of the channel topography, the model fit very well the timing and 

magnitude of the measured flow peaks. 

 Model Performance 30
 Figure 16. Hydrographs for March 22-23, 1976, at Georgia Highway 120 for measured and 

computed water-surface elevations for 2- and 60-second time-step sizes 

Figure 17 shows the simulated water-surface profile for a 2-second time step at 1, 

500, 1,000, 1,500, 1,750, and 2,200 elapsed minutes. The error caused by ignoring the 

momentum of the incoming flow is not evident in the water-surface profiles because the 

reach is much longer than the flume simulations. The proportion of the inflow momentum 

to the total momentum of the system is very small for this simulation. 

 Model Performance 31
Figure 17. Simulated water-surface profiles of Chattahoochee River at 30, 510, 990, 1,500, 1,740, 

and 2,190 minutes on March 22-23, 1976

 Model Performance 32
Figure 17. Simulated water-surface profiles of Chattahoochee River at 30, 510, 990, 1,500, 1,740, 

and 2,190 minutes on March 22-23, 1976—Continued 

Snoqualmie River Simulation 

The Snoqualmie River simulation is for a reach of the Snoqualmie River that is 

between Snoqualmie Falls and the Carnation Farm Road crossing of the Snoqualmie 

River. The simulated reach is about 27 km long and includes portions of the Raging and 

Tolt Rivers. The simulated reach covers an area that is 13,200 m by 8,240 m. The base 8-

meter-resolution topography grid used is derived from 2-meter-resolution Lidar elevation 

data. The resulting model mesh is a 515-by-825-node mesh with a 16-meter resolution. 

Flow is not simulated for nodes that are located above a topography elevation of 45 m. 

 Model Performance 33
The Snoqualmie simulations use discharges measured at USGS streamflow-gaging 

stations Raging River near Fall City (12145500), Tolt River near Carnation (12148500), 

and Snoqualmie River near Snoqualmie (12144500) and the measured water-surface 

elevation at Snoqualmie River near Carnation (12149000) as the boundary conditions. 

Time-step size used for the simulations is 1.5 seconds. This time-step size allowed the 

model to propagate the flow across dry areas and maintain model stability for the flow 

events tested. Snoqualmie River simulations ran at about one half real time for the 1.5- 

second time step. 

Two flow events were used to calibrate and verify the model, the flood events of 

November 24, 1986, and of December 3, 1975. Hourly inflow discharges for the 

calibration were computed from USGS stage hydrographs and stage-discharge ratings 

curves at the stream-gaging stations Snoqualmie near Snoqualmie (12144500), Raging 

River near Fall City (12145500), and Tolt River near Carnation (12148500). Measured 

data were not available for comparing the computed flooded areas with actual flood 

areas. 

The 1986 flood event was used to determine (or calibrate) a Manning's roughness for 

the modeled reach. The event had a single peak (fig. 18). The calibration simulation 

period is from when the flow was contained to the low-water channel at the most 

upstream stream gage 12144500 (November 23, 1986, 1:00 a.m.) until the peak stage 

occurred at the most downstream gage 12149000 (November 24, 1986, noon). The 

measured peak discharge was 1,645 m3

/s at the most downstream gage, 12149000, and 

1,617 m3

/s at the most upstream gage, 12144500. During the 1986 event, the Tolt River 

peaked before the flows at the most upstream and downstream gages 12144500 or 

12149000. 

 Model Performance 34
Figure 18. Hydrographs for the November 24, 1986, flood event for Snoqualmie River streamgaging stations 

Hydrographs begin on November 23, 1986, 1:00 a.m. and end on November 25, 1986, 11:00 p.m. 

(U.S. Geological Survey streamflow-gaging stations 12149000-Snoqualmie River near Carnation, 

12144500-Snoqualmie River near Snoqualmie, 12148500-Tolt River near Carnation, and 12145500-

Raging River near Fall) 

The inflow discharge hydrograph for the Tolt River includes only the flow for the 

South Fork of the river. Flow for the North Fork of the Tolt was not available for the 

simulation. Because the flow contributed by the North Fork is small relative to the total 

Snoqualmie River flow, this error was not considered significant. 

During calibration, the roughness value was adjusted until a reasonable match was 

obtained between the measured and simulated peak water-surface elevations. Manning's 

roughness values were varied from 0.035 to 0.120. Roughness values are constant in the 

Snoqualmie application regardless of flow depth. Simulated peak water-surface 

elevations for the 36-hour calibration hydrograph are shown in figure 19. Peak watersurface elevations were determined by selecting from the hourly model output the 

maximum elevation determined at each node over the simulation. Peak water-surface 

elevations were compared with peak elevations measured by the COE (U.S. Army Corps 

of Engineers, written commun., January 9, 2001) for the flood at five profile points 

located in the application reach (fig. 19). The accuracy of the peak water-surface 

 Model Performance 35
elevations measured by the COE is unknown. The best fit was for a Manning's roughness 

value of 0.12. Table 4 lists the measured and calibrated water-surface elevations for the 

five profile points. Errors for the calibrated peak water-surface elevations vary from 

 -0.53m to +0.67m. Peak water-surface elevations are overestimated at the upstream 

profile points and are underestimated at the downstream profile points. The average error 

for the COE profile points is -0.11m and the average absolute error is 0.41m. 

Figure 19. Simulated peak water-surface elevations on the Snoqualmie River for the November 

24, 1986, flood event for a Manning's roughness value of 0.12 

Peak elevations were measured by the U.S. Army Corps of Engineers at the five profile points 

 Model Performance 36
Table 4. Comparison of measured and simulated peak water-surface elevations 

for calibration of the Snoqualmie River model, using the flood event of November 

24, 1986, and a Manning’s roughness value of 0.12 

Peak water-surface elevation 

(meters) COE profile point 

number 

(See fig. 19 for 

location) Measured Simulated

Error 

(meters) 

14 23.04 22.51 -0.53

16 23.77 23.04 -0.73

17 24.07 24.03 -0.04

18 25.35 25.29 -0.06

22 33.95 34.62 +0.67

The 1975 flood event was used for verification. The verification simulation period is 

from December 1, 1975, at 1:00 a.m. (while the flow was still confined by the low-water 

channel) and until December 3, 1986, at 11:00 a.m. (when the peak arrived at the most 

downstream gaging station 12149000). The flood event has multiple peaks at the most 

downstream gage, 12144500 (fig. 20). The peak discharge was 1,475 m3

/s at the most 

downstream gage and was 1,467 m3

/s at the most upstream gage, 12144500. Peak flow at 

gage 12149000 for the 1975 flood event was sustained for a much longer period than for 

the 1986 flood event used for the calibration. The flow at the Tolt initially peaked at 

about the same time as at the Snoqualmie gage. 

 Model Performance 37
Figure 20. Hydrographs of the December 3, 1975, flood event for Snoqualmie River stream-gaging 

stations 

Hydrographs start on December 1, 1975 1:00 a.m. and end on December 6, 1:00 a.m. (U.S. 

Geological Survey streamflow-gaging stations 12149000-Snoqualmie River near Carnation, 

12144500-Snoqualmie River near Snoqualmie, 12148500-Tolt River near Carnation, and 12145500-

Raging River near Fall City.) 

Simulated peak water-surface elevations for the 59-hour verification hydrograph are 

shown in figure 21. Peak water-surface elevations were compared with those measured 

by the COE for the flood at five profile points located in the modeled reach. The accuracy 

of peak water-surface elevations measured by the COE is unknown. Table 5 lists the 

measured and simulated water-surface elevations for the five profile points for the 

verification event. Simulation errors range from -0.37 m to +0.68 m. The average error 

for the COE profile points is -0.09 m and the average absolute error is 0.36 m. Watersurface elevations are overestimated at the upstream profile point and are underestimated 

at the other profile points but are reasonable. The simulated water surface is flatter than

the measured water surface. This may be because the mesh cell size or the time-step size 

is too large or because of errors in the measured peak water-surface elevations. No effort 

was made to decrease cell size or time step because of the excessively increased 

computational times needed for the available desktop computer system. 

 Model Performance 38
Figure 21. Simulated peak water-surface elevations on the Snoqualmie River for the December 

3, 1975, flood event for a Manning's roughness value of 0.12 

Peak elevations were measured by the U.S. Army Corps of Engineers at the five profile points 

 Model Performance 39
Table 5. Comparison of measured and simulated peak water-surface elevations 

for verification of the Snoqualmie River model, using the flood event of December 

3, 1975, and a Manning’s roughness value of 0.12 

COE profile point Peak water-surface elevation (meters) 

number (See fig. 

21 for location) Measured Simulated

Error 

(meters) 

14 22.94 22.73 -0.22

16 23.46 23.09 -0.37

17 24.22 24.01 -0.21

18 25.60 25.26 -0.34

22 33.80 34.48 +0.68

SUMMARY AND CONCLUSIONS 

The numerical computer flow model Transient Inundation Model for Rivers – 2 

Dimensional (TrimR2D) can be used successfully to model realistic river systems, as 

demonstrated by the Chattahoochee and Snoqualmie River simulations. It can model 

these systems accurately with less user intervention than some other two-dimensional 

depth-averaged flow models because it can simulate changes in flow regimes and cell 

wetting and drying. TrimR2D is best used in reaches where roughness does not 

significantly vary areally. For river systems, the solution should not be used near the 

inflow boundary, and the inflow boundary should be located several channel widths 

upstream of the area of interest because of the omission of inflow momentum. 

Roughness values used in the uniform flow simulations do not precisely duplicate 

the results expected from Manning’s equation. The model tends to over estimate uniformflow depths. This behavior makes calibration of the roughness values used in an 

application important. Failure to calibrate the roughness values to measured data may 

result in over estimation of the water-surface elevations in the model. 

Roughness parameters cannot be varied horizontally in the model. This limits model 

application to locations where areal roughness changes are not significant. Changes in 

areal roughness can dramatically affect flow patterns and water elevations. Flood plains

where land use varies from agriculture to heavy woodland are an example where flood 

patterns are significantly affected by roughness changes. TrimR2D permits roughness 

values to change as a function of flow depth. In some situations, this may allow a higher 

roughness value to be used in the floodplain. In some situations, the current version of 

TrimR2D cannot accurately simulate flow where large changes in land use occur within 

the floodplain. 

To simulate horizontal inflows appropriately requires the careful use of source terms. 

The inability to impart flows with momentum can result in errors for flows near the 

inflow boundary and in reaches of short length. For typical river systems, the error in the 

computations caused by ignoring the inflow momentum is minimized for stream

 References 40
locations farther away from the inflow location. Alternately, inflow momentum can be 

imparted to the flow by using a pool (containing the source nodes) and weir at the inflow 

end of the mesh. 

The Chattahoochee River data are typical for moderate-sized rivers. A good match of 

the hydrograph was found for the Chattahoochee River data at a location that was several 

kilometers downstream of the inflow location. This finding demonstrates that the model 

results can be used successfully when the locations of interest are sufficiently far from the 

inflow location. The Snoqualmie simulations demonstrated the ability of TrimR2D to

propagate flow across a floodplain. Results for this simulation were reasonable. 

TrimR2D is capable of handling large local changes in water surface and velocity 

without diverging from the solution. This feature permits the model to compute steadystate or initial conditions without user intervention and to run for cases that are not 

possible for models such as RMA2 and FESWMS. The normal flow cases demonstrated 

that the model can easily be started from a dry bed and converge to a solution, even with 

the presence of sharp water-surface gradients in the flow. This is in contrast to flow 

models such as RMA2 and FESWMS. Modeling and serving on the Internet of real-time 

flood inundation requires a model that will run with very little, if any, user intervention.

TrimR2D requires less user intervention than other two-dimensional flow models. It is a 

robust model that is capable of giving accurate results for many types of river systems. 

REFERENCES 

Bates, P.D., and DeRoo, A.P.J., 2000, A simple raster-based model for flood 

inundation simulation: Journal of Hydrology, v. 236, p.54-77. 

Casulli, V., 1990, Semi-implicit finite difference methods for the two-dimensional 

shallow water equations: Journal of Computational Physics, v. 86, p.56-74. 

Casulli, V., and Cattani, E., 1994, Stability, accuracy and efficiency of a semiimplicit method for three-dimensional shallow water flow: Computers Math. 

Applications, v. 27, no. 4, p.99-112. 

Cheng, R.T., Casulli, V., and Gartner, J.W., 1993, Tidal, residual, intertidal mudflat 

(TRIM) model and its applications to San Francisco Bay, California: Estuarine, 

Coastal and Shelf Science, v.36, p.235-280. 

Chow, T.V., 1959, Open-channel hydraulics: New York, McGraw-Hill, p. 89-156. 

Conte, S.D., and deBoor, C., 1980, Elementary numerical analysis: an algorithmic 

approach: New York, McGraw-Hill, 

p. 356. 

Davidian, J., 1984, Computation of water-surface profiles in open channels: 

Techniques of Water-Resources Investigations of the United States Geological 

Survey, Book 3, Chapter A15, p. 17-18. 

Faye, R.E., and Cherry, R.N., 1980, Channel and dynamic flow characteristics of the 

Chattahoochee River, Buford Dam to Georgia Highway 141: U.S. Geological 

Survey Water-Supply Paper 2063, 65 p. 

 References 41
Fulford, J.M., 1998, Evaluation and comparison of four one-dimensional unsteady 

flow models: U.S. Geological Survey Water-Resources Investigations Report 

98-4242, 33p. 

Froehlich, D.C., 1989, Finite element surface-water modeling system: Twodimensional flow in a horizontal plan - Users manual: Federal Highway 

Administration RD-88-177, 285 p. 

Ostenaa, D.A., Levish, D.R., Klinger, R.E., and O’Connell, D.R.H., 1999, Phase 2 

paleohydrologic and geomorphic studies for the assessment of flood risk for the 

Idaho National Engineering and Environmental Laboratory, Idaho: Bureau of 

Reclamation Report 99-7, 117 p. 

Perry, C.A., March 2000, Significant floods in the United States during the 20th 

Century - USGS measures a century of floods: U.S. Geological Survey Fact 

Sheet 024-00, 4 p. 

Schaffranek, R.W., Baltzer, R.A., and Goldberg, D.E., 1981, A model for

simulations of flow in singular and interconnected channels: Techniques of 

Water-Resources Investigations of the United States Geological Survey, Book 7, 

Chapter C3, 110 p. 

Schmidgall, T., and Strange, J.N., November 1961, Floods resulting from suddenly 

breached dams conditions of high resistance: Miscellaneous Paper No. 2-374, 

Report 2, U.S. Army Corps of Engineers Waterways Experiment Station, 117 p. 

U.S. Army Corps of Engineers Waterways Experiment Station, 1996, Users guide to 

RMA2 version 4.3: U.S. Army Corps of Engineers, Waterways Experiment 

Station, Hydraulic Laboratory, 227 p. 

Walters, R.A., and Denlinger, R.P., 1999, Appendix C- Description of flood 

simulation models: Bureau of Reclamation Report 99-7, 12 p. 

 References 42
Appendix A. File Specifications for Model Input and 

Output Files 

The data needed for program execution and the results produced by the program are 

grouped into several files by type. The files containing data for program execution (input 

files) are designated by the file-name extensions .inp, .grd, .pts, and .win and the file 

name timebc.dat. They are free-format ASCII files and can be created with a text editor. 

Values must be separated by a space. Files containing program results (output files) are 

designated by the file-name extension .dat and tsdata.dat and the file name restart.dat.

The restart.dat file is a binary file (unformatted) and can be used to supply initial 

conditions. The other files with .dat extension and the file tsdata.dat are column-justified 

ASCII files. The file contents and variable definitions of each file type are described in 

this appendix. 

Program Control File 

The program control file has a file-name extension of .inp and contains indices used 

in the numerical solution, filenames for files containing input data, roughness parameters, 

initial conditions, and output file parameters. The TrimR2D program queries the user

with a file browser window for the name of this file when the program is executed. 

File contents 

Line number Variables 

1 - 2 comment lines describing simulation 

3 comment line written to *.dat ASCII output file 

4 FilenamesGRD 

5 FilenamesPTS

6 FilenamesWIN 

7 Outfilename 

8 comment line containing names for parameters on following line 

9 nskip, maxstp, nbegin, maxrex, nrstrt, noptype 

10 comment line containing names for parameters on following line 

11 epsi, qpsi, delt, theta, timebc 

12 comment line containing names for parameters on following line 

13 nlayer, wind-x, wind-y, zmin, etaobc 

14 comment line containing names for parameters on following line 

15 (n(k),k=1,nlayer+1) 

16 comment line containing names for parameters on following line 

17 (zn(k),k=1,nlayer+1)

18 comment line containing names for parameters on following line 

 Appendix A 43
19 nu-x, nu-y 

20 comment line containing names for parameters on following line 

21 ntsdata, ntsskip, tsscale 

22 comment line containing names for parameters on following line 

23 (irow(l),jcol(l),idf(l)), l=1,ntsdata 

Variable definitions 

FilenamesGRD – file name of grid file with .grd as the file extension 

FilenamesPTS – file name of the discharge file with .pts as the file extension 

FilenamesWIN – file name of the wind file with .win as the file extension

Outfilename – file name of the ASCII output file that contains program results 

nskip – integer number of time steps to skip writing program results to the file 

Outfilename, starting at time step nbegin 

maxstp – integer number of time steps executed in the simulation, elapsed time = maxstp 

x delt 

nbegin – integer time step (as an integer) at which to start output of program results 

maxrex – maximum integer number of iterations for the conjugant-gradient matrix solver 

nrstrt - integer flag for the restart file read/write option. The restart file can be used as 

initial conditions for a simulation. Allowable flags are: 

• 1 - an existing unformatted output file (restart.dat) is read and used as the initial 

conditions to initialize the dependent variables. 

• 2 - an unformatted file (restart.dat) is written at the end of the computer run. The 

file contains the results for the last executed time step. This file can be used as an 

initial condition file. 

• 3 - an existing unformatted file (restart.dat) is read and used as initial conditions 

and is written over with the results for the last executed time step. 

• 5 - an existing unformatted file (restart.dat) is read and used as initial conditions. 

A free-format file, rstij.dat, that contains the i,j indices and the values for the 

dependent variables for the active nodes is written. 

• 6 - the file rstij.dat is read and used to initialize the problem. 

noptype – integer flag for output type: 

• 1 – binary output files written at the requested time interval. The files have the 

name trimrout.nnn where nnn is the sequence number starting at 000. 

• 2 – The ASCII output file with .dat extension where the output at each requested 

time interval is written in a separate zone. This provides a simple means for 

animation of the results in the commerical graphics program TecplotTM by Amtec 

Engineering. 

 Appendix A 44
epsi – convergence tolerance for water-surface elevation, entered as a real number 

qpsi – not used in two-dimensional model 

delt – time-step size in seconds, entered as a real number 

theta – time-weighting factor (values range from 0.5 to 1.0), entered as a real number. A 

typical value is 0.6. 

timebc – integer flag for boundary-condition options: 

• 0 - boundary conditions are time-independent (steady). The water-surface 

elevation in the boundary cells are set to the value etaobc and the discharge (or 

flow) values are given in the *.pts input file. 

• ≠ 0 – file timebc.dat is read for the boundary condition parameters. 

nlayer – integer number of layers with separate roughness coefficients 

wind-x – initial wind velocity in x coordinate direction (east-west), in m/s, entered as a 

real number 

wind-y – initial wind velocity in y coordinate direction (north-south), in m/s, entered as a 

real number 

dzmin – minimum flow depth or depth at which a node is considered dry, in meters, 

entered as a real number 

etaobc – initial water-surface elevation, in meters, at the boundary cells (see fig. 4 in 

text), entered as a real number 

n(k) – Manning’s roughness coefficient for the kth layer, entered as a real number 

zn(k) – the max depth, in meters, for which the Manning’s roughness coefficient for the 

kth layer is applied, entered as a real number. The kth roughness coefficient is applied 

from the zn(k-1) to the zn(k) depths.

nu-x – horizontal eddy viscosity in the x-axis direction, in m2

/s, entered as a real number 

nu-y – horizontal eddy viscosity in the y-axis direction, in m2

/s, entered as a real number 

ntsdata – integer number of computational cells at which to write time series of 

simulation results 

ntsskip – integer number of time steps to skip before writing simulation results for the 

selected cells to the time series file 

tsscale – time scale for time-series output, entered as a real number. For instance, 

tsscale=60 gives a time scale in minutes. 

irow(l) – x-axis (east-west) cell position (elevation node; see fig. 3 in text) of the lth timeseries computational mesh location 

jcol(l) – y-axis (north-south) cell position (elevation node; see fig. 3 in text) of the lth 

time-series computational mesh location 

idf(l) – integer flag indicating the variable to output at the lth time-series location: 

• 1 – water-surface elevation 

• 2 - u component of velocity; (uj-1/2 + uj+1/2)/2 

 Appendix A 45
• 3 - v component of velocity; (vi-1/2 + vi+1/2)/2 

• 4 - speed; 

2 2

i-1/2 i+1/2 j-1/2 j+1/2 v v u u

2 2

⎛ ⎞ + ⎛ ⎞ +

⎜ ⎟ + ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

Model Mesh (Grid) File 

The model mesh (grid) file has a file-name extension of .grd and contains header 

data that describe the resolution (or cell dimension) and size (number of cells) of the 

regular rectangular mesh that is used by the model. It also contains the topography or 

ground elevations for the topography mesh. Ground elevations are specified for the user 

mesh at locations that are at half the spacing of the model grid. This eliminates the need 

to enter elevations for a staggered mesh. Elevations are referenced to a common vertical 

datum that is usually chosen to be the undisturbed water surface (or initial water-surface 

elevation) at the outlet. Elevations are positive in the downward direction and negative in 

the upward direction. The elevation values are read from upper left to lower right on a 

row-by-row basis. 

File contents 

Line number Variables 

1 nrow, ncol

2 dx, dy, x0, y0 

3 hland, hsea, latitude, zscale 

4 - last (depth(i,j),i=1,nrow),j=1,ncol) 

Variable definitions 

nrow - integer number (or index) of user-supplied ground-elevation locations in the xaxis direction 

ncol – integer number (or index) of user-supplied ground-elevation locations in the y-axis 

direction 

dx – model mesh spacing, in meters, for the x-axis direction, entered as a real number. 

The value is the spacing between elevation locations in the user-supplied mesh or onehalf the x-length of a mesh cell. 

dy – model mesh spacing, in meters, for the y-axis direction, entered as a real number. 

The value is the spacing between elevation locations in the user-supplied mesh. The 

computational cells are square, so dy equals dx. 

x0 – x-axis coordinate origin of the mesh, in meters, entered as a real number. This value 

is not used in the program. 

y0 – y-axis coordinate origin of the mesh, in meters, entered as a real number. This value 

is not used in the program. 

 Appendix A 46
hland – elevation of permanently dry land, in meters, entered as a real number. Cells with 

elevations higher than hland are removed from the calculations. This value is compared 

to depth(i,j) x zscale. 

hsea – elevation, in meters, used to detect water-surface-elevation boundary conditions

within the mesh. This value is compared to depth(i,j) x zscale. Set to a value higher than 

hland and enter as a real number. 

latitude – latitude, in degrees, used to compute the Coriolis force, entered as a real 

number. Because Coriolis forces are very small for riverine systems, set latitude = 0. 

zscale – elevation scale, entered as a real number. All depth(i,j) are multiplied by this 

value. Use a value of –1.0 to convert elevation data to coordinate system used by model. 

(depth(i,j),i=1,nrow),j=1,ncol) –ground elevation, in meters, listed starting at the upper 

left of the grid (1,1) and read row by row, ending at the lower right (nrow,ncol). Ground 

elevations are entered for all cells, including inactive cells. 

Time-Dependent Boundary Condition File 

The time-dependent boundary condition file, timebc.dat, is a free-format file. The 

file must be named timebc.dat and may be created using a text editor. TrimR2D can be 

compiled to use different or customized file-reading routines. Two versions of boundary 

file-reading routines are documented in this section. Both versions allow time-varying

discharges to be specified for the cell locations identified in the intial-discharge file (file 

extension .pts). The first version (time_bc.for) is for locations where the boundary cells 

water-surface elevation can be represented by either a constant or harmonic function. The 

second version (timebc_interp.for) is for locations where a time series of water-surface 

elevations is used for the boundary cells water-surface elevation. 

Harmonic Boundary Cell Water-Surface Elevation 

This version can be used to set water-surface elevations in the boundary cells, 

etaobch

, with the following equation: 

( ) numhar

h

m=1

etaobc = etaobc + ⎡amp(m)sin 2π (freq(m)ts-phase(m) ⎤ ∑ ⎣ ⎦

File contents 

Line number Variables 

1 numhar, numqbcs, numqsteps 

2-(numhar*2 +1) freq(m), amp(m), phase(m),m=1,numhar 

(numhar*2 +2) – last yr(ts), mon(ts), day(ts), min(ts),(ampq(s),s=1,numqbcs) 

Variable definitions 

numhar – number of harmonics in tidal water-surface-elevation boundary conditions

numqbcs – integer number of discharge-boundary locations, entered in the initial 

discharge file 

 Appendix A 47
numqsteps – number of discharges that are specified 

freq(m) – frequency of tidal constituent, in cycles per second, for the mth harmonic, 

entered as a real number 

amp(m) – amplitude of tidal constituent, in meters, for the mth harmonic, entered as a 

real number 

phase – phase angle for tidal constituent, in degrees, for the mth harmonic, entered as a 

real number 

yr(ts) – the year for the ts location in the time series, entered as an integer. The time 

series spans from the start to the end of the simulation. 

mon(ts) – the month for the ts location in the time series, entered as an integer 

day(ts) – the day for the ts location in the time series, entered as an integer 

min(ts) - the minute for the ts location in the time series, entered as integer 

(ampq(s),s=1,numbbcs) – the discharge values, in cubic meters per second, applied at 

time specified by yr(ts), mon(ts), day(ts), min(ts). The values are applied at the cell 

locations entered in the initial discharge file and are listed in the order specified in the 

initial discharge file. 

Time-Series Boundary Cell Water-Surface Elevation 

File contents 

Line number Variables 

1 netas, nqnodes, cunits, zconstant, timestep 

2-last wsel(t), cellq(s),s=1,nqnodes

Variable definitions 

netas – integer number of elevations. Always set to 1 

nqnodes – integer number of inflow cells 

cunits – multiplicative factor used to convert units in time-dependent boundary condition 

file to meters, entered as real number 

zconstant – constant, in meters, to add to Wsel to shift elevations in the vertical, entered 

as a real number 

timestep – time interval between entries in the file, in seconds, entered as a real number 

wsel(ts) – outflow water-surface elevation at time ts, entered as a real number 

(cellq(s),s=1,nqnodes) – discharge values, in cubic meters per second, at cell location s 

applied at time ts. The values are for the cell locations entered in the initial discharge file 

and are listed in the order specified in the initial discharge file. 

 Appendix A 48
Initial Discharge File 

The initial-discharge file, with a file name extension of .pts, is read by free format 

and may be created with a text editor. The file specifies the cell locations at which 

discharge boundary conditions are supplied and the initial discharge value in m3

/s. The 

number of lines in the file is equal to the number of cells at which discharge boundary 

conditions are supplied, plus one. 

File contents 

Line number Variables 

1 names of parameters on following line 

2 - last ics, jcs, cellQ(ics,jcs) 

Variable definitions 

 ics – east-west index in model coordinates of discharge node 

 jcs – north-south index in model coordinates of discharge node 

 cellQ(ics,jcs) – initial discharge, in cubic meters per second, for the ics, jcs cell 

Wind-Data File 

The wind-data file, with a file name extension of .win, is not usually used when river 

systems are simulated. It applies a wind velocity for a specified time over the entire 

computational grid. 

File contents 

Line number Variables 

1 Names for parameters on following line. 

2 - last hour, wx, wy

Variable definitions 

hour – the elapsed hour at which to apply wind velocity. The elapsed hour is computed 

from the start of a simulation. 

wx – wind velocity in the x-axis direction, in meters per second 

wy – wind velocity in the y-axis direction, in meters per second 

ASCII Output File 

The ASCII output file, with a file name extension of .dat, is written when the 

variable noptype in the program control file (file extension .inp) is set to 2. The time is

not written because a constant time interval is used for the file. Results are written at the 

time interval specified in the program control file (file extension .inp) by nskip. The 

output file name is entered in the program control file. The output file name must include 

the .dat extension. The file contains the results of the simulation in a format that can be 

 Appendix A 49
easily imported to a commerical graphics program, TecplotTM by Amtec Engineering. All 

values are reported at the computational cell centers in order to simplify visualization. 

File contents 

Line number Variables 

1 Title information listed on line 3 of the *.inp file 

1 ZONE I = inum,J= jnum 

2 - last i,j, h(i,j), eta(i,j,t), u(i,j,t), v(i,j,t), H(i,j,t)

Variable definitions 

inum – number of cells in x-axis direction of grid

jnum – number of cells in y-axis direction of grid

i – x-axis cell location (or index) 

j – y-axis cell location (or index) 

h(i,j) – averaged ground elevation, in meters, referenced to the vertical datum at cell i, j 

eta(i,j,t) – water-surface elevation, in meters, referenced to the vertical datum at cell i,j 

for time step t 

u(i,j,t)- averaged x-axis velocity at cell i,j for time step t 

v(i,j,t) – averaged y-axis velocity at cell i,j for time step t 

H(i,j,t) – water depth in meters at cell i,j for time step t 

Binary Output Files 

Multiple binary output files are written when the variable noptype in the program 

control file (file extension .inp) is set to 1. Binary files are written at the time interval

specified in the program control file (*.inp) by nskip. The files are named trimrout.nnn 

where nnn is the sequence number for the requested time interval. File contents are the 

same as the ASCII output file but omit the first 2 lines. 

Time-Series Output File 

The time-series file, tsdata.dat, contains the time series data requested in the program

control file (file extension .inp). Results are written at the time interval specified in the 

program control file (file extension .inp) by ntsskip. This file is always named tsdata.dat. 

File contents 

Line number Variables 

1-maxstp/nskip row, itsnodes(k,t),k=1,ntsdata

Variable definitions 

row – the line number 

 Appendix A 50
itsnodes(k,t) – the requested variable for time step t at the kth grid location specified in 

the program control file 

Binary Restart File 

The binary (or unformatted) restart file, restart.dat, can be either an input file or an 

output file. It contains data that can be used to supply initial conditions to the program.

Whether this file is written or used for initial conditions is specified in the program

control file (file extension .inp) by the nrstrt variable. 

ASCII Restart File 

The ASCII restart file, rstij.dat, can either be an input or output file. It contains data 

that can be used to supply initial conditions to the program and is named rstij.dat. 

Whether this file is written or used for initial conditions is specified in the program

control file (file extension .inp) by the nrstrt variable. 

File contents 

Line number Variables 

1-number of cells i,j, eta(i,j), u(i,j), v(i,j) 

Variable definitions 

i – east-west grid index 

j – north-south grid index

eta(i,j) – water-surface elevation, in meters, referenced to the vertical datum at cell i,j for 

the last time step 

u(i,j)- averaged east-west velocity at cell i,j for the last time step 

v(i,j) – averaged north-south velocity at cell i,j for the last time step 

 Appendix A 51
52