TRIGRS—A Fortran Program for Transient 

Rainfall Infiltration and Grid-Based Regional 

Slope-Stability Analysis, Version 2.0

By Rex L. Baum, William Z. Savage, and Jonathan W. Godt

Open-File Report 2008–1159

U.S. Department of the Interior 

U.S. Geological Survey
U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director 

U.S. Geological Survey, Reston, Virginia: 2008

Revised and reprinted: 2009

About USGS Products

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

About this Product

For more information concerning this program publication contact:

Jill McCarthy, USGS Central Region Geologic Hazards Team

Box 25046

Denver Federal Center

MS 966

Denver, CO 80225-0046.

(303) 273-8579

Or visit the Central Region Geologic Hazards Team Website at:

http://geohazards.cr.usgs.gov/

Publishing support provided by:

Denver Publishing Service Center, Denver, Colorado

Manuscript approved for publication, April 17, 2008

Suggested citation:

Baum, R.L., Savage, W.Z., and Godt, J.W., 2008, TRIGRS—A Fortran program for transient 

rainfall infiltration and grid-based regional slope-stability analysis, version 2.0: U.S. Geological 

Survey Open-File Report, 2008-1159, 75 p.

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

endorsement by the U.S. Government.

Although this report is in the public domain, permission must be secured from the individual 

copyright owners to reproduce any copyrighted material contained within this report.
iii

Contents 

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

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

System Requirements................................................................................................................2

Installation..................................................................................................................................2

Features........................................................................................................................................3

Changes from Previous Version................................................................................................. 3

Theoretical Basis of the Program ..................................................................................................4

Infiltration Models for Wet Initial Conditions ...............................................................................4

Infiltration Model for Unsaturated Initial Conditions ....................................................................7

Governing Equations...............................................................................................................7

Basic Solution for Constant Surface Flux................................................................................8

Early-Time Solutions for Constant Surface Flux......................................................................9

Later-Time Solutions for Variable Surface Flux.....................................................................10

Water Table Rise ..................................................................................................................11

Pressure Diffusion Below the Water Table...............................................................................15

Infinite-Depth Basal Boundary, Variable Applied Head .........................................................15

Finite-Depth Basal Boundary, Variable Applied Head...........................................................16

Early-Time Formula............................................................................................................16

Later-Time Formula............................................................................................................16

Infiltration, Runoff, and Flow Routing .......................................................................................17

Topographic Indexing............................................................................................................18

Slope Stability ..........................................................................................................................19

Limitations................................................................................................................................20

Program Operation......................................................................................................................22

TRIGRS Initialization File .........................................................................................................22

Input File Format ......................................................................................................................24

Output File Formats .................................................................................................................24

TopoIndex Initialization File......................................................................................................25

Troubleshooting .......................................................................................................................25

Comparison of Computed Results with Published Analyses.......................................................25

Acknowledgments .......................................................................................................................26

References Cited ........................................................................................................................26

Appendix 1. Listing of Files Included in This Distribution ...........................................................59

Appendix 2. Tutorial ...................................................................................................................62

Using TopoIndex to Prepare Data for Runoff Routing..............................................................62

Preparing the TRIGRS Initialization File for an Analysis ..........................................................64

Using TRIGRS to Analyze Slope Stability during a Simple Storm............................................67
iv

Figures 

1. Sketch showing local xyz coordinate system. ........................................................................28

2. Sketch showing the definition of the vertical coordinate, += zxZ cossin  , used to compute 

the depth at an arbitrary location. .....................................................................................29

3. Schematic cross sections showing boundary conditions assumed in deriving the two 

infiltration formulas (equations 1 and 2)............................................................................30

4. Hypothetical rainfall history of varying duration and intensity.................................................31

5. Example of transient water-table rise and pressure-head changes in response to precipitation 

at a cell where steady upward flow occurs but the water table is initially beneath the 

ground surface. ................................................................................................................32

6. Shallow ground-water conditions in hillside soils....................................................................33

7. Vertical soil water-content profiles from the ground surface down to an arbitrary depth, Z, 

showing 

the water-table rise as a result of infiltration. ....................................................................34

8. Plot of example input and output of unsaturated zone infiltration model. ...............................35

9. Graph showing example of Gardner (1958) models for soil water characteristics. ................36

10. TRIGRS grid cell numbering scheme...................................................................................37

11. Diagram of an elevation grid file, comparing values of (A) elevation and (B) topographic 

index.................................................................................................................................38

12. Diagram comparing numeric codes for D8 flow directions used by the utility program 

TopoIndex 

with those used by ESRI Arc/Info. ....................................................................................39

13. Example weighting factors for runoff distribution from an arbitrary cell. ...............................40

14. Vertical spacing for computing infiltration and factor of safety..............................................41

2-1. Screen image of the file "dem.asc" with grid dimensions highlighted.................................. 70

2-2. Screen image of the file "TopoIndexLog.txt" with array length parameters highlighted.......71

2-3. Sample analysis of data from the folder “TRIGRS/data/tutorial”. ........................................72

Tables

1. Weighting factors for runoff distribution..................................................................................42

2. Sample listing of the initialization file "tr_in.txt".......................................................................43

3. Definitions of input variables, excluding file path names, found in the initialization file and 

related input files. .............................................................................................................45

4. Output files generated by TRIGRS.........................................................................................48

5. Example listings of a grid file for a 10-row by 10-column grid. ..............................................49

6. Partial listing of the file "TIdscelList_tutorial.txt", which contains the cell numbers of 

downslope receptor cells. ................................................................................................. 51

7. Partial listing of the file " TIwfactorList_tutorial.txt ", which contains the weighting factors of 

downslope cells. ...............................................................................................................52

8. Partial listing of the file " TIcelindxList_tutorial.txt ", which contains the cell index for a grid of 

100 cells taken from the tutorial data set..........................................................................53

9. Partial listing of the TRIGRS output file " TRlist_z_p_fs_XXXXXXXX.txt ", depth profile of 

pressure head and factor of safety for each cell. ..............................................................54

10. Example listing of entries in the initialization file for the program TopoIndex, "tpx_in.txt".....55

11. Names of TopoIndex output files..........................................................................................56

12. Listing of GridMatch initialization file, “gm_in.txt.” ................................................................57
v

13. Common problems and solutions for users of the TRIGRS and TopoIndex programs.........58

2-1. Contents of the initialization file “tpx_in.txt” for tutorial problem, listed by line number. ......73

2-2. Contents of initialization file “tr_in.txt” for tutorial problem, listed by line number. ...............74

Conversion Factors 

Inch/Pound to SI

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

Area

square foot (ft2) 0.09290 square meter (m2) 

Flow rate

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

Mass

pound, avoirdupois (lb) 0.4536 kilogram (kg) 

Pressure

pound per square foot (lb/ft2) 0.04788 kilopascal (kPa) 

pound per square inch (lb/in2) 6.895 kilopascal (kPa) 

Density

pound per cubic foot (lb/ft3) 16.02 kilogram per cubic meter 

(kg/m3

) 

Hydraulic conductivity

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

SI to Inch/Pound

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

Area

square meter (m2) 10.76 square foot (ft2) 

Volume

cubic meter (m3) 35.31 cubic foot (ft3) 

cubic meter (m3) 1.308 cubic yard (yd3) 

Flow rate

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

Mass

kilogram (kg) 2.205 pound avoirdupois (lb)
vi

Pressure

kilopascal (kPa) 20.88 pound per square foot (lb/ft2) 

Density

kilogram per cubic meter 

(kg/m3

) 

0.06242 pound per cubic foot (lb/ft3

) 

Hydraulic conductivity

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

TRIGRS—A Fortran Program for Transient 

Rainfall Infiltration and Grid-Based Regional 

Slope-Stability Analysis, Version 2.0

By Rex L. Baum, William Z. Savage, and Jonathan W. Godt

Abstract

The Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Model

(TRIGRS) is a Fortran program designed for modeling the timing and distribution of shallow, 

rainfall-induced landslides. The program computes transient pore-pressure changes, and attendant 

changes in the factor of safety, due to rainfall infiltration. The program models rainfall infiltration, 

resulting from storms that have durations ranging from hours to a few days, using analytical 

solutions for partial differential equations that represent one-dimensional, vertical flow in isotropic, 

homogeneous materials for either saturated or unsaturated conditions. Use of step-function series 

allows the program to represent variable rainfall input, and a simple runoff routing model allows 

the user to divert excess water from impervious areas onto more permeable downslope areas. The 

TRIGRS program uses a simple infinite-slope model to compute factor of safety on a cell-by-cell 

basis. An approximate formula for effective stress in unsaturated materials aids computation of the 

factor of safety in unsaturated soils. Horizontal heterogeneity is accounted for by allowing material 

properties, rainfall, and other input values to vary from cell to cell. This command-line program is 

used in conjunction with geographic information system (GIS) software to prepare input grids and 

visualize model results.

Introduction 

The Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Model

(TRIGRS) is a Fortran program for computing transient pore-pressure changes, and attendant 

changes in the factor of safety, due to rainfall infiltration. The original version (Baum and others, 

2002) was based on the method outlined by Iverson (2000), with implementation of complex storm 

histories, an impermeable basal boundary at finite depth, and a simple runoff routing scheme. In 

version 2, we have retained the features of the original version. We have, however, expanded the 

model to address infiltration into a partially unsaturated surface layer above the water table by

using an analytical solution of the Richards equation for vertical infiltration (given by Srivastava 

and Yeh, 1991). The analytical transient flow model is one-dimensional and represents the vertical 

infiltration in homogeneous isotropic materials that occurs during a storm. To improve usability, 

we have implemented property zones to simplify model input and have made numerous small 

changes. The program operates on a gridded elevation model of a map area and accepts input from 

a series of ASCII text files. Infiltration, hydraulic properties, and slope stability input parameters 

are allowed to vary over the grid area thus making it possible to analyze complex storm sequences 

over geologically complex terrain. The optional routing scheme achieves mass balance between 

rainfall input, infiltration, and runoff over the entire grid by allowing excess water to flow to 
2

downslope cells that are receiving less direct precipitation than they are able to absorb. The 

program saves output to a series of text files that can be imported to GIS software for display or 

further analysis. 

Analyses using TRIGRS are applicable to areas that are prone to shallow precipitationinduced landslides and that satisfy other model assumptions reasonably well. These assumptions 

include a well-documented initial water table and steady background flux, and relatively isotropic, 

homogeneous hydrologic properties. The soil-water characteristic curves are approximated by an 

exponential function (Gardner, 1958). However, we have not evaluated the magnitude of such 

possible errors. Model results are very sensitive to the initial conditions, particularly the steady 

component of the flow field and initial depth of the water table. Consequently, the model may 

produce questionable results where the initial water table depth is poorly estimated. 

The TRIGRS program runs from the command line or in a simple input/output window with 

relatively little user interaction. The user controls an analysis by means of an initialization file that 

contains the names of all other input and output files as well as other parameters needed to run the 

program. The following discussion will briefly describe system requirements for running the 

program, installation, and features and limitations of the program, and will provide a detailed 

description of the initialization file. An included tutorial helps familiarize the user with program 

features and operation. 

System Requirements

The TRIGRS program will run on computing platforms that support Fortran 77 and Fortran 

90/95. The main program is written in standard Fortran 90 and uses Fortran 90 and Fortran 77 

subroutines.1

 This distribution includes the source code and executable2

Installation

files that will run on 

Windows and Mac OS computers. The Windows version requires Windows 2000/XP and a 

Pentium-class processor. The Mac OS version requires Mac OS X (version 10.2 or later) and a 

PowerPC G3 (or higher) processor; TRIGRS runs from the command line, inside the Terminal 

application. The amount of RAM required depends on the problem size and the operating system. 

Analysis of a 24-period storm over a typical 30-m DEM for a USGS 7-1/2 minute quadrangle, 

which has about 175,000 cells, requires about 44 MB of dedicated RAM (in addition to what is 

required by the operating system). Analysis of the same storm over a 10-m DEM (about 1,575,000 

cells) requires about 400 MB of dedicated RAM. The executable programs, source code, and 

tutorial/test data sets occupy about 2.2 MB of disk space. Disk space for data storage requirements 

depend on the problem size and output options, but can become large (hundreds to thousands of 

MB). 

The software can be downloaded from the U.S. Geological Survey publications web server 

at http://pubs.usgs.gov/of/2008/1159/. Check http://landslides.usgs.gov/research/software.php for 

updates and bug fixes to the code. After decompressing the archive file (OF20081159_W.zip, for 

Windows or OF20081159_M.zip, for Mac OS X), copy the folder TRIGRS to any desired folder on 

your Windows PC, or Macintosh computer. Appendix 1 contains a listing of files included in this 

distribution. Appendix 2 gives a brief tutorial to acquaint the user with program operation.

 

1 Some of the subroutines make use of widely implemented, but non-standard, Fortran 77 extensions for string 

manipulation (“adjustl()”) and type declaration (“implicit none”).

2

Portions of this program include material copyrighted © by Absoft Corporation 1988–2008. 
3

Features 

As mentioned previously, TRIGRS implements analytical solutions similar to those 

described by Iverson (2000) to compute both transient pore pressure due to vertical infiltration and 

attendant changes in factor of safety in a digital landscape. We have extended Iverson's analysis by 

adding solutions for an impermeable basal boundary at a finite depth, an optional unsaturated zone 

above the water table, and a simple runoff-routing scheme to disperse excess water from cells 

where the application rate (rainfall intensity) exceeds the infiltrability (or infiltration capacity, 

Hillel, 1982). The program allows the following input parameters to vary from cell to cell 

throughout the model: precipitation intensity, slope, soil depth, initial water-table depth, saturated 

vertical hydraulic conductivity, hydraulic diffusivity, a three-parameter soil-water characteristic 

curve, cohesion for effective stress, angle of internal friction for effective stress, and total unit 

weight of soil. The precipitation rate may vary through time and may include water from other 

sources, if desired. The user can observe the changes in the pore pressures and factors of safety as 

the storm progresses (after program execution is completed). The program reads ASCII grid files 

in standard formats that can be output by recent versions of ArcGIS (version 9.2 and earlier) and 

GRASS GIS (version 6.3 and earlier). Program output is in the same ASCII grid format as the 

input except for an optional file that contains a listing of pore pressure and factor of safety at 

various depths for each grid cell. The TRIGRS program restricts the factor of safety to a maximum 

of 10 to facilitate plotting and eliminate division by zero in computing the factor of safety on flat 

slopes. To speed computation, the user can specify a minimum slope angle for computations; cells 

that are flatter than the minimum slope are assigned a factor of safety of 11 to distinguish them 

from the others.

The main program was written in Fortran 90 to take advantage of its ability to allocate 

arrays dynamically. Consequently, TRIGRS can be used to analyze grids of various sizes and 

storms having an arbitrary number of periods of different rainfall intensity without the need to 

recompile the program. The upper limit of the problem size that can be analyzed depends on 

compiler options, the operating system, and hardware being used. The interested user can find the 

source code in the TRIGRS folder. 

Changes from Previous Version

The addition of the previously mentioned new features, primarily the unsaturated infiltration 

model and the introduction of property zones and a property-zone grid in place of individual grids 

for each physical property parameter, required major restructuring of the initialization file. Users 

of the previous version should consult the section on "Program Operation" for information about 

the new organization of the initialization file, new options, and related topics. The Mac OS version 

now runs from the command line, inside the Terminal application, and all input and output files 

now use UNIX line endings. 

In the original version of TRIGRS each run analyzed the data for a single specified elapsed 

time in the storm sequence. To observe how pore pressures and factors of safety change during a 

storm, it was necessary to run TRIGRS for different values of elapsed time, t . Addition of an 

unsaturated zone above the water table necessitated restructuring of TRIGRS to solve for multiple 

elapsed times during a single program run. This restructuring now makes it possible to compute 

pore pressures and factors of safety for multiple elapsed times during a storm sequence; results are 

saved to a series of grid files for different elapsed times. 
4

Theoretical Basis of the Program

TRIGRS combines models for infiltration and subsurface flow of storm water, routing of 

runoff, and slope stability to calculate the effects of storms on the stability of slopes over large 

areas. The following is a brief description of the basic models and formulas used by TRIGRS to 

represent these processes. Formulas used in TRIGRS refer to the coordinate system depicted in 

figures 1 and 2. 

Infiltration Models for Wet Initial Conditions

The infiltration models in TRIGRS for wet initial conditions are based on Iverson’s (2000) 

linearized solution of the Richards equation and our own extensions to that solution (Baum and 

others, 2002; Savage and others, 2003, 2004). Iverson’s solution consists of a steady component 

and a transient component. The steady seepage component depends on the initial depth to the 

water table and a steady infiltration rate. The steady infiltration rate, the saturated hydraulic 

conductivity, and the slope angle determine the steady (initial) flow direction (within the x-z plane). 

The steady vertical hydraulic gradient is a function of the slope, the initial (steady or background) 

infiltration rate, and the saturated hydraulic conductivity. Model results are very sensitive to the 

steady-seepage initial condition. Consequently, accurate initial conditions are probably the most 

significant input for the model. Obtaining accurate initial conditions requires a significant number 

of field observations possibly supplemented by steady-flow modeling to extrapolate between 

observations. In the absence of accurate initial conditions, use of TRIGRS is limited to modeling 

hypothetical scenarios. 

The transient component assumes one-dimensional, vertical downward flow, with a simple 

specified time-varying flux of fixed duration and intensity at the ground surface and a zero flux 

condition for times greater than the starting time at an infinitely deep basal boundary (fig. 3). 

These infiltration models apply to saturated or tension-saturated initial conditions, so that flow is in 

the linear range for Darcy’s law and the hydraulic diffusivity is approximately constant. We have 

used a series of Heaviside step functions to implement Iverson’s (2000) suggested summation of 

his original solution for rainfall of constant intensity to represent a general time-varying sequence 

of surface fluxes of variable intensities and durations. This generalized solution used in TRIGRS is 

given: 

( )

1

2

1 1

1 2

1

1

2

11 1 1

1 2

1 1

,( )

2 ( )[ ( )] ierfc

2[ ( )]

2 ( )[ ( )] ierfc

2[ ( )]

N

nZ

n n

n S

n

N

nZ

n n

n S

n

Zt Z d

I Z Ht t Dt t

K Dt t

I Z Ht t Dt t

K Dt t

 

=

+ +

=

+

= −

        + −−           −  

        − −−           −  

∑

∑

(1A) 

In equation 1A, 

 is the ground-water pressure head; 

t is time; 
5

= zZ cos , where Z is the vertical coordinate direction (positive downward) and depth 

below the ground surface, z is the slope-normal coordinate direction (also 

positive downward), and  is the slope angle (figs. 1 and 2); 

d is the steady-state depth of the water table measured in the vertical direction; 

2 cos ( / ) ZLT S   = − I K ; 

KS is the saturated hydraulic conductivity in the Z direction; 

ZLT I is the steady (initial) surface flux;

InZ is the surface flux of a given intensity for the n

th time interval; 

2

1 0 D D= / cos  , where D0 is the saturated hydraulic diffusivity ( SKD SS / 0 = , where KS is 

the saturated hydraulic conductivity and SS is the specific storage); 

N is the total number of time intervals; and 

H t − tn ( ) is the Heaviside step function and nt is the time at the th n time interval in the 

rainfall infiltration sequence. 

The function ierfc is of the form 

( ) ( ) ( ) 1 2 ierfc exp erfc   

 = −− , (1B) 

where erfc( )  is the complementary error function. The first term in equation 1A represents the 

steady part of the solution and the remaining terms represent the transient part. Figure 4 illustrates 

the kind of storm history that equation 1A is capable of modeling. 

As an alternative to the solution with an infinitely deep basal boundary, we have added and 

implemented in TRIGRS a solution for pore pressure in the case of an impermeable basal boundary 

at a finite depth LZ d (fig. 3). This solution is given: 

1

2 1 1 1

1 1 2 2 1 1

1

2 11 1

1

( ,) ( )

(2 1) ( ) (2 1) ( ) 2 ( )[ ( )] ierfc ierfc

2[ ( )] 2[ ( )]

(2 1) 2 ( )[ ( )] ierfc

N

nZ LZ LZ LZ LZ

n n

n m S

n n

N

nZ LZ

n n

n S

Zt Z d

I md d Z md d Z Ht t Dt t

K Dt t Dt t

I m d Ht t Dt t

K

 

∞

= =

+ +

=

= −

         − −− − +−    + −−   +          − −     

− − − −−

∑ ∑

∑ 1 1

1 2 2 1 1 1 1

( ) (2 1) ( ) ierfc

2[ ( )] 2[ ( )]

LZ LZ LZ

m

n n

d Z md d Z

Dt t Dt t

∞

=

+ +

         − − +−      +          − −     

∑

 (2) 

Except for the addition of dLZ, the depth of the impermeable basal boundary measured in the Zdirection, the notation in equation 2 is the same as that in equation 1. The infinite series in 

equations 1 and 2 converge rapidly so that in computation, only the first few terms of the series are 

needed. The version of TRIGRS described in this report tests these series for convergence so that 

the minimum number of terms needed to achieve accurate, consistent results is used at each grid 

cell.

Equations 1 and 2 represent two different subsurface conditions (fig. 3). Equation 1 applies 

where hydraulic properties are uniform with depth and equation 2 applies where a well-defined 
6

decrease in hydraulic conductivity exists at a finite depth. For example, equation 2 would apply 

where loose relatively permeable slope deposits overlie relatively impermeable bedrock. We 

considered allowing the user to specify a flux at the basal boundary; however, this would add 

another unknown to the problem without adding significantly to the program’s analytical 

capabilities. Note that only one solution (either equation 1 or 2) is applied to all cells in an 

analysis. 

Following Iverson (2000) we impose an additional physical limitation that pressure head 

under downward gravity-driven flow cannot exceed that which would result from having the water 

table at the ground surface, that is: 

 ),( ≤ ZtZ  (3A) 

Thus, TRIGRS automatically reduces pressure head to Z anywhere that equation 1 or 2 has 

computed a value that is higher than Z . Setting the background flux, IZLT , to zero results in a 

limiting pore pressure distribution equal to what occurs during slope-parallel flow. As a result of 

user requests, we have somewhat relaxed this requirement by allowing the user the option to 

specify a value of  that represents either hydrostatic ( =1) or slope-parallel ( = cos2 ) flow 

conditions rather than using the more general definition used by Iverson (2000). 

TRIGRS applies a modified form of equation 3A at cells where the initial flow direction is 

upward. If the steady water table is at the ground surface, then exfiltration occurs and any 

precipitation runs off along with the exfiltrating water. If the initial water table is below the ground 

surface and flow is upward, then TRIGRS allows precipitation to infiltrate until it raises the water 

table to the ground surface. This adjustment is accomplished by adding the transient and steady 

pressure-head components above the water table as at other cells by applying either equation 1 or 2 

and by maintaining the steady upward hydraulic gradient below the water table as indicated in 

figure 5. The following formula applies below the water table in such cells.

 −= dZtZ t)(),(  ; Z > dt (3B) 

In equation 3B, dt is the depth to the transient water table.

Both equations 1 and 2 are generalizations of analogous heat-flow problems (Carslaw and 

Jaeger, 1959). Equations 1 and 2 are valid for saturated and nearly saturated soils for which the 

infiltrability and hydraulic conductivity approximately equal the saturated hydraulic conductivity. 

The solutions are strictly valid only where transient infiltration is vertically downward, but, in 

application, computed pore pressures will likely be accurate enough if the transient lateral flow 

component is relatively small. However, because the steady-flow component may have any 

direction consistent with possible values of  , superposition of the transient and steady 

components can result in flow in an arbitrary direction within the x-z plane. As a result of 

approximations made in linearizing the Richards equation, the formulas are valid for homogeneous 

soils and for depths that are shallow compared to the square root of the upslope contributing area 

(Iverson, 2000). 

The formulas for both solutions of the infiltration model (equations 1 and 2) include the 

complementary error function. The program, TRIGRS, uses an approximation of the error function 

that was developed by W.J. Cody of Argonne National Laboratory (Cody, 1969). The original 
7

Fortran source code, also by W.J. Cody, is archived online at http://www.netlib.org/specfun/erf

(accessed 10/17/2006). 

Infiltration Model for Unsaturated Initial Conditions

To make the TRIGRS model applicable to a wider range of initial conditions, we have 

added an analytic solution for unsaturated flow as an option for estimating infiltration at the ground 

surface. This option treats the soil as a two-layer system consisting of a saturated zone with a 

capillary fringe above the water table overlain by an unsaturated zone that extends to the ground 

surface (fig. 6). The unsaturated zone absorbs part of the water that infiltrates the ground surface 

and the remaining water passes through the unsaturated zone and accumulates at the base of the 

unsaturated zone above the initial water table (fig. 7A and 7B). The unsaturated zone acts like a 

filter that smoothes and delays the surface infiltration signal at depth (fig. 8).

This model uses four parameters (r , s , and KS

 ) to approximate the 

soil-water characteristic curve for wetting of the unsat

, as defined below, and

urated soil (Gardner, 1958) and approximates 

the infiltration process as one-dimensional, vertical flow (Srivastava and Yeh, 1991; Savage and 

others, 2004). The infiltrating water accumulates at the base of the unsaturated zone and thus raises 

the water table. The load of water accumulating above the original water table propagates 

downward as a diffusive pressure wave or waves, thus raising pore pressure at depth. For relatively 

thin saturated zones, like those common in hillsides, the pressure waves propagate to the basal 

boundary very rapidly. A separate analytical solution models propagation of these pressure waves 

through the saturated zone. In the following paragraphs we first describe formulas used for 

computing flux, hydraulic conductivity, and pore pressure distribution in the unsaturated zone. 

Next, we describe the procedure for estimating water-table rise from the unsaturated zone 

computations and then describe how pore pressure changes are computed below the original water 

table. 

Governing Equations

A one-dimensional form of the Richards equation describes infiltration at the ground 

surface and vertical flow through the unsaturated zone (Freeze and Cherry, 1979). We use the 

coordinate transformation described by Iverson (2000) to account for the effects of a sloping 

ground surface: 

( ) 2

1 1

cos

K

t Z Z

           = −         (4A) 

As shown by Srivastava and Yeh (1991), the exponential hydraulic parameter model of 

Gardner (1958) can be used to linearize equation 4A and provide an analytic solution for transient 

infiltration through the unsaturated zone. Gardner’s (1958) unsaturated conductivity model is 

generally considered appropriate for coarse-grained soils. In this model, dependence of hydraulic 

conductivity, K( )  , and fluid content,  , on pressure head in the Richards equation is given by the 

following formulas:

( ) ( )  = KK S exp  * (4B) 
8

( )  *)exp( −+= rsr (4C) 

In equations 4B and 4C,  is the pressure head, * = −0 , where 0 is a constant defined below, 

KS is the saturated hydraulic conductivity, K( )  is the hydraulic conductivity function,  is the 

volumetric water content, r is the residual water content, and s is the water content at saturation. 

The parameter  is obtained by fitting equation 4B to a characteristic curve for soil such as that 

shown in figure 9. As we recall from the discussion of figure 6, the hydraulic conductivity is 

constant, independent of pressure head, and equal to KS within the capillary fringe and below the 

water table. Also as shown by Gardner (1958), 1  represents the (vertical) height of the capillary 

fringe above the water table and, thus, at the top of the capillary fringe the pressure head is

0 −= 1  . By substituting 0 = −1  into the formula for  *, equations 4B and 4C reduce to 

K( )  = KS and =  S at the top of the capillary fringe, and we obtain the soil-water characteristic 

curve shown by the solid line in figure 9. 

Substitution of equations 4B and 4C into the one-dimensional version of the Richards 

equation 4A, leads to a linear partial differential equation with the solutions given below 

(Srivastava and Yeh, 1991). Note that for the coordinate transformation for the sloping ground 

surface we use the following:  2

1 = cos . 

Basic Solution for Constant Surface Flux

The flux at the base of the unsaturated zone, obtained by solving the linearized Richards 

equation, is given: 

1

1 2

1 1 2

1

4( )exp exp 2 4

( ,) sin( ) exp

1 2

2

u

Z Z ZLT

u m mu

m

m u

m u

d t I II D

qd t d D t d d













∞

=

      − −   −        

=   Λ Λ   ⋅   −Λ  

+ +Λ  

∑ (5A) 

in which 

1

( )

S

s r

K D



  = − (5B) 

is the soil-water diffusivity (Freeze and Cherry, 1979), D t  is equivalent to the non-dimensional 

time used by Srivastava and Yeh (1991), du is the vertical depth to the top of the capillary fringe 

(fig. 2), 1 u  d is equivalent to the non-dimensional depth used by Srivastava and Yeh (1991), and 

the values of Λm are the roots of the characteristic equation shown below. All other notation is 

identical to that used previously. 

The hydraulic conductivity in the unsaturated zone varies with time and depth according to 

the following formulas. The initial condition, K Z( ) ,0 , is given: 
9

[ 1 0 ] [ ] 1 ( ,0) exp( ) exp ( ) KZ I I K ZLT ZLT S u =− −  − − d Z (5C) 

and the hydraulic conductivity solution for the unsaturated zone is given: 

1 0 1

1

1 1 2

1 1 2

1

( , ) [ exp( )]exp[ ( )]

4( )exp exp 2 4

sin[ ( )]sin( ) exp

1 2

2

nZ nZ S u

nZ ZLT

m u m u

m

m u

m u

KZt I I K d Z

Z t I I D

dZ d D t d d





 



 





∞

=

=− − − −

    − −   −      

Λ− Λ ⋅   −Λ  

+ +Λ

∑

(5D) 

The values of Λmare the positive roots of the pseudoperiodic characteristic equation, 

1 tan( ) 2 0 u Λ + Λ=  d . (6) 

In TRIGRS, we solve for values of Λm by a combination of bracketing, bisection, and NewtonRaphson iteration (Press and others, 1986). 

Early-Time Solutions for Constant Surface Flux

The F

D t 0.05  ≤

ourier series in equations 5A and 5D converge slowly for very early times, and for 

values of , we substitute the following expression for equation 5A: 

1

1 1

1 1

( )exp 2

exp erfc ( ,) 2 4 2 4

(2 ) (2 ) exp erfc

2 4 2 4

nZ ZLT

u

u u

Z I I

Z Z t t D D qd t D t

d Z t t d Z D D

D t

 



 





 

 

      −    

        + +   =            

    − −   −+ +              

(7A) 
10

( )

1

1

1

1 1

1

1 1

( )exp[ ( )] 2( )exp 2

erfc

2

exp erfc , 2 4 2 4

(2 ) erfc

2

(2 ) (2 ) exp erfc

2 4 2 4

ZLT ZLT s u nZ ZLT

u

u u

Z I I K dZ I I

Z

D t

Z Z t t D D K Zt D t

d Z

D t

d Z t t d Z D D

D t



 





 









 



 

  − − − −+ −    

       

    −+ +   =      

  − −      

    − − ++ +         

 

 

 

 

 

 

    

(7B) 

Later-Time Solutions for Variable Surface Flux

To accommodate the time-varying input at the upper boundary, such as would occur during 

a rainstorm with varying intensity, we have put the solutions for unsaturated zone infiltration into 

Heaviside series expansions. 

We use the following formula for the basal flux arriving at the base of the unsaturated zone; 

it is a generalized form of Srivastava and Yeh’s (1991) solution given above as equation 5A. 

1

1 2

1

1 1 2

1

( ) 4( )exp exp 2 4

( ,) ( ) sin( ) exp ( )

1 2

2

u n

nZ nZ ZLT

N

u n m mu n m n

m u

m u

d tt I II D

qd t Ht t d Dtt d d













∞

=

=

     − − − −         

= −   Λ Λ   ⋅   −Λ −  

+ +Λ  

∑ ∑

1 1

1 1 2

1 1

1 1 2

1

( ) 4( )exp exp 2 4

( ) sin( ) exp ( )

1 2

2

u n

nZ nZ ZLT

N

n m mu n m n

m u

m u

d tt I II D

Ht t d Dtt d d













+

∞ +

= +

=

     − − − −         

− −   Λ Λ   ⋅   −Λ −  

+ +Λ  

∑ ∑ (8A) 

Computation of K Z( ) ,t is based on a generalized form of the solution of Srivastava and Yeh 

(1991): 
11

1 0 1

1

1

1 1 2

1 1 2

1

[ exp( )]exp[ ( )]

( ) ( , ) ( ) 4( )exp exp 2 4

sin[ ( )]sin( ) exp ( )

1 2

2

nZ nZ s u

N

n

n nZ ZLT

n

m u m u

m n

m u

m u

I I K d Z

Z t t KZt Ht t I I D

dZ d Dtt d d





 



 





=

∞

=

 

 

− − − −

    − = −− −     −       

Λ− Λ ⋅   −Λ −  

+ +Λ  

∑

∑

1 0 1

1 1

1

1

1 1 2

1

1 1 2

1

[ exp( )]exp[ ( )]

( ) ( ) 4( )exp exp 2 4

sin[ ( )]sin( ) exp ( )

1 2

2

nZ nZ s u

N

n

n nZ ZLT

n

m u m u

m n

m u

m u

I I K d Z

Z t t Ht t I I D

dZ d Dtt d d





 



 





+

+

=

∞

+

=

 

 

− − − −

    − − − −−     −       

Λ− Λ ⋅   −Λ −  

+ +Λ  

∑

∑

(8B) 

Pressure head in the unsaturated zone is obtained from equations 4B and 8B. 

0

1

cos ( , ) ( , ) ln

s

KZt Z t

K



 



 

= +    

(9) 

The user may specify one of two values for 0 , either 0 = 0, for computation at the water table in 

agreement with the solutions of Srivastava and Yeh (1991), or 0 −= 1  , for computation at the 

top of the capillary fringe as suggested by Savage and others (2004). When using 0 = 0, the depth 

of the water table, d, is automatically substituted in place of du in all formulas presented here. As 

a result of the greater depth involved, using 0 = 0 results in a slower rise time for the water table. 

The volumetric water content is computed by substituting the results of equation 9 into equation 

4C. 

Water Table Rise

Water-table rise occurs in response to vertical infiltration when the amount of infiltrating 

water reaching the water table exceeds the maximum amount that can be drained by gravity. For 

the one-dimensional geometry assumed by our infiltration models, the long-term vertical hydraulic 

gradient is  −1, the maximum rate of drainage, qZ max , is 
12

Z max Kq S ( )  −= 1 (10A) 

and excess flux, qZex , is determined by the following formulas:

qZex = 0, if ( ) ( ) max , u d Z ZLT qd t c q I ≤ − (10B) 

or 

( ) ( ) max , Zex u d Z ZLT q qd t c q I =− − , if ( ) ( ) max , u d Z ZLT qd t c q I > − (10C) 

where the flux exiting from the base of the unsaturated zone, ( )tdq u , , is given by equation 8A. The 

constant, d c , is 1 for the case where the basal no-flow boundary is at infinite depth, and is 

arbitrarily set to 0.1 when the basal no-flow boundary is at finite depth. If the user selects a slopeparallel final pressure distribution, steady background flux is ignored, I ZLT = 0 ; in addition, 

selecting hydrostatic-pressure distribution prevents all basal drainage, so that qZ max = 0 in 

equations 10A–10C. 

The equations governing transient water-table rise are nonlinear and we are unaware of any 

techniques for solving the equations analytically to produce a continuous solution for pore pressure 

above and below the water table. Such problems are usually solved by numerical methods (see 

Freeze and Cherry, 1979; Huyakorn and Pinder, 1983). Consequently, we use approximate 

formulas to estimate the rate of water-table rise and pressure head in the vicinity of the rising water 

table.

To compute the water-table rise, we compare the excess flux accumulating at the water 

table (for 0  = 0 ) or top of the capillary fringe (for 0   = −1/ ) to available pore space directly 

above the water table or capillary fringe (fig. 7A–C). Although some air remains trapped during 

infiltration, we assume complete saturation below the rising water table in our computations that 

follow. We note here that added uncertainty results for the case 0   = −1/ , because the capillary 

fringe is imperfectly saturated and its height is dynamic during infiltration. We integrate the excess 

flux arriving at the water table (equation 10B or 10C) over time to compute the volume of water 

(per plan-view unit area) accumulating at the base of the unsaturated zone, ( ) V t A . The limits of 

integration are from t = 0 to any time, t , that is of interest during the simulated infiltration event. 

Numerical integration of qZex by repeating Simpson’s three-point rule over successive time steps 

yields VA ( )t ; numerical integration is necessary because the analytical solution for the time integral 

of equation 5A converges poorly. We compute the fillable volume (available pore space) per unit 

(plan-view) area, Vf( ) Z,t , above the water table by integrating the fillable porosity  f

(  f s = − ) from the initial top of the saturated zone up to various equally spaced heights above it 

in the unsaturated zone. (All pore space is considered to be filled in the capillary fringe directly 

above the water table.) The result is shown: 
13

1 ( ,) ( ) ( ,)

u

Z

f sr u

s d

V Zt d Z K ZtdZ

K

 

 

= − −+  

  ∫ (11A) 

For convenience, we define, 

1 ( ,) ( ,)

u

Z

s d

T Zt K ZtdZ

K = ∫ . (11B) 

In equation 11B, TZt ( ,) has units of length [ L ] and ( ,) KT Zt S is conceptually similar to the 

transmissivity of a layer between the top of the saturated zone and some depth Z in the overlying 

unsaturated zone. In figure 7C, the roughly triangular areas between the water-content curve 

labeled t

1, the dashed vertical line labeled s , and the dashed horizontal lines labeled Zk or Zk+1

represent the volume Vf( ) Z,t in equation 11A. We obtained formulas for T Z( ) ,0 and T Z( ) ,t by 

integrating ( ) ,0 / KZ KS and ( ) , / K Zt KS with respect to Z, from the top of the saturated zone, 

u d , to Z: 

1

1 0

1

1 1 exp[ ( )] ( ,0) ( ) [ exp( )] u

ZLT u ZLT s

s

d Z TZ I d Z I K

K



 



  −− − = −− −    

 (12A) 
14

( )

1

1 0

1

1

2

1

1 2 2 1

1

1 1

1 exp[ ( )] ( ) ( ) [ exp( )]

1 2( ) ( ) ( , ) exp 4

sin( )exp ( )

1 1 2

4 2

2 exp exp 2 2

u

n nZ u nZ s

nZ ZLT n

s

mu m n

m u

m m u

u

m

d Z Ht t I d Z I K

I I t t TZt D

K

d Dtt

d d

d Z







 











 

∞

=

  −− − − −− −    

  − −   =− −        

Λ −Λ −

     Λ+ + +Λ      ⋅    Λ −   ⋅  

∑

1

1 1 [sin( ( )) 2 cos( ( ))]

N

n

mu m mu   d Z d Z

=

 

 

 

 

 

 

        

  ⋅ Λ − +Λ Λ −    

∑

( )

1

1 1 0

1

1

1

2

1 1

1 2 2 1

1

1 1

1 exp[ ( )] ( ) ( ) [ exp( )]

1 2( ) ( ) exp 4

sin( )exp ( )

1 1 2

4 2

2 exp exp 2 2

u

n nZ u nZ s

nZ ZLT n

s

mu m n

m u

m m u

u

m

d Z Ht t I d Z I K

I I t t D

K

d Dtt

d d

d Z







 











 

+

+

∞ +

=

  −− − − −− −    

  − −   −− −        

Λ −Λ −

     Λ+ + +Λ      ⋅    Λ −   ⋅  

∑

1

1 1 [sin( ( )) 2 cos( ( ))]

N

n

mu m mu   d Z d Z

=

 

 

 

 

 

 

        

  ⋅ Λ − +Λ Λ −    

∑

 (12B) 

No early-time equivalent of equation 12B is available. We use numerical integration of the results 

of equation 8B for rare cases in which early-time results are needed for computation of the watertable rise. 
15

The water-table rise is first bracketed between the two nearest depths, Zk and Zk+1, at which 

the available volume has been computed, such that Vf Zk ( ) ,t ≥ VA ( )t ≥ Vf Zk+1 ( ) ,t , and then it is 

estimated by linear interpolation between Zk and Zk+1(fig. 7C). Due to the slope of the water 

content curve directly above the capillary fringe (fig. 7B), the method outlined in figure 7C may 

slightly overestimate the actual rise of the water table, ∆h. The error can be reduced by reducing 

the depth increment, Zk +1 − Zk . 

Pressure Diffusion Below the Water Table

Once the water-table rise has been computed at a grid cell for time increment n, we apply 

the load (as measured by the computed water-table rise, hn ) at the initial top of the saturated zone 

and compute the pressure head using formulas for a constant head boundary condition that have 

been adapted from analogous heat-flow problems (Carlslaw and Jaeger, 1959) and generalized for 

time-varying water-table rise, hn . For consistency with equation 3A, pressure head, hn , applied 

by the newly accumulated water above the initial water table is hn n =  h . 

Infinite-Depth Basal Boundary, Variable Applied Head

The computer program TRIGRS uses pressure-diffusion formulas for computing head 

changes below a rising water table. For relatively homogeneous hillsides, a model having an 

impermeable basal boundary at infinite depth can approximate the pressure rise.

1

1 1 1 1 1

( , ) ( )erfc ( )erfc 2( ) 2( )

N N

w w

w hn n hn n

n n n n

Z Z Z t Ht t Ht t

Dt t Dt t    +

= = +

   =− −−    − −    ∑ ∑ (13) 

In equation 13, Zw is the vertical distance below the initial water table, Z Zd w = − . 
16

Finite-Depth Basal Boundary, Variable Applied Head

For hillsides that have a strong permeability contrast at shallow depth, such as many 

hillsides that have a colluvial mantle only a few meters thick, a model that has an impermeable 

boundary at finite depth can estimate pressure rise. 

Early-Time Formula

As with the infinite depth model, a solution that provides accurate results for early and 

intermediate times is based on the complementary error function: 

1

1 1

1 1

1 1

(2 1) ( 1) erfc 2( )

( ,) ( )

(2 1) ( 1) erfc 2( )

m LZ w

N m n

w hn n

n m LZ w

m n

md Z

Dt t

Z t Ht t

md Z

Dt t

 

∞ +

=

= ∞ +

=

    − −   −   −     = −  

    − +

+ −   −    

∑

∑

∑

1

1 1 1

1

1 1

1 1 1

(2 1) ( 1) erfc 2( )

( )

(2 1) ( 1) erfc 2( )

m LZ w

N m n

hn n

n m LZ w

m n

md Z

Dt t

Ht t

md Z

Dt t



∞ +

= +

+

= ∞ +

= +

    − −   −   −     − −  

    − +

+ −   −    

∑

∑

∑

(14) 

In equation 14, dLZ denotes the vertical height or thickness of the saturated layer, from the initial 

water table down to the impermeable boundary. 

Later-Time Formula

At later times, we use a formula that is based on a Fourier series solution.

2 2

1 1

2

1

1

4 1 (2 1) ( ) 1 ( 1) exp 21 4

( ,) ( )

cos (2 1) 1 2

m n

N m LZ

w hn n

n w

LZ

m Dt t

m d

Z t Ht t

Z

m

d





 



∞ −

=

=

    − −   −− −  −   = −  

      ⋅ −−          

∑

∑

2 2

1 1 1

2

1

1

1

4 1 (2 1) ( ) 1 ( 1) exp 21 4

( )

cos (2 1) 1 2

m n

N m LZ

hn n

n w

LZ

m Dt t

m d

Ht t

Z

m

d









∞ − +

=

+

=

    − −   −− −  −   − −  

      ⋅ −−          

∑

∑ (15) 

The step-function series in equations 13–15 contain N successive time steps to represent 

water-table rise from the initial base of the unsaturated zone. The time steps correspond to those in 
17

the unsaturated infiltration series. All time steps are completed at one grid cell before proceeding 

to the next cell. 

Infiltration, Runoff, and Flow Routing

The program, TRIGRS, uses a simple method for routing of surface runoff from cells that 

have excess surface water to adjacent downslope cells where it can either infiltrate or flow farther 

down slope.3

Ks

 We assume that runoff occurs when the precipitation and runoff supplied to a cell 

exceed its infiltrability. The saturated hydraulic conductivity, , generally equals the infiltrability, 

i, for saturated and tension-saturated soils (Hillel, 1982; Iverson, 2000). The purpose of routing 

the surface runoff is to prevent the loss of excess precipitation that cannot infiltrate at the cell of 

origin and to improve the performance of the model in urbanized or other areas where pavement or 

other impervious surfaces exist. However, the user may bypass the runoff routing, if desired. 

Running TRIGRS without routing would be appropriate where a complex system of storm drains 

complicates runoff routing. 

We compute the infiltration, I, at each cell as the sum of the precipitation, P, plus any 

runoff from upslope cells, Ru , with the limitation that infiltration cannot exceed the saturated 

hydraulic conductivity, Ks: 

I = P+ Ru , if ≤+ KRP su (16A) 

or 

= KI s , if >+ KRP su (16B) 

At each cell where P + Ru exceeds Ks the excess is considered runoff, Rd , and is diverted to 

adjacent downslope cells.

Rd = P + Ru − Ks, if KRP su ≥−+ 0 (17A) 

or 

= 0 Rd , if <−+ 0 KRP su (17B) 

Overland flow between adjacent cells is assumed to occur instantaneously; we have made 

no attempts to model the rate of overland flow. Consequently, individual storm periods should be 

long enough to allow surface water to flow to adjacent cells. The routing method enforces mass 

balance for each time step through the storm; the total precipitation at all the cells must equal the 

water that infiltrates at all the cells and the water that flows to edges of the problem domain (or any 

closed depressions) without infiltrating. Thus, excess precipitation that cannot infiltrate at its cell 

of origin has the opportunity to infiltrate at some point down slope. However, TRIGRS does not 

carry runoff over from one time step to the next, nor does it track water that enters storm drains. In 

other words, water that runs off a cell during any given time step will either infiltrate at another cell 

or reach the edge of the model (or some other sink) within that time step. 

 

3

Runoff routing is applied only to the transient part of the solution. If TRIGRS detects that the steady infiltration rate 

at a cell exceeds the saturated hydraulic conductivity, it reduces the value of Is Ks to a value that makes  = 0 at the 

cell, with the result that pressure head above the user-specified steady water table is zero. Cells where this occurs are 

identified in the log file, so that the user can go back and correct the input data.
18

In addition to runoff from cells where precipitation exceeds the infiltrability, water is 

assumed to run off from any cells where the water table is initially at the ground surface and the 

initial (steady) rate of infiltration is negative. Water is exfiltrating from such cells and runs to the 

adjacent downslope cell(s). The TRIGRS program tracks exfiltrating water in the mass balance 

calculations.

By computing the infiltration and runoff at cells in order from the topographically highest 

cell(s) to the lowest, the program TRIGRS avoids the necessity of iterating to reach a solution that 

satisfies mass balance. Mass-balance calculations for each storm time step verify that all water 

applied as input is accounted for as either infiltration or runoff. Mass balance will be achieved only 

if the topographic data have been correctly indexed. A digital elevation model (DEM) can be 

indexed correctly only if it has been adjusted in GIS software to be hydrologically consistent. This 

adjustment is a selective smoothing process accomplished by raising the elevation of single-cell 

closed depressions to match surrounding cells, and slightly raising or lowering the elevation of cells 

in flat areas to produce flow directions that are consistent with surrounding topography. The 

TRIGRS program writes the results of the mass-balance calculations to a log file, named 

"TrigrsLog.txt", in the same folder as the application. 

Topographic Indexing

We include a companion utility program, TopoIndex (Topographic Index). TopoIndex 

prepares a group of data files for the runoff-routing calculations. These output files contain a list of 

cell numbers ranked in order from highest to lowest elevation and a grid of cell numbers indicating 

the neighboring cell on the steepest downslope path. Figure 10 illustrates the grid cell numbering 

scheme. The TRIGRS program uses these data to control the order of runoff routing calculations as 

described above. If any of these files are missing, then TRIGRS will compute infiltration based on 

the rainfall input without routing any excess surface water. The output files of TopoIndex also 

contain a list of grid cells and their downslope neighbors, along with a list of weighting factors that 

determine what proportion of the runoff is transferred to each neighboring downslope cell.

The program TopoIndex uses a two-step process to index the cells from highest to lowest 

elevation. The first step, sorting the cell elevations, consists of defining an index using the 

“Heapsort” algorithm (Press and others, 1986). The second step corrects the index by comparing 

the indices of neighboring cells to their relative positions along downslope paths. If the indices of 

two neighboring cells that have the same apparent elevation (due to numerical rounding or 

truncation) are out of order, they are swapped. This iterative step continues until the number of 

corrections on the last pass equals zero. Figure 11 compares grid cell elevation (fig. 11A) and cell 

index (fig. 11B) for a sample grid. We have successfully tested TopoIndex on a grid of a 

approximately 2 million cells. The TRIGRS program uses output files from TopoIndex to 

implement the runoff-routing scheme described in the previous paragraphs. Our tutorial (appendix 

2 of this report) describes how to operate the utility program and prepare data for the analysis. 

Figure 12 shows the designation of D8 flow directions used by TopoIndex in contrast to the ESRI 

direction codes. Successful operation of TopoIndex requires the use of a hydrologically consistent 

DEM as input. TopoIndex does not account for storm drains and similar artificial drainage 

diversions that are not expressed in the topography. 

TopoIndex allows the user considerable latitude in determining how excess water, Rd , will 

be distributed among the downslope cells. Table 1 outlines the available choices of distribution 

patterns, which (in this implementation) apply only to square cells. Many of the patterns depend on 

an exponent, , and the slopes between the current grid cell and its downslope neighbor cells 

similar to the method of Freeman (1991). The same exponent is used in computing weighting 
19

factors throughout the grid. Figure 13 illustrates how the choice of distribution patterns affects the 

output. Figure 13A and 13B show cell numbers and elevations from part of the sample grid. The 

default case ( > 20 , fig. 13C) is to divert excess water only to the adjacent cell on the steepest 

downslope path (the D8 method of O'Callaghan and Mark, 1984). Letting  = 0 distributes the 

excess water evenly among all adjacent downslope cells; this pattern has the greatest dispersion of 

all the choices (fig. 13E). If  =1 (fig. 13F), then flow is distributed among all downslope cells in 

proportion to the slope, which is equivalent to the method of Quinn and others (1991). Tarboton 

(1997) showed that this method has significant dispersion. Choosing higher exponents decreases 

the dispersion and diverts increasing proportions of the flow to the cell on the steepest downslope 

(D8) path (figs. 13G, H). Negative exponents are not implemented because they would concentrate 

flow on flatter slopes. 

If the user enters a negative exponent in the TopoIndex initialization file, it acts as a flag to 

select a method that is similar to Tarboton’s (1997) D-infinity method (fig. 13D). The D-infinity 

method assumes that water flows down the steepest slope, computes the direction of steepest slope 

and attempts to direct flow in that direction by partitioning the flow between the two cells nearest 

to the steepest slope direction. In this case, TopoIndex uses the D8 direction as a starting point and 

compares the slopes of the two neighboring directions to identify the triangular facet that contains 

the steepest slope. If the steepest slope coincides with the D8 flow direction, then all flow is 

diverted to the cell in that direction. Otherwise flow is distributed between the two cells that 

straddle the direction of the steepest slope; the weighting factors are proportional to the angles 

between the grid directions and the dip direction of the facet as described by Tarboton (1997). 

Each of the runoff distribution patterns has its peculiar characteristics. The D8 method is 

the simplest and has the lowest dispersion; however, the resulting runoff distribution is irregular 

and somewhat unrealistic. We prefer the D-infinity method because it is physically more realistic 

than any of the other methods implemented in TopoIndex. The D-infinity method is used in the 

steady-state regional slope stability program SINMAP (Pack and others, 1998). The remaining 

methods, in which the weighting factors are proportional to the slope raised to an exponent, are 

provided to allow the user to achieve greater dispersion or smoothing of the runoff distribution if 

desired. 

Slope Stability 

Following Iverson (2000), we model slope stability using an infinite-slope stability analysis

(Taylor, 1948). In his analysis, failure of an infinite slope is characterized by the ratio of resisting 

basal Coulomb friction to gravitationally induced downslope basal driving stress. This ratio, called 

the factor of safety, FS, is calculated at a depth Z by

tan ( , ) tan ( ,) tan sin cos

W

S

S

c Zt F Zt

Z

 

  

′ ′ ′ − = + (18) 

where c ′ is soil cohesion for effective stress,  ′ is the soil friction angle for effective stress, W is 

unit weight of groundwater, and S  is soil unit weight. The pressure head,  ( ) ,tZ , in equation 18 is 

subject to the limitation imposed by equation 3 and is obtained from either equations 1, 2, 9, or 13–

15 depending on the particular conditions modeled. 
20

Factors of safety, FS , are calculated by equation 18 for transient pressure heads at multiple 

depths Z (fig. 14). Failure is predicted when FS <1 and stability holds where FS ≥1. The state of 

limiting equilibrium exists where = 1 FS . Thus, the depth Z where FS first drops below 1 will be the 

depth of landslide initiation. This initiation depth depends on soil properties and the time and depth 

variation of pressure head, which, in turn, depends on rainfall history. 

In the unsaturated zone we use a simple approximation for Bishop’s (1959) effective stress 

parameter,  , in computing the factor of safety. The approximation ( ) ( ) −−=  rsr , which 

was suggested by Vanapalli and Fredlund (2000), has application to a generalized effective stress 

law and represents a simplified form of the suction-stress characteristic curve (see Lu and Likos, 

2004). To compute the factor of safety above the water table, the matric suction, W tZ ),(  , is 

multiplied by  . In TRIGRS, we compute a depth-averaged soil unit weight for use in equation 18 

to account for the lower unit weight of partially saturated soil above the water table.

Limitations

All numerical and analytical models of ground water flow and slope stability are subject to 

limitations imposed by simplifying assumptions, approximations, and other shortcomings in the 

underlying theories. The TRIGRS program is based on transient vertical infiltration models 

combined with a simple slope stability model and is subject to the same limitations that are inherent 

in the original models (Srivastava and Yeh, 1991; Iverson, 2000; Savage and others 2003): 

• The model assumes flow in homogeneous, isotropic soil; thus, application to areas with marked 

soil anisotropy or heterogeneity in hydrologic properties may cause errors in the solution. 

• The TRIGRS program models only one-dimensional vertical infiltration, which is dominant in

determining pore pressure change early in a storm. During longer storms and periods between 

storms, lateral flow also has a significant influence on the magnitude and distribution of pore 

pressure. Estimates of pore pressure using the TRIGRS model may be in error once lateral flow 

begins to respond to the storm. Consequently, TRIGRS is not suitable for modeling long-term 

effects of alternating periods of rainfall and evapotranspiration, and choosing the correct initial 

conditions for a given storm is critical to obtaining accurate results. The most suitable method 

for handling lateral flow within the TRIGRS framework is currently under investigation. Users 

may estimate “worst-case” pressure-head conditions resulting from late-storm slope-parallel 

flow by setting the initial water table at the ground surface and setting the steady background 

flux to zero or by using the slope-parallel option (tables 2 and 3) to obtain a pore pressure 

distribution equivalent to slope-parallel flow. 

• The results are very sensitive to the initial conditions. Ideally, the user should know the initial 

water-table depth and something about the initial flow field throughout the study area. The user 

must specify the initial depth of the water table and steady infiltration or exfiltration rate at each 

cell or estimated values for the entire area. The initial water-table depth and rate allow the user 

to represent fairly general initial flow conditions. In urban and suburban areas, data from 

geotechnical borings and monitoring wells for environmental remediation projects may provide 

estimates of the initial conditions. In the absence of detailed information, simple ground-water 

flow models that can be adapted to hill-slope conditions might provide reasonable estimates 

(Beven and others, 1984; Beven, 2004). 

• The program uses a simple method for routing surface runoff from cells where precipitation 

exceeds the infiltrability (Ks

 for saturated and tension-saturated soils) to neighboring downslope 
21

cells. This method does not account for evapotranspiration, which is expected to be low during 

most storms. Furthermore, the runoff routing model cannot be used to estimate stream flow 

because runoff from cell to cell is assumed to happen instantaneously and the runoff receives no 

contributions from subsurface flow, except for cells where steady exfiltration is specified in the 

initial conditions. 

• Procedures for correcting and preprocessing digital elevation models (filling pit artifacts or onecell depressions, and computing flow directions) are not included in TRIGRS or TopoIndex. 

However, such procedures are available in modern GIS software. TopoIndex and TRIGRS 

require a DEM that has been corrected to remove noise, such as one-cell depressions. 

• Model results are sensitive to the engineering properties of geologic units in the model area, 

some of which may vary over orders of magnitude. For any geologic material the hydraulic 

conductivity and diffusivity may vary over several orders of magnitude. 

• The slope-stability model is based on an infinite slope analysis, which assumes uniform slope, 

thickness, physical properties, and that pore-water pressure is a function of only depth and time. 

Applying this method over a grid that represents three-dimensional topography may result in 

errors where slope, thickness, physical properties or pore pressure changes abruptly. Nonuniform lateral stresses (or inter-cell forces) and non-planar failures may occur in areas where 

the infinite slope assumptions are violated. Thus, the model is more accurate where conditions 

favor planar failure, which typically occur at depths that are shallow compared to the width of 

the grid cells and in areas where slopes are smooth and have a relatively uniform mantle of 

weathered material. Based on published comparisons of various methods of slope stability 

analysis, we expect the model to compute lower factors of safety than more rigorous methods of 

stability analysis would in areas where failures tend to be non-planar (Duncan, 1996). 

• The TRIGRS model is not recommended for use on slopes steeper than 60°. For slopes steeper 

than about 50°–60°, depending on the depth, cohesion, and other parameters, equation 18 

predicts increasing, rather than decreasing, values of FS for constant depth. The safety

Z cos

 factor 

increases for slopes steeper than 50°–60° because the thickness of the sliding mass, , 

decreases as the slope becomes steeper, with the result that the proportion of shearing resistance 

contributed by cohesion increases significantly over that contributed by friction. Although this 

difficulty can be overcome by dividing a depth grid by cos to effectively transform it to a 

thickness grid and using the transformed grid with TRIGRS, this transformation will not 

overcome the uncertainty resulting from use of a one-dimensional model to represent the 

infiltration process on slopes steeper than 60°–70°. 

• The program does not check for consistent units of measurement between the different input 

parameters. However, we have provided a simple utility program called UnitConvert that can 

be used to convert data contained in grid files between different systems of units. 

• A lack of information about how conditions change once a slope has failed limits the ability of 

TRIGRS to forecast what happens after the factor of safety has fallen below unity (failure) at 

any particular grid cell. Although it is feasible to modify TRIGRS to update the properties and 

hill-slope geometry after failure occurs, the main difficulty lies in predicting what the updated 

properties and geometry should be. In many instances, the scars of shallow failures expose 

more competent material, and considerable time (years to decades) must pass before sufficient 

weathered material accumulates to allow failure at the same location. However, depending on 

cell size, successive failures might occur at neighboring locations within the same cell during a 
22

storm or later storms. Consequently, for computational purposes, we have treated cells as if 

they are unchanged by failures during a given simulation—in other words, TRIGRS does not 

modify the soil thickness or properties when FS < 1. 

Program Operation 

The Windows versions of the programs TRIGRS, TopoIndex, Gridmatch, and UnitConvert 

are launched by double clicking their icons. These programs have a command line interface, and 

run in their own console windows. Mac OS versions are started from the command line in the 

Terminal application.

4

TRIGRS Initialization File

 Upon starting, TRIGRS opens a user-supplied initialization file found in the 

present working directory, and reads a list of file names and parameters needed to compute the pore 

pressures and factors of safety in each cell. The program then opens and reads files of input data, 

computes the results, and saves the output to a series of files. The program exits automatically 

when it is finished. The TRIGRS program writes a brief log file that may be useful in 

troubleshooting or in keeping a record of successive analyses. The log file includes a record of the 

starting and ending time, any errors that occurred, a copy of data from the initialization file, 

statistics on convergence of series solutions, and other information useful in interpreting program 

output. The utility programs, TopoIndex, GridMatch, and UnitConvert, operate in the same 

manner. Successful operation of the programs requires that the initialization file and other input 

files be organized and formatted correctly. Therefore, we describe the format and contents of these 

files before presenting a brief tutorial (appendix 2). 

The initialization file (table 2) controls input and output. Use of an initialization file 

accelerates analysis when compared to the alternative of repeatedly choosing input files from 

interactive dialog boxes. The initialization file is in ASCII text format and contains a list of header 

lines, input file names, and constants needed to complete the analysis. To change problem 

parameters, the user edits the initialization file using a text editor. Consequently, most runtime 

errors result from mistakes in editing the initialization file or errors in the input grids. To reduce 

errors in editing the file, a descriptive or explanatory heading precedes each line (or in one case a 

group of related lines) of data in the initialization file. The first line of data contains a project title. 

It is followed by several lines of numeric data including program parameters, properties for each 

zone, rainfall infiltration rates and times; the second group of lines contains names of input files; 

and the third group of lines specifies which output files will be saved, as well as their file names 

and paths. The fourth and final group contains user options that affect how the analysis is 

performed. Table 2 shows the organization of the initialization file and defines the input file names 

using partial path names (starting at the current folder). The next few paragraphs describe the 

structure of the initialization file and general characteristics of the data. We have provided a more 

detailed explanation of the initialization data needed to perform a simulation in the tutorial 

(appendix 2). 

Numeric data.—The file contains several lines of numeric data needed to determine array 

sizes, model parameters, precipitation history, initial conditions, and physical properties. An 

 

4 To run TRIGRS from the command line, start Terminal, type cd and drag the (main) TRIGRS folder icon to the 

Terminal window. After the path name of the folder appears on the command line, press return. Drag the TRIGRS 

program icon (in the folder TRIGRS/bin) to the terminal window and press return. TRIGRS will launch immediately 

and look for the initialization file in the present working directory, which in this case is the folder selected using the cd

command. 
23

explanatory heading precedes each line of data. The file is organized with a project title and 

heading on the first two lines, numerical input data and headings on the next six lines, followed by 

two lines of headings, one line of data for each property zone, and two more lines each for rainfall 

rate, corresponding times, and their respective headings (table 2). The user may set physical 

parameters, such as soil and water table depths, background flux, and rainfall (infiltration) rates to a 

constant value over the entire grid by specifying a positive value. A negative value for any of these 

parameters prompts the program to read an input grid for spatially variable data. The user 

determines the number of property zones and assigns physical properties (soil strength parameters 

and hydraulic properties) to each zone (one line per zone). A single property-zone grid maps all 

these properties to the model. Table 3 defines the numeric variables listed in table 2 and indicates 

appropriate ranges of values. 

Input files.—The user must specify a file name for each input file. The TRIGRS program 

automatically reads a name for each potential input file; therefore the user must provide a 

placeholder for any input file that either does not exist or will not be used for a given simulation. 

The TRIGRS program can accept either full or partial path names5

In table 2, path names of input files are on the next group of lines after the numeric data 

described previously. Each file name occupies a single line in the file. The total number of lines in 

the file depends on the number of property zones and number of periods of different rainfall 

intensity. Be sure to include a file path name, rifil(), for each period, but only one heading for this 

group of file names. In this example, the rainfall rates during the storm are set constant over the 

grid. Two storm periods are specified, so two lines are needed for the rainfall rate input file names

(placeholders in this case). 

as long as 255 characters; 

however, the rules for specifying path names vary between operating systems. If the paths are 

typed incorrectly, the program will report an error. To minimize errors in specifying file path 

names, the user may place all the input data files needed for a given simulation into a specific 

folder. On Windows systems, placing a copy of the executable program (TRIGRS.exe) and 

initialization file into the folder where the data reside will further simplify operation. By doing so, 

the user needs only to specify the file names and does not need to worry about the paths. 

Output files.—Table 4 lists the names and describes the contents of the output files that 

TRIGRS can generate. A series of eight logical values (true or false) determines which output grid 

files are saved to disk. These are listed directly after the input file names (table 2). A numeric flag 

at the end of the list of output file choices determines whether the cell-by-cell list file is saved and 

whether it contains optional data. Note that this file can become very large, particularly if the 

analysis includes multiple time steps. This option should be reserved for use with small test grids 

to prevent filling the disk drive. The user may specify an identification code to group files from a 

particular simulation. Users may also specify a separate folder where the output files can be saved. 

Certain output grids can be saved at selected times during the simulated storm. 

Model options.—The last few lines of the initialization file (table 2) allow the user to select 

options that control how the model performs certain functions in the unsaturated zone (appendix 2). 

These options allow the user to compute the fillable porosity either analytically or numerically, 

compute positive pressures above the original water table as water accumulates at the base of the 

unsaturated zone, and determine how the capillary fringe is treated. The final option allows the 

user to select the slope of the beta line, which determines maximum pore pressure at each depth. 

On start-up, TRIGRS will look for the default initialization file, "tr_in.txt", in the same 

folder where the program file resides. If it is not present, then TRIGRS prompts the user to enter 

 

5 The full path is the concatenation of the drive name, the folder name(s), and the file name. For example, C:\\Program 

Files\TRIGRS\TRIGRS.EXE is a full path name and TRIGRS\TRIGRS.EXE is a partial path name. Parent directories 

can be specified by typing "..\" in Windows or ".." in Mac OS and UNIX.
24

the name of the initialization file. TopoIndex, GridMatch, and UnitConvert follow a similar 

convention. 

Input File Format

Most input files contain gridded data and all grid files use the same basic format, consisting 

of a short list of file parameters followed by the grid data. File parameters can be entered in 

formats that support either ESRI ArcGIS or the open-source GRASS GIS. The files are in spacedelimited ASCII text format and can be viewed or edited with a text editor. 

The first six lines of the ArcGIS ASCII grid file contain parameters defining the grid. Table 

5A shows an elevation grid from the tutorial data folder. The first two lines define the number of 

columns and rows in the grid. The third and fourth lines indicate the coordinates of the lower-left

(southwest) corner of the grid, the fifth line indicates the grid spacing, and the sixth line defines the 

value used to identify cells where no data exists. Each parameter is preceded by a text label up to 

14 characters long (table 5A). Text labels for the ArcGIS ASCII grid format should appear exactly 

as shown in table 5A. The remaining lines (line 7 and beyond) contain the value (topographic 

elevation in this example) of each grid cell. The grid data begins at the northwest corner cell in the

upper left-hand corner of grid and the cells are listed from left to right (west to east), and row-byrow from top to bottom (north to south). The grid data may be written in either decimal or 

scientific floating-point format. Single-precision arrays hold the grid data. Spaces, commas, or 

tabs may be used as delimiters. The default "no data" (null) value for cells to be excluded from the 

analysis is -9999. 

The header for a GRASS GIS two -dimensional ASCII grid format also lists six lines (table 

5B). The first four lines list the north, south, east, and west boundaries of the grid in UTM or 

similar coordinates, and the last two lines contain the number of rows and columns in the grid. 

Grid values are arranged and formatted the same way as for the ArcGIS files. The GRASS GIS 

software allows the user to specify the value for no-data cells, but the header does not specify this 

value, so we suggest setting the no-data (null) value to -9999 to agree with the default value in 

TRIGRS.

The only exceptions to the grid-file format are files used to control the routing of excess 

water. These text files include three lists as shown in tables 6, 7, and 8. Table 6 illustrates the 

format of a file containing a list of cells downslope from each grid cell. Table 7 follows a similar 

format and contains the weighting factors for the corresponding cells listed in table 6. Table 8 

illustrates the format of the index file that indicates the order in which cells are analyzed to 

correctly partition infiltration and runoff. The input routing method also uses a grid file that 

indicates which neighboring cell is on the steepest downslope path. 

Output File Formats

The TRIGRS program automatically generates a log file and can generate optional output 

files (if specified in the initialization file, table 2). The log file (always named "TrigrsLog.txt") 

contains information about program operation and repeats the data from the initialization file. By 

renaming the file (add the identification code applied to the other output files from the simulation), 

the user can save a record of the simulation to keep with the other output files. Most of the optional 

output files follow the same grid format as the input files. Three of the files contain (1) the 

minimum factor of safety, (2) the pressure head, and (3) the depth at which the minimum occurs. 

In case of runoff, TRIGRS saves a file showing runoff at the end of each storm period. 

Another optional output file contains a cell-by-cell listing of the pore pressure and factor of 

safety at each depth in the profile (fig. 14 and table 9). This file can become quite large because its 
25

size increases with the user-specified number of vertical increments (fig. 14) and with the number 

of cells in the grid. For example, on a grid with 20,522 cells, the profile data using five vertical 

increments (six points per cell) occupied about 5 MB, where the ASCII grid output files occupied 

368 KB each. The suggested use of this file is for tracking time histories of a few representative 

cells organized into a small grid for model calibration.

TopoIndex Initialization File

TopoIndex requires less input than TRIGRS, but its initialization file is similar in design. 

Table 10 shows an example initialization file for TopoIndex. The first two lines contain numeric 

data needed to control program operation and establish array sizes. The next two contain the names 

of input files and the remaining lines control output. Table 11 lists the output files. TRIGRS uses 

four of the seven output files that TopoIndex is capable of generating (indicated by * in table 11). 

The other three are provided as a convenience for the user in visualizing the data.

Troubleshooting

The following suggestions will help to ensure that TRIGRS and associated programs work 

as expected.

• Use GIS software to hydrologically correct (remove pit artifacts from) the digital elevation 

model (DEM).

• Make sure that all grid files (1) are in one of the specified formats, (2) have the same 

numbers of rows and columns, and (3) have the same number of no-data cells in the same 

locations (fig. 10). We have provided a utility program “GridMatch” to help the user verify 

that all input grids have the same structure. It checks the number and locations of no-data 

cells as well as verifying that the number of rows and columns are the same. Table 12 is a 

sample listing of the initialization file for GridMatch. Results of the file comparison are 

saved in the file gridmatchLog.txt. Sample mismatching grid files, used to demonstrate the 

operation of GridMatch, are found in the folder Data/gridmatch/ within the TRIGRS 

distribution. 

• Check the path names of all files listed in the initialization file for correct location and 

spelling. 

• Be sure that the initialization file is in the same folder as the program file on Windows 

systems or in the present working directory on Mac OS X and other UNIX-based systems. 

• Save and close the initialization file and all other input files before attempting to run 

TRIGRS or one of the associated utility programs 

A list of common problems and solutions appears in table 13. 

Comparison of Computed Results with Published Analyses

We tested TRIGRS against Iverson’s (2000) examples from the Minor Creek landslide and 

from an experiment at the USGS debris-flow flume and were able to reproduce his results for pore 

pressure and factor of safety profiles. Thus, we were able to verify that the computer code 

accurately implements Iverson’s (2000) formulas for pore pressure and factor of safety. After 

publishing the original version of TRIGRS, students of Ricardo Rigon (University of Trento, 

Trento, Italy) pointed out an error in the coordinate transformation in Iverson’s (2000) original 

paper. Correction of that simple error causes our numerical results to differ somewhat from 

Iverson’s (2000) published values, though preserving their basic form. Folders named 
26

“MinorCreek” and “flume” contain files with necessary data to reproduce these results. By varying 

elapsed time, t, in the initialization file, the user can compute pressure head and factor of safety 

profiles like those of figures 7, 8, 10, and 11 of Iverson (2000). Note that the computed profiles 

conform to the limit imposed by equation 3. We have also verified that the new version of 

TRIGRS described in this report reproduces the published results of Srivastava and Yeh (1991). 

Data in the folder “sy91” produce profiles that allow the user to reproduce the curves plotted in 

figure 3 of Srivastava and Yeh (1991). 

Acknowledgments 

Mark Reid and Jon McKenna (both USGS) provided critical reviews of the original 

program and documentation. John Michael (USGS) prepared GIS grids used in development and 

testing of the programs. Bill Schulz (USGS) tested the tutorial and software and provided helpful 

suggestions. Diana Salciarini (University of Perugia, Italy) provided suggestions for improving the 

user interface, and Silvia Simoni (University of Trento, Italy) made us aware of minor errors in the 

original code. Carlo Sica (University of Salerno, Italy) and Dianne Brien (USGS) provided critical 

reviews of the present version.

References Cited

Baum, R.L., Savage, W.Z., and Godt, J.W., 2002, TRIGRS—A Fortran program for transient 

rainfall infiltration and grid-based regional slope stability analysis: U.S. Geological Survey 

Open-File Report 02-0424, 35 p., 2 appendices. 

Beven, K.J., 2004, Rainfall-runoff modeling: The Primer: Chichester, United Kingdom, Wiley, 372 

p. 

Beven, K.J., Kirkby, M.J., Schofield, N., and Taqq, A.F., 1984, Testing a physically-based flood 

forecasting model (TOPMODEL) for three UK catchments (Crimple Beck, Hodge Beck, Wye): 

Journal of Hydrology, v. 69, p. 119–143. 

Bishop, A.W., 1959, The principle of effective stress: Teknisk Ukeblad, v. 39, p. 859–863. 

Carslaw, H.S., and Jaeger, J.C., 1959, Conduction of heat in solids: New York, Oxford University 

Press, 510 p. 

Cody, W.J., 1969, Rational Chebyshev approximations for the error function: Mathematics of 

Computation, v. 23, p. 631–638. 

Duncan, J.M., 1996, Soil slope stability analysis, in Turner, A.K., and Schuster, R.L., Landslides—

Investigation and mitigation: Washington, D.C., National Academy Press, Transportation 

Research Board Special Report 247, p. 337–371. 

Freeman, T.G., 1991, Calculating catchment area with divergent flow based on a regular grid: 

Computers and Geosciences, v. 17, no. 3, p. 413–422. 

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: Englewood Cliffs, N.J., Prentice Hall, 604 p. 

Gardner, W.R., 1958, Some steady-state solutions of the unsaturated moisture flow equation with 

application to evaporation from a water table: Soil Science, v. 85, p. 228–232. 

Hillel, Daniel, 1982, Introduction to soil physics: San Diego, Calif., Academic Press, 364 p.

Huyakorn, P.S., and Pinder, G.F., 1983, Computational methods in subsurface flow: New York, 

Academic Press, 473 p. 

Iverson, R.M., 2000, Landslide triggering by rain infiltration: Water Resources Research, v. 36, no. 

7, p. 1,897–1,910. 

Lu, Ning, and Likos, W.J., 2004, Unsaturated soil mechanics: New York, Wiley, 584 p. 
27

O'Callaghan, J.F., and Mark, D.M., 1984, The extraction of drainage networks from digital 

elevation data: Computer Vision, Graphics and Image Processing, v. 28, p. 328–344. 

Pack, R.T., Tarboton, D.G., and Goodwin, C.N., 1998, The SINMAP approach to terrain stability 

mapping, in International Congress of the International Association of Engineering Geology and 

the Environment, 8th, Vancouver, British Columbia, Canada, September 21–25, 1998, 

Proceedings: Rotterdam, Netherlands, A.A. Balkema, v. 2, p. 1,157–1,165. 

Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical recipes, the 

art of scientific computing: New York, Cambridge University Press, 818 p.

Quinn, P., Beven, K., Chevallier P., and Planchon, O., 1991, The prediction of hillslope flow paths 

for distributed hydrological modeling using digital terrain models: Hydrological Processes, v. 5, 

p. 59–80. 

Savage, W.Z., Godt, J.W., and Baum, R.L., 2003, A model for spatially and temporally distributed 

shallow landslide initiation by rainfall infiltration, in Rickenmann, D., and Chen, C., eds., Debrisflow hazards mitigation—mechanics, prediction and assessment: Rotterdam, Millpress, p. 179–

187. 

Savage, W.Z., Godt, J.W., and Baum, R.L., 2004, Modeling time-dependent aerial slope stability, 

in Lacerda, W.A., Erlich, M., Fontoura, S.A.B., and Sayao, A.S.F., eds., Landslides—Evaluation 

and stabilization, Proceedings of the 9th International Symposium on Landslides: London, A.A. 

Balkema Publishers, v. 1, p. 23–36. 

Srivastava, R., and Yeh, T.-C. J., 1991, Analytical solutions for one-dimensional, transient 

infiltration toward the water table in homogeneous and layered soils: Water Resources Research 

v. 27, p. 753–762. 

Tarboton, D.G., 1997, A new method for the determination of flow directions and contributing 

areas in grid digital elevation models: Water Resources Research, v. 33, no. 2, p. 309–319. 

Taylor, D.W., 1948, Fundamentals of Soil Mechanics: New York, Wiley, 700 p. 

Vanapalli, S. K., and Fredlund, D.G., 2000, Comparison of different procedures to predict 

unsaturated soil shear strength, in Shackelford, C.D., Houston, S.L., and Chang N.-Y., eds., 

Advances in unsaturated geotechnics. (Proceedings of Geo-Denver 2000, Denver, Colo., August 

5-8, 2000): Reston, Va., American Society of Civil Engineers, Geotechnical Special Publication 

99, p. 195–209. 
δ

z

x y

Figure 1. Sketch showing local xyz coordinate system with x tangent to the local surface slope, y

tangent to the local topographic contour and z normal to the x-y plane; the slope angle, δ , is with 

respect to the horizontal. 

28
Ground surface

δ

z cos(δ)

x sin(δ)

z

Z(x,z)

x

Figure 2. Sketch showing the definition of the vertical coordinate, = δ + zxZ cossin δ , used to 

compute the depth at an arbitrary location. If Z and z have a common origin, then Z

(after Iverson, 2000). 

= z cos/ δ

29
Ground surface Ground surface

dLZ

Water table Water table

A. B.

Uniform KZ

lower KZ

higher KZ

dLZ

d

d

Figure 3. Schematic cross sections showing boundary conditions assumed in deriving the two 

infiltration formulas (equations 1 and 2). The dimension Zmax corresponds to the input variable 

zmax, and the initial depth to the water table, d , corresponds to the input variable cdep (table 3). A, 

The infinite-depth boundary condition, which is assumed in equation 1, is appropriate for areas 

where the vertical hydraulic conductivity is relatively uniform with depth, B, The finite-depth 

boundary condition, which is assumed in equation 2, is appropriate where a more permeable 

surficial layer overlies a less permeable substrate, such as regolith over bedrock. 

30
t

1 t3 t2 ti ti+1 tN tN+1

Time

IZi /KZ

IZ1 /KZ

IZ2 /KZ

IZi /KZ

IZN /KZ

Figure 4. Hypothetical rainfall history of varying duration and intensity. In practice, the initial time, 

is set to zero. Times t , ..., t are stored in the array capt() (see tables 2 and 3). The first 

rainfall period begins at time and ends at time t ; the i

th rainfall period begins at time t and ends 

at time t . Rainfall intensities for each period are stored either in the array cri(), if constant over 

the study area, or in the files rifil(1), ..., rifil(nper) if intensity during each period varies over the 

study area (table 2). 

1t 1 N+1

1t 2 i

i+1

31
0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

-0.8 -0.4 0 0.4 0.8 1.2

DEPTH BELOW 

GROUND 

SURFA

CE, IN METER

S

PRESSURE HEAD, IN METERS

t=0 s

t=50 s

t=200 s

t=400 s

Final water table

Initial

water

table

Figure 5. Example of transient water-table rise and pressure-head changes in response to 

precipitation at a cell where steady upward flow occurs but the water table is initially beneath the 

ground surface. The TRIGRS program maintains the initial upward hydraulic gradient below the 

water table and adds the transient infiltration solution to the initial condition above the water table. 

In the example shown here: KS=1 x 10-4 m/s, =1 x 10-3 m2

/s, δ = 31°, =-1.0 x 10-5, =5 x 10-5 m/s 

(see equations 1A, 2, and 3B). 

D0 I I ZLT nZ

32
Capillary

fringe

Water table

δ

du

d

Z

max

Unsaturated

 zone

Saturated

zone

Ground 

surface

Basal

boundary

Figure 6. Shallow ground-water conditions in hillside soils. The unsaturated zone above the water 

table has depth . The capillary fringe is between the unsaturated zone and the water table at 

depth, 

du

d. The lower boundary, which is treated as impervious in this model, is at depth Z . max

33
Capillary Fringe

Fillable

pore space

Depth, 

Z

Water content, θ

θr θs

Water table

Unsaturated Saturated

Water-filled

pore space Capillary Fringe

Fillable

pore space

Depth, 

Z

Water content, θ

θr θs

Water table

Unsaturated Saturated

Water-filled

pore space

Capillary Fringe

Fillable

pore space

Depth, 

Z

Water content, θ

θr θs

Water table

Unsaturated Saturated

VA

Zk

Zk+1

t0 t1

t0 t0 t1

hn

Δh

A B C

d

hn

d dt

t0

t0

t0

t0

t1

t1

34

Figure 7. Vertical soil water-content profiles from the ground surface down to an arbitrary 

depth, Z, showing the water-table rise as a result of infiltration. A, Initial conditions at time 

t0; soil water content increases with depth above the top of the capillary fringe (assumed 

to be saturated in our model) as indicated by width of shaded area. B, Conditions at a 

later time, t

1, showing water-table rise of amount ∆h. The shaded area between dotted 

line t0 and solid line t

1 indicates the water absorbed by the soil. C, Approximations used 

in estimating water-table rise at the later time. The shaded area between dashed line t0

and solid line t

1 represents the computed volume of water absorbed by the soil (change in 

water content from t0 to t1 , see equations 4B, 4C, and 8B). The triangular area labeled 

VA represents the cumulative volume of water (minus drainage, equations 10B, 10C) 

computed to reach the water table or base of the unsaturated zone between t0 and t

1. 

The volume of fillable pore space above the top of the capillary fringe (or above the water 

table if 0 = 0 ) is computed at evenly spaced depths Z1,Z2,...,Zk,...(equations 10–12). The 

water table rise, ∆h , is estimated by linear interpolation on Vf between depths Zk and 

Zk+1 to estimate the height filled by volume VA .
0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0

200400600800

10001200140016001800

NORMALIZED FLOW RATE, q/Ks

TIME, IN SECONDS

Surface influx

Basal outflux

Figure 8. Plot of example input and output of unsaturated zone infiltration model. Partial absorption 

of water infiltrating at the surface results in damping and smoothing of the signal at the base of the 

unsaturated zone. 

35
0.0 -0.5 -1.0 -1.5 -2.0

PRESSURE HEAD, IN METERS

0.0

0.1

0.2

0.3

0.4

0.5

VOLUMETRI

C WATER 

CONTENT

Gardner (1958)

data

Figure 9. Graph showing example of Gardner (1958) models for soil water characteristics. The 

model was fit to soil-water retention data for sandy colluvium derived from glacial outwash sand 

near Everett, Wash. Gardner parameters: α = 3.5, θS = 0.41, θ = 0.05. R

36
1 2

7

6543

111098

13

16

12 14

17

15

18 19

222120

N

Figure 10. TRIGRS grid cell numbering scheme. Numbering begins at upper left (northwest) corner 

of a grid and increases from left to right, row by row. No-data cells (shaded) are not numbered. 

37
103 102 102 101 100 99 97 95 94 92

100 100 100 99 97 96 94 92 90 88

97 97 96 95 93 92 91 88 85 82

94 93 92 90 89 87 84 82 79 76

90 89 87 86 84 81 78 75 74 72

86 84 83 80 77 75 72 70 69 69

81 79 76 74 72 69 68 67 68 71

76 73 71 69 68 66 67 68 71 73

71 69 68 66 65 66 68 71 74 76

68 66 65 65 66 68 71 74 77 81

A

1 2 3 4 5 9 13 17 20 26

6 7 8 10 14 15 19 25 30 34

11 12 16 18 23 27 28 35 40 46

21 22 24 29 32 36 41 45 51 57

31 33 37 38 42 47 53 60 65 68

39 43 44 50 55 61 70 77 82 79

49 52 59 62 69 78 90 91 88 72

56 66 71 80 89 94 92 86 74 67

76 81 85 95 98 96 84 75 64 58

83 93 100 99 97 87 73 63 54 48

B

Figure 11. Diagram of an elevation grid file, comparing values of (A) elevation and (B) topographic 

index. 

38
TopoIndex direction codes ESRI direction codes

1

1

987

64

2 6 3 32 4

16

842

128

N

Figure 12. Diagram comparing numeric codes for D8 flow directions used by the utility program 

TopoIndex with those used by ESRI Arc/Info. 

39
100

1.0

272625

765

15 16 17

919293

949697

9799

.122

.312

.250 .250 .250

.250

.303

.432

.343.182

.172

.338

.108

.922 .078

power-law weighting factor

(exponent=20)

slope squared weighting factor 

(exponent=2)

slope proportional weighting factor

 (exponent=1)

uniform weighting factor

 (exponent=0)

D-infinity weighting factor

 (exponent=-1)

D8 weighting factor

 (exponent>20)

.688

Grid cell number Grid cell elevation

A B

C D

E F

G H

Figure 13. Example weighting factors for runoff distribution from an arbitrary cell (no. 16 of the 

tutorial grid file). A, Grid cell numbers. B, Elevation at cell 16 (shaded) and its nearest neighbor 

cells. C, Weighting factor for D8 method. D, Weighting factors for the D-infinity method. E–H, Slope 

power-law weighting factors for various exponents. Note that weighting factors of flatter cells 

approach zero as the exponent increases. E, Uniform distribution to all adjacent downslope cells. F, 

Distribution to all adjacent downslope cells, with weighting factor proportional to slope. G, 

Distribution to adjacent downslope cells, with weighting factor proportional to square of slope. H, 

Distribution to adjacent downslope cells, with weighting factor proportional to slope raised to 20th

power. Weighting factors for cells 17 and 25 are zero. 

40
zmax

(or czmax)

zmin

Ground surface

(zmax-zmin)

nzs

nzs

equal 

steps

Grid cell centerline

Figure 14. Vertical spacing for computing infiltration and factor of safety. The user specifies the 

upper limit, zmin, the lower limit, zmax, and the number of steps, nzs, (table 3). To prevent 

division by zero, the upper limit, zmin, should be a small, finite depth below the ground surface. 

The TRIGRS program will compute the pore pressure and factor of safety at equally spaced depths 

between zmin and zmax. 

41
42

Table 1. Weighting factors for runoff distribution. 

[The weighting factors, wij , are a function of the slope, sij , between grid cell i and each neighboring downslope cell, j; 

the number of neighboring grid cells, n, that are downslope from cell i, and an exponent selected by the user,. The 

program automatically rounds all weighting factors to three decimal places and adjusts them so that the sum of 

weighting factors at each cell equals 1. The default pattern is to divert flow only to the adjacent cell on the steepest 

downslope path. Weighting factors computed by TopoIndex are valid only for square cells (distance between rows 

equal to distance between columns)]

Exponent Weighting factor Distribution pattern

 < 0

1 4 4

wi d

      = −      

2 4

wi d

    =    ; if two cells

Narrow; flow is diverted only to the one or two 

neighboring downslope cells, j and k, that straddle 

the steepest downslope path. If the steepest 

downslope path is aligned with a grid direction or 

diagonal, then flow is diverted only to the 

neighboring cell in that direction. The weighting 

factors depend on the angle d  between the D8 

flow direction and the steepest slope. Here the 

exponent,, serves only as a flag. 

wij = 1; if only one cell

 = 0

wij = 1

n

Fan-shaped; uniform distribution among adjacent 

downslope cells regardless of slope. 

0 <  ≤ 20

wij = sij

 sij



j =1

j =n

∑ Fan-shaped; the positive value of  results in 

weighting factors that favor steeper flow paths. The 

pattern becomes more strongly skewed toward the 

steepest downslope path as slope angle increases.

 > 20 wij = 1; cell on steepest path Narrow; flow is diverted only to the neighboring 

cell that is on the steepest downslope path (D8 

method, O'Callaghan and Mark, 1984) wij = 0; other downslope

cells
43

Table 2. Sample listing of the initialization file "tr_in.txt". [An explanatory heading precedes each line of data, and each heading occupies a single line. Some of the longer 

headings appear to occupy two lines in the table below; however they occupy only one line in the actual file. These 

long headings are distinguished by a hanging indent. The total number of lines in the file depends on the number of 

property zones and number of periods of different rainfall intensity. Units in this example, time, s; length, m; force, N 

(kg-m/s2

); pressure and stress, Pa (=N/m2); unit weight, N/m3] 

Name of project (up to 255 characters)

TRIGRS, version 2, Tutorial

imax, row, col, nwf, tx, nmax

100, 10, 10, 154, 1, 30

nzs, mmax, nper, zmin, uww, t, zones

10, -100, 2, 0.001, 9.8e3, 216000, 2

zmax, depth, rizero, Min_Slope_Angle (degrees)

-3.001, -2.4, -1.0e-9, 0.

zone, 1

cohesion,phi, uws, diffus, K-sat, Theta-sat,Theta-res,Alpha

3.5e+03, 35., 2.2e+04, 6.0e-06, 1.0e-07, 0.45, 0.05, -0.5

zone, 2

cohesion,phi, uws, diffus, K-sat, Theta-sat,Theta-res,Alpha

8.0e+03, 31., 2.2e+04, 8.0e-4, 1.0e-04, 0.45, 0.06, -8.

cri(1), cri(2), ..., cri(nper)

3.e-7, 9.e-5 

capt(1), capt(2), ..., capt(n), capt(n+1)

0, 172800, 216000

File name of slope angle grid (slofil) 

Data/tutorial/slope.asc

File name of property zone grid (zonfil)

Data/tutorial/zones.asc

File name of depth grid (zfil) 

Data/tutorial/zmax.asc

File name of initial depth of water table grid (depfil)

Data/tutorial/depthwt.asc

File name of initial infiltration rate grid (rizerofil)

Data/tutorial/rizero.asc

List of file name(s) of rainfall intensity for each period, (rifil()) 

Data/tutorial/ri1.asc

Data/tutorial/ri2.asc

File name of grid of D8 runoff receptor cell numbers (nxtfil)

Data/tutorial/TIdscelGrid_tutorial.asc

File name of list of defining runoff computation order (ndxfil)

Data/tutorial/TIcelindxList_tutorial.txt

File name of list of all runoff receptor cells (dscfil)

Data/tutorial/TIdscelList_tutorial.txt

File name of list of runoff weighting factors (wffil)

Data/tutorial/TIwfactorList_tutorial.txt

Folder where output grid files will be stored (folder)

Data/tutorial/

Identification code to be added to names of output files (suffix)

itx1

Save grid files of runoff? Enter T (.true.) or F (.false.)

T 

Save grid of minimum factor of safety? Enter Enter T (.true.) or F (.false.)

T 

Save grid of depth of minimum factor of safety? Enter Enter T (.true.) or F (.false.)

T 

Save grid of pore pressure at depth of minimum factor of safety? Enter Enter T (.true.) or F 

(.false.)

T 

Save grid files of actual infiltration rate? Enter T (.true.) or F (.false.)

T 

Save grid files of unsaturated zone basal flux? Enter T (.true.) or F (.false.)

F 

Save listing of pressure head and factor of safety ("flag")? (Enter -2 detailed, -1 normal, 0 none)

-2 

Number of times to save output grids

1 

Times of output grids 
44

Table 2. Sample listing of the initialization file "tr_in.txt".—Continued 

216000.

Skip other timesteps? Enter T (.true.) or F (.false.)

F 

Use analytic solution for fillable porosity? Enter T (.true.) or F (.false.)

T 

Estimate positive pressure head in rising water table zone (i.e. in lower part of unsat zone)? 

Enter T (.true.) or F (.false.)

T 

Use psi0=-1/alpha? Enter T (.true.) or F (.false.) (False selects the default value, psi0=0)

F 

Log mass balance results? Enter T (.true.) or F (.false.)

T 

Flow direction (enter "gener", "slope", or "hydro")

slope

Add steady background flux to transient infiltration rate to prevent drying beyond the initial 

conditions during periods of zero infiltration?

F 
45

Table 3. Definitions of input variables, excluding file path names, found in the initialization 

file and related input files. 

[See table 2 for location of each variable in the initialization file, "tr_in.txt". Integers are assumed to be 4-byte 

integers]

Program Control

Variable 

name

Type Description Range of values

imax integer Number of grid cells (excluding no-data 

cells). 

imax>0

row integer Number of grid rows. 0<row<imax

col integer Number of grid columns. 0<col<imax

nwf integer Length of arrays to store flow routing 

locations and weighting factors. Set equal to 

total number of runoff weighting factors 

computed by TopoIndex. 

imax≤nwf≤8*imax

tx integer A multiplier that, in combination with nper,

determines how many time steps are used in 

the computations. Computational effort and 

storage requirements increase directly with 

tx.

tx≥1

nmax integer Maximum number of roots, Λn , (equation 

6) and corresponding terms in series 

solutions for unsaturated zone. 

nmax≥1; the series normally 

converges with a few tens of 

terms or less.

nzs integer Number of vertical increments. nzs≥1 

mmax integer Maximum number of terms in series solution 

for impermeable base at finite depth, 

maximum value of m in equation 2; also a 

flag to choose between a finite (eq. 2) or 

infinite depth (eq. 1) solution. 

mmax≥1; typically 

10<mmax<25; setting 

mmax<1 chooses the original 

Iverson solution (equation 1) 

for infinite depth 

zones integer Number of physical property zones used to 

model the problem domain. An integer 

property zone grid must be supplied if more 

than one zone is used.

zones≥1 

flag integer Output flag used to specify how depth 

profiles are saved.

flag = -1 or -2; output depth 

profile at all cells, any other 

value, no profile output

slomin real Minimum slope angle; no computations of 

pore pressure or factor of safety on slopes 

lower than slomin.

slomin≥0 

flowdir character, 

(length=5)

Assumed flow direction for computing slope 

of the beta line, which determines the 

maximum pore pressure.

Also determines maximum drainage rate 

computed by equations 10A–10C. 

flowdir=“gener,” general, 

maximum pressure head 

determined by equation 3A; 

or

“slope,” slope parallel, 

maximum pressure head 

equals that for slope-parallel 

flow; or 

“hydro,” hydrostatic 

conditions, maximum 

pressure head equals depth
46

Table 3. Definitions of input variables, excluding file path names, found in the initialization 

file and related input files.—Continued

Mechanical and Hydrologic Properties and Dimensions

Variable 

name

Type Description Range of values

zmin real Shallowest depth at which to compute pore 

pressures, Zmin

zmin>0; typically 

.01≤zmin≤0.1

czmax real Maximum depth at with to compute pore 

pressures and factor of safety, Zmax

czmax> zmin; read values 

from grid file if negative

cdep real Initial depth of water table, d 0<cdep≤zmax; read values 

from grid file if negative

crizero real Steady, pre-storm infiltration rate, ZLT I crizero>0; read values from 

grid file if negative

uww real Unit weight of water,  w uww>0

c() real array Soil cohesion, c c≥0 

phi() real array Angle of internal friction,  phi>0; read values from grid 

file if negative

uws() real array Total unit weight of soil,  s uws>0; read values from grid 

file if negative

dif() real array Hydraulic diffusivity, D0 dif>zmin; read values from 

grid file if negative

ks() real array Vertical hydraulic conductivity of saturated 

soil, KS

ks>0; read values from grid 

file if negative

ths() real array Soil saturated volumetric water content, s ths>0

thr() real array Soil residual volumetric water content, r ths>thr≥0 

alp() real array A fitting parameter for soil size distribution, 

 , approximately equal to the inverse of the 

height of capillary rise.

alp>0; negative value causes 

TRIGRS to treat the zone as 

tension saturated
47

Table 3. Definitions of input variables, excluding file path names, found in the initialization 

file and related input files.—Continued

Storm Properties

Variable 

name

Type Description Range of values

nper integer Number of periods of different rainfall rates 

during a storm, N

nper≥1 

t real Elapsed time since start of storm, t t>0

cri() real array Rate of application (rainfall intensity) during 

each storm increment, InZ

cri()>0; read values from grid 

file if negative

capt() real array Cumulative duration of storm at end of each 

increment, tn

capt(1)=0; capt(1)<capt(2)<, 

..., capt(n)<capt(n+1)

Runoff-Routing Parameters (found only in data files, not in initialization file)

Variable 

name

Type Description Range of values

nxt() integer 

array

Cell that is downslope along the steepest path 

from the present cell

1≤ nxt ≤imax

indx() integer 

array

Index of cell numbers from (topographically) 

highest to lowest

1≤ indx ≤imax

dsctr() integer 

array

Counter array to keep track of downslope 

cells

dsctr≤nwf

dsc() integer 

array

Cell numbers of downslope cells dsc≤nwf

wf() real array Weighting factors for distributing runoff to 

downslope cells

0≤wf≤1 

Map Grid Data

Variable 

name

Type Description Range of values

zmax() real array Depth of basal boundary for each cell zmax>0

slo() real array Slope angles slo≥0 

depth() real array Initial depth of water table depth≥0 

rizero() real array Steady, background infiltration rate All real values

ri() real array Storm infiltration rate ri>0

zo() integer 

array

Property zone numbers zo>0
48

Table 4. Output files generated by TRIGRS. [The abbreviation "TR" at the beginning of a file name identifies a file generated by the program TRIGRS; 

XXXXXXXX denotes the user-defined identification code (as many as eight characters) used to group files from the 

same run. Numerals in parentheses indicate the order in which optional files are selected in the initialization file (table 

2)]

File Name Description

TrigrsLog.txt Log file; TRIGRS overwrites this file every time it runs. 

TRrunoffPer###XXXXXXXX.asc Grid of computed runoff during a given storm period; ### 

designates the storm period ending at time capt(###) (1). 

TRfs_min_XXXXXXXX.asc Grid of the minimum computed factor of safety at each cell 

for time, t (2). 

TRz_at_fs_min_XXXXXXXX.asc Grid of depths at which the minimum factor of safety was 

computed for time t (3). 

TRp_at_fs_min_XXXXXXXX.asc Grid of pressure head at the depth corresponding to minimum 

factor of safety for time t (4). 

TRinfilratPer###XXXXXXXX.asc Grid of actual infiltration rates during a given storm period; 

### designates the storm period ending at time capt(###) (5). 

TRunszfluxTS###XXXXXXXX.asc Grid or unsaturated zone basal flux for a given storm period; 

### designates the storm period ending at time capt(###) (6). 

TRlist_z_p_fs_XXXXXXXX.txt Cell by cell list of pressure head and factor of safety at each 

depth increment (see table 9). Use output flag to select 

abbreviated, expanded, or no output (7). 

TRnon_convrg_UZ_XXXXXXXX.asc

TRnon_convrg_SZ_XXXXXXXX.asc

Grid files locating non-convergent cells. Generated 

automatically if infinite series computations at any cells do 

not converge. UZ denotes unsaturated zone, SZ denotes 

saturated zone.
49

Table 5. Example listings of a grid file for a 10-row by 10-column grid.

A, Heading is for ESRI’s ASCII grid format. 

ncols 10

nrows 10

xllcorner 563435

yllcorner 5258305

cellsize 10

NODATA_value -9999

103 102 102 101 100 99 97 95 94 92 

100 100 100 99 97 96 94 92 90 88 

97 97 96 95 93 92 91 88 85 82 

94 93 92 90 89 87 84 82 79 76 

90 89 87 86 84 81 78 75 74 72 

86 84 83 80 77 75 72 70 69 69 

81 79 76 74 72 69 68 67 68 71 

76 73 71 69 68 66 67 68 71 73 

71 69 68 66 65 66 68 71 74 76 

68 66 65 65 66 68 71 74 77 81 
50

Table 5. Example listings of a grid file for a 10-row by 10-column grid.—Continued

B, Heading is for GRASS GIS two-dimensional raster ASCII format. 

north: 5258405

south: 5258305

east: 563535

west: 563435

rows: 10

cols: 10

103 102 102 101 100 99 97 95 94 92 

100 100 100 99 97 96 94 92 90 88 

97 97 96 95 93 92 91 88 85 82 

94 93 92 90 89 87 84 82 79 76 

90 89 87 86 84 81 78 75 74 72 

86 84 83 80 77 75 72 70 69 69 

81 79 76 74 72 69 68 67 68 71 

76 73 71 69 68 66 67 68 71 73 

71 69 68 66 65 66 68 71 74 76 

68 66 65 65 66 68 71 74 77 81 
51

Table 6. Partial listing of the file "TIdscelList_tutorial.txt", which contains the cell numbers 

of downslope receptor cells. 

[The number -9999 indicates the beginning of a new cell, the next line contains the number of the grid cell, and the 

following lines list as many as eight adjacent cells that are downslope from the grid cell designated at the beginning of 

the list. The sequence repeats for each cell in the grid. If no other cells are downslope from a particular cell, such as 

cells 93 and 94, its cell number is indicated as the only downslope cell. Text in the column to the right of the numbers 

identifies each line of data, but does not appear in the actual output files]

-9999 separator

 1 cell number

 11 neighboring downslope cell number

 -9999 separator

 2 cell number

 12 neighboring downslope cell number

 -9999 separator

 3 cell number

 13 neighboring downslope cell number

 14 neighboring downslope cell number

 -9999 separator

 4 cell number

 15 neighboring downslope cell number

 ... 

-9999 separator

 92 cell number

 93 neighboring downslope cell number

 -9999 separator

 93 cell number

 93 neighboring downslope cell number

 -9999 separator

 94 cell number

 94 neighboring downslope cell number

 -9999 separator

 95 cell number

 94 neighboring downslope cell number

 ... 
52

Table 7. Partial listing of the file " TIwfactorList_tutorial.txt ", which contains the weighting 

factors of downslope cells. In this example, the weighting factors are based on the Dinfinity method (Tarboton, 1997). 

[The number -9999. indicates the beginning of a new cell, and the separator is followed by a decimal point, because the 

weighting factors are floating point values. The next line contains the number of the grid cell, and the following lines 

list weighting factors corresponding to cells listed in table 6. The sequence repeats for each cell in the grid. If no other 

cells are downslope from a particular cell, such as cells 93 and 94, it receives a weighting factor of 1.000. Text in the 

column to the right of the numbers identifies each line of data, but does not appear in the actual output files]

-9999. separator

 1 cell number

1.000 weighting factor for neighboring downslope cell

-9999. separator

 2 cell number

1.000 weighting factor for neighboring downslope cell

-9999. separator

 3 cell number

0.410 weighting factor for neighboring downslope cell

0.590 weighting factor for neighboring downslope cell

-9999. separator

 4 cell number

1.000 weighting factor for neighboring downslope cell

 ... 

-9999. separator

 92 cell number

1.000 weighting factor for neighboring downslope cell

-9999. separator

 93 cell number

1.000 weighting factor for neighboring downslope cell

-9999. separator

 94 cell number

1.000 weighting factor for neighboring downslope cell

-9999. separator

 95 cell number

1.000 weighting factor for neighboring downslope cell

 ... 
53

Table 8. Partial listing of the file "TIcelindxList_tutorial.txt," which contains the cell index 

for a grid of 100 cells taken from the tutorial data set. 

[The first column indicates the index to each cell number listed in the second column. Cell numbers are assigned 

beginning at the upper left corner of the grid from left to right across rows and from top to bottom (figs. 10 and 11). 

No-data cells are unnumbered. Refer to section of the text titled "Topographic Indexing" for further information] 

 1 1

 2 2

 3 3

 4 4

 5 5

 6 11

 7 12

 8 13

 9 6

 10 14

 11 21

 12 22

 ... 

 88 69

 89 75

 90 67

 91 68

 92 77

 93 92

 94 76

 95 84

 96 86

 97 95

 98 85

 99 94

 100 93
54

Table 9. Partial listing of the TRIGRS output file " TRlist_z_p_fs_XXXXXXXX.txt ", depth 

profile of pressure head and factor of safety for each cell. 

[Values appearing in the actual file may differ from those listed here, depending upon the user's choice of input 

parameters. The file has a three-line header to identify the numeric data that follow. The first line after the header

contains the cell number and the slope angle in degrees; the next six lines contain six columns of fixed-width output. 

The first column lists the depth, Z; the next three list the pore pressure, as pressure head at depth Z (column two, total 

pressure head; column three, initial pressure head; column four, transient pressure head); column 5 lists the steady flow 

pressure cutoff; and the sixth lists the factor of safety. The sequence repeats for the next cell starting on the eighth line

after the header. The number of lines of (six-column) output for each cell is one plus the number of vertical increments 

(1+nzs, see fig. 14). Output shown below results when flag is set to -2; if flag is set to -1, then the third, fourth and 

fifth columns are omitted]

TRIGRS depth profiles at each cell

Cell Number, Slope angle

Z P Pzero Ptran Pbeta FS

1 16.7000

0.10000E-02 0.91741E-03 -0.91650 5.5674 0.91741E-03 10.000

0.48080 0.44109 -0.47632 5.1009 0.44109 0.86418

0.96060 0.88127 -0.36146E-01 4.6608 0.88127 0.84360

 1.4404 1.3214 0.40403 4.2470 1.3214 0.83672

 1.9202 1.7616 0.84420 3.8590 1.7616 0.83329

 2.4000 2.2018 1.2844 3.4965 2.2018 0.83122

 2 11.3000

0.10000E-02 0.96160E-03 -0.96063 5.6999 0.96160E-03 10.000

0.48080 0.46233 -0.49926 5.2331 0.46233 1.2947

0.96060 0.92371 -0.37887E-01 4.7921 0.92371 1.2652

 1.4404 1.3851 0.42349 4.3767 1.3851 1.2553

 1.9202 1.8465 0.88486 3.9867 1.8465 1.2504

 2.4000 2.3078 1.3462 3.6216 2.3078 1.2474

 3 12.0000

0.10000E-02 0.95676E-03 -0.95581 5.6856 0.95676E-03 10.000

0.48080 0.46001 -0.49675 5.2188 0.46001 1.2174

0.96060 0.91907 -0.37696E-01 4.7779 0.91907 1.1895

 1.4404 1.3781 0.42136 4.3627 1.3781 1.1802

 1.9202 1.8372 0.88041 3.9728 1.8372 1.1755

 2.4000 2.2962 1.3395 3.6080 2.2962 1.1728
55

Table 10. Example listing of entries in the initialization file for the program TopoIndex, 

"tpx_in.txt"

[File path names are limited to 255 characters; “Rows”, “Columns”, and “Number of iterations” are positive integers. 

“Exponent” is a real number as defined in table 1. The default folder for storing output is the folder where the program 

TopoIndex resides]

Name of project (up to 255 characters)

Tutorial data for TopoIndex analysis

Rows, Columns, flow-direction numbering scheme (ESRI=1, TopoIndex=2)

10, 10, 2

Exponent, Number of iterations

25, 10

Name of elevation grid file 

data/tutorial/dem.txt

Name of direction grid 

data/tutorial/directions.txt

Save listing of D8 downslope neighbor cells? Enter T (.true.) or F (.false.) 

T 

Save grid of D8 downslope neighbor cells? Enter T (.true.) or F (.false.)

T 

Save cell index number grid? Enter T (.true.) or F (.false.)

T 

Save list of cell number and corresponding index number? Enter T (.true.) or F (.false.) 

T 

Save flow-direction grid remapped from ESRI to TopoIndex? Enter T (.true.) or F (.false.)

T 

Name of folder to store output? Leave blank to use default folder

data/tutorial/

Identification code for output files? (8 characters or less)

tutorial
56

Table 11. Names of TopoIndex output files. 

[The abbreviation "TI" at the beginning of the file name identifies a file generated by the program TopoIndex; 

XXXXXXXX denotes the user-defined identification code used to group files from the same run. Numerals in 

parenthesis indicate the order in which optional files are selected in the initialization file (table 10)]

File Name Description

TopoIndexLog.txt Log file; TopoIndex overwrites this file every time it runs. 

TIdsneiList_XXXXXXXX.txt List of downslope receptor cells in D8 flow direction (1). 

TIdscelGrid_XXXXXXXX.asc *Grid of downslope receptor cells in D8 flow direction (2). 

TIcelindxGrid_XXXXXXXX.asc Grid of cell index numbers that specify computation order for 

runoff routing in TRIGRS (3). 

TIcelindxList_XXXXXXXX.txt *List of cell numbers and corresponding index number 

(computation order) (4). 

TIflodirGrid_XXXXXXXX.asc Flow-direction grid in terms of the TopoIndex numbering 

scheme shown in figure 12 (5). 

TIdscelList_XXXXXXXX.txt *Cell by cell list of downslope cells for which non-zero 

weighting factors have been computed. 

TIwfactorList_XXXXXXXX.txt *Cell by cell list of weighting factors for downslope receptor 

cells. 

*Required by TRIGRS for runoff routing. 
57

Table 12. Listing of GridMatch initialization file, “gm_in.txt.” 

[This listing shows an example of three files to be compared. A line of descriptive text precedes each line of numeric 

data and the list of input files. The first file listed is always treated as the master file and all subsequent files are 

compared to the master file.]

Number of grid files to test

3 

Rows, columns

10,10

Names of input files (text string, up to 255 characters per line; one file per line, list master 

grid first)

data/gridmatch/master.asc

data/gridmatch/test1.asc

data/gridmatch/test2.asc
58

Table 13. Common problems and solutions for users of the TRIGRS and TopoIndex 

programs

Problem Things to check

The program displays the message 

“Cannot locate default initialization 

file, type name of initialization file and 

press RETURN to continue.” 

1. Make sure the initialization file is in the working directory. 

2. If you have been editing the initialization file, save and close the file. 

3. Make sure the initialization file is named correctly.

TopoIndex does not converge or exits 

without generating the expected output 

files.

Check the log file “TopoIndexLog.txt” to see if the corrections converged 

rapidly during the first several iterations and then stuck at a constant number 

of corrections during the remaining iterations until the program terminated. 

If so, then, 

1. Use GIS software to verify that the DEM was hydrologically corrected by 

filling single-cell depressions. (Genuine depressions such as lakes may be 

included in the analysis if adequate bathymetric data exist or they can be 

excluded from the analysis by designating such areas as no-data cells in all 

grids). 

2. Verify that the number and locations of no-data cells in the DEM and flow 

direction grids match. 

3. Verify that similar values are used to represent no-data cells and integer 

and floating-point grids.

4. Verify the format and spelling in the headings (check spelling and 

capitalization of labels) of grid cells.

TRIGRS displays an error message 

referencing a file name. 

1. Check the path name of the file, as entered in the initialization file, to be 

sure it is typed correctly and that the folder exists. 

2. Confirm that the file is not locked or open in another application.

TRIGRS or TopoIndex displays a 

message that says "Error reading file 

filename." 

1. Check file format; grid files should follow the format shown in tables 5A

and 5B. 

2. Check end of file. Open the offending data file and position the cursor at 

the end of the file. If the end of the file is directly after the last value, press 

return/enter to type an “end-of line” character, then save and close the data 

file. (If no return character exists at the end of the last line data, an error 

occurs when the programs try to read the last line). 

3. Check file type (line endings). Different operating systems use different 

line endings for text files. Use a text editor to confirm that the data files are 

of the correct type for your operating system (DOS for PC or UNIX for 

MacOS X, Linux, and other UNIX systems).

TRIGRS displays an error message 

referencing a program line number. 

This may result from errors in reading or writing a file. 

1. Identify and open the corresponding file, correct any errors, and try again. 

2. In case the program reports array boundary errors, check the initialization 

file for incorrect values (tables 2, 3, and 4) and for input or output files listed

in the wrong order. 

3. Check the input grids for mismatched no-data cells. 

The arrangement of cells in the output 

grids does not match that in the 

original DEM. 

1. Check the input grids for mismatched no-data cells. 

2. Check GIS coverages for differences in boundaries used to assign 

properties or to define no-data cells of mismatching grids. For example, if 

using geologic map coverage to generate grids of physical properties and 

applying a shoreline from a different coverage to define no-data cells in the 

elevation grid, a mismatch may occur along the shoreline. In other words, 

one grid may contain no-data cells at the same (row-column) locations as data 

cells in the other grid. 
59

Appendix 1. Listing of Files Included in This Distribution

Folder names appear below in boldface type and slashes indicate folder hierarchy. File 

suffixes for make files and compiled programs vary between operating systems; the “.exe” and 

“.amk” suffixes used below correspond to the Windows convention. 

TRIGRS/

gm_in.txt 

memoryReq.xls 

tpx_in.txt 

tr_in.txt

tr_in_flume.txt

tr_in_mc.txt

tr_in_sy.txt 

uc_in.txt 

TRIGRS/bin/ 

GridMatch.exe 

TopoIndex.exe 

TRIGRS.exe

unitconvert.exe 

TRIGRS/data/flume/ 

slope.asc

TIcelindxList_flume.txt

TIdscelGrid_flume.asc

TIdscelList_flume.txt

TIwfactorList_flume.txt

zones_flume.asc

TRIGRS/data/gridmatch/ 

master.asc

test1.asc

test2.asc

TRIGRS/data/MinorCreek/ 

slope.asc

TIcelindxList_MC.txt 

TIdscelGrid_MC.asc

TIdscelList_MC.txt

TIwfactorList_MC.txt 

zones_MC.asc

TRIGRS/data/sy91/ 

slope.asc

zones.asc
60

TRIGRS/data/tutorial/ 

dem.asc

depthwt.asc 

directions.asc 

ri1.asc 

ri2.asc 

rizero.asc

slope.asc

TIcelindxList_tutorial.txt

TIdscelGrid_tutorial.asc

TIdscelList_tutorial.txt

TIwfactorList_tutorial.txt

zmax.asc

zones.asc

TRIGRS/doc/

Readme.txt 

TRIGRS/src/GridMatch/ 

GridMatch.amk

gridmatch.f90 

srdgrd.f 

TRIGRS/src/TopoIndex/ 

isvgrd.f 

mpfldr.f

nxtcel.f90 

rdflodir.f

sindex.f 

slofac.f

srdgrd1.f 

TopoIndex.amk 

tpindx.f90 

TRIGRS/src/TRIGRS/ 

calerf.f 

dbsct.f 

derfc.f

dsimps.f95 

flux.f90 

grids.f95 

input_file_defs.f95 

input_vars.f95 

irdgrd.f 

irdswm.f

isvgrd.f 

iverson.f95 

ivestp.f95 
61

model_vars.f95 

pstpf.f95 

pstpi.f95 

rnoff.f95

roots.f 

satfin.f95 

satinf.f95

savage.f95 

smallt.f95

srdgrd.f 

srdswm.f 

ssvgrd.f 

steady.f95 

svgstp.f95 

TRIGRS.amk

trigrs.f90

trini.f95

unsfin.f95 

unsinf.f95 

unsth.f95 

TRIGRS/src/UnitConvert/ 

srdgrd.f 

ssvgrd.f 

UnitConvert.amk 

unitconvert.f90 
62

Appendix 2. Tutorial 

We have provided a data set for a 100-cell (10 cells wide by 10 cells long) grid to 

demonstrate use of TRIGRS, along with separate single-cell grids for reproducing Iverson's (2000) 

flume and Minor Creek landslide analyses and Srivastava and Yeh’s (1991) infiltration analysis. 

These data sets are found in separate folders within the folder “data” and include a complete set of 

files for topographic, precipitation, and physical properties data, along with necessary data for 

routing runoff, where applicable (appendix 1). The following paragraphs describe use of the 

programs TopoIndex and TRIGRS and how a user controls an analysis by editing the initialization 

files. Data for this exercise are found in the folder "tutorial," within the "data" folder. By 

following the instructions in this tutorial, the user will analyze a small sample problem that 

illustrates many of the features of TRIGRS. The sample data set includes spatially variable 

geometry and physical properties data, a spring, and two-stage (spatially uniform) rainfall. The 

tutorial analysis illustrates development of local runoff, gradual build-up of pore pressure as a 

result of infiltration, and local slope failure at the end of the storm. Before starting, make sure that 

your system preferences are set to show file extensions1

Using TopoIndex to Prepare Data for Runoff Routing

 and that you have a text editor installed on 

your system. 

Use the program TopoIndex at the beginning of a project to prepare data needed for runoff 

routing in TRIGRS. Once a topographic data set and routing method are selected for a study area, 

the output of TopoIndex can be reused for multiple analyses of that area by TRIGRS. 

First, we will examine the necessary input files for TopoIndex. In the folder TRIGRS, open 

the file “tpx_in.txt”; the contents should resemble table 10A in the main text. Table 2-1 lists the 

contents of the file with line numbers to aid discussion. A line of explanatory text precedes each 

line of data. For example, the text on line 1 describes the data on line 2. The project title appearing 

on line 2 (tables 10 and 2-1) may include letters, numbers, and spaces. The succeeding paragraphs 

describe the other information appearing in the file “tpx_in.txt.” and the names of input files. The 

user may also select certain output files, enter an identification code that will help group the files 

from a given run, and may specify the folder where the files will be stored. 

The first two values on line 4 specify the number of rows and columns in the elevation grid 

file. The number of rows and columns appears in the heading of the grid files. Open the elevation 

grid file "dem.asc" in the folder "TRIGRS/data/tutorial" to verify that the number of rows and 

columns listed in "tpx_in.txt" agrees with corresponding numbers in the grid file (fig. 2-1 and table 

2-1, lines 3–4). If necessary, correct the number of rows and columns in the initialization file and 

then close the grid file. 

The numbering scheme of the (input) direction grid (fig. 12) is the third value listed on line 

4 (table 2-1); the example data uses the TopoIndex numbering scheme. Open the direction grid 

file, "directions.asc" (in the same folder), to determine which grid-direction numbering scheme is 

used in the file. Compare the values listed in the grid to figure 12. If the grid values range from 1 

 

1

If Windows is not set to display file extensions, select “Folder Options” from the “View” menu (Windows XP or the 

“Tools” menu in Windows 2000) at the top of the “TRIGRS” folder window. Click on the “View” tab in the “Folder 

Options” window and uncheck the box next to “Hide file extensions for known file types” and then click “OK” to save 

the change. To show extensions in Mac OS X, click the pointer on a vacant spot on the desktop to make the “Finder” 

active and then select “Preferences” from the “Finder” menu. Check the box next to “Always show file extensions” 

and then click the red button in the upper left corner of the “Finder Preferences” window to save the change and close 

the window.
63

to 9, then the initialization file should contain the numeral 2 at the end of its second line of data to 

indicate that the grid uses the TopoIndex numbering scheme. If the numbers in the grid are powers 

of two ranging from 1 to 128, then the initialization file should contain the numeral 1 at the end of 

the second line of data (the project title constitutes the first line of data) to indicate that the grid 

uses the ESRI numbering scheme. When finished checking the numbering scheme, close the 

direction-grid file and enter the correct value on line 4 of the initialization file. 

The first value on line 6 is an exponent (, table 1) that determines the flow distribution 

weighting factors as described previously (fig. 13). The exponent selects one of the runoff routing 

methods defined in table 1. Figure 13 shows example results for various exponents. Return to the 

initialization file “tpx_in.txt” and look at the exponent, on the third line of data (table 2-1, lines 5–

6). We suggest starting with any negative exponent, which selects the D-infinity method. If 

greater dispersion is needed to distribute water more smoothly across the slope, then choose a small 

positive exponent. Use a large positive exponent to minimize the dispersion. Enter -1 as the initial 

value for the exponent in the initialization file (table 2-1, line 6). 

The second value on line 6 specifies the maximum number of iterations for indexing the 

cells. As we described previously, TopoIndex uses an iterative process to determine the order in 

which runoff calculations should proceed. The maximum number of iterations needed to complete 

the indexing process varies with the size of the grid. Our limited experience with grids of various 

sizes indicates that the following formula can estimate the maximum number of iterations.

iterations ≤1 +

rows + columns

20 

Verify that the maximum number of iterations listed in the initialization file (table 2-1, lines 5–6) is 

at least as large as indicated by the above formula. For a 10x10 grid, allow at least two iterations. 

The purpose in setting the maximum number of iterations is to allow TopoIndex to stop after a 

reasonable time in case the DEM is not hydrologically consistent and indexing cannot converge. 

TopoIndex will stop after the index solution converges or the maximum number of iterations is 

reached. Consequently, in certain cases the user may need to try a larger number of iterations than 

suggested above in order to achieve convergence. 

Text on line 8 and 10 of the initialization file identifies input files. Verify the names of the 

input files that are listed in the initialization file (table 2-1, lines 7–10). The tutorial folder 

("TRIGRS/data/tutorial") should contain files named "dem.asc" and "directions.asc". In the 

initialization file, the name (path) consists of the folder name followed by a separator and the file 

name. The separator is a slash (/) for Mac OSX2

Letters on lines 12, 14, 16, 18, and 20 determine which of the five possible output files are 

saved. The letter "T" stands for "true", and answers as "yes," while "F" stands for "false" and 

answers as "no." Verify that the necessary output files are selected for output. A "T" should 

appear on the line below the question "Save grid of D8 downslope receptor cells?" and the question 

"Save list of cell number and corresponding index number?" (table 2-1, lines 14, 18). These two 

files are needed to run TRIGRS later in the tutorial. For now, an "F" should appear below the other 

three questions that ask about saving files. 

 or a backslash (\) for Windows. The combination 

of folder and file name can be as long as 255 characters. To avoid using folder or path names, the 

user can copy the executable file and initialization file to the folder that contains the input data 

files. Please consult the documentation for your computer's operating system to learn more about 

rules for specifying file path names. 

 

2 TRIGRS and related programs use the POSIX path names on Mac OSX. Note also that the files use UNIX line 

endings, rather than traditional Mac OS line endings.
64

Text on line 22 of the initialization file identifies the output folder. Verify that the 

specified output folder is "data\tutorial" (second to last line of data in the initialization file). The 

folder name should be typed "data\tutorial\" for a Windows PC or "data/tutorial/" for a Mac OS 

computer. Type an identification code of eight characters or less on the last line of the initialization 

file (in this case we have used "tutorial", see line 24, table 2-1). Be sure to use only characters that 

your computer's operating system allows for file and folder naming. 

Close the file "tpx_in.txt" and double-click the program icon for TopoIndex (see 

instructions in the section “Program Operation” if using Mac OSX). The program will run and then 

exit, closing its window. Open the folder “Tutorial” and then open the files 

"TIdscelList_tutorial.txt", "TIwfactorList_tutorial.txt", and "TIcelindxList_tutorial.txt" to examine 

their contents. Partial listings of these files appear in tables 6, 7, and 8. The file 

"TIcelindxList_tutorial.txt" lists cells in their correct order from highest to lowest elevation; the left 

column contains the index, or pointer, and the right column contains the cell number for each index 

value (table 8). Cell numbers increase from left to right, row by row from top to bottom, skipping 

any no-data cells (fig. 10). Cell numbers in the file "TIdscelList_tutorial.txt" (table 6) correspond 

to weighting factors in the file "TIwfactorList_tutorial.txt" (table 7 and fig. 13). When finished 

studying the output, close the files. At this point, either experiment with calculating the weighting 

factors by different methods (by changing the value of the exponent in “tpx_in.txt”) or proceed 

with the tutorial. 

Preparing the TRIGRS Initialization File for an Analysis

We will now perform a sample analysis of slope stability using TRIGRS for the same grid 

we processed using TopoIndex. Return to the folder TRIGRS and double click the icon “tr_in.txt” 

to view its contents; it should resemble table 2. Notice that an explanatory heading precedes each 

line of data, except that the names of the rainfall-intensity files are listed on succeeding lines under 

a single heading. Tables 2 and 3 define the entries in the initialization file and indicate permissible 

ranges of values. Table 4 identifies output files, most of which are optional and can be selected in 

the initialization file. Table 2-2 contains the same information as table 2, with the addition of line 

numbers to aid the reader in locating some of the entries in the table. Note that changing the 

number of property zones or the number of rainfall intensity periods will change the number of 

lines in the initialization file, but the relative locations of the entries remain the same.

Array sizes.—Compare the array parameters on line 4 (second line of input data) of the 

initialization file against values that were computed by the program TopoIndex (figure 2-2 and 

table 2-2). In the TRIGRS folder, open the file TopoIndexLog.txt and find the line of text "Data 

cells, Rows, Columns, Downslope cells" near the bottom of the list. Compare the corresponding 

numerical values in the line following with the four values on the corresponding data in "tr_in.txt" 

(table 2-2, lines 3–4). If any of the numbers differ, change the values in "tr_in.txt" to agree with 

their corresponding values in "TopoIndexLog.txt". The TRIGRS program uses these numbers to 

allocate storage and the program will quit spontaneously or report an array-boundary error if these 

numbers are too small. 

Simulation parameters.—Several other parameters that control the nature of a given 

simulation are included near the beginning of the initialization file (table 2-2, lines 3–6). These 

parameters are defined in table 3. The number of vertical increments (nzs, fig. 14) determines how 

detailed the simulation is with depth. Choose the smallest number of increments that gives 

sufficient vertical detail when analyzing large grids, because computation time increases directly 

with the number of vertical increments. 
65

The sign of mmax (table 2-2, lines 5–6) determines whether pore pressures are computed 

using equation 1 or equation 2 for saturated infiltration. Setting mmax to a negative value selects 

Iverson’s (2000) solution for infinite depth (equation 1) and a positive value selects the solution for 

finite depth (equation 2). A positive value of mmax determines how many terms of the infinite 

series in equation 2 will be used in computing the pore pressures. Our tests indicate that the series 

converges rapidly and that 20 terms yield sufficient accuracy for most purposes. If convergence 

does not occur at one or more cells, then mmax can be increased. A corresponding parameter, 

nmax, (table 2-2, lines 3–4) controls series solutions for unsaturated infiltration. The parameter 

nmax controls the number of roots found for equation 6 and the maximum number of terms used in 

computing any infinite series where the roots appear. A value of nmax=30 is usually more than 

adequate to insure convergence the series. 

Both saturated and unsaturated formulations are able to compute pore pressures for multistage storms. The rainfall periods, nper, correspond to individual bars in the rainfall histogram of 

figure 4. Two values of rainfall intensity are listed on line 16 of the initialization file (table 2-2), 

corresponding to the two periods in this example (nper=2). The rainfall intensities are in meters per 

second. More details are provided in a later paragraph titled "Rainfall Intensities."

Set the minimum depth, zmin, to a small positive value (0.001 m in this example, table 2) to 

investigate transient infiltration and slope-stability from just below the surface downward to the 

maximum depth. Enter a value for the unit weight of water, uww, to establish the units of 

measurement for your other data. In this case, we are using meter-kilogram-second units, and the 

unit weight of water is approximately 9800 kg-m/s2

 m3(or N/m3). Finally, enter the elapsed time, t, 

since the beginning of the storm. We have entered 216,000 s to compute conditions at the end of 

the storm. Enter the time as "216000" without the comma, because commas are used to delimit 

values in the input file. 

The value tx (lines 3–4) determines how many intermediate time steps are used in the 

calculations. When using the saturated infiltration solutions (equations 1 and 2), tx may be set to 1 

to allow TRIGRS to compute a solution as rapidly as possible for a particular elapsed time. To 

produce a time series, tx can be set to a larger number to achieve the desired number of time steps; 

however, increasing the number of time steps does not affect accuracy of the solution. The number 

of time steps is the product of tx and nper. When using unsaturated infiltration solutions, tx can be 

set to a small positive number (usually 10 or less) to compute pore pressure response. Increasing 

the time steps for unsaturated infiltration may increase accuracy of the results because numerical 

integration is employed in solving for the water-table rise. 

Physical parameters.—Many of the input parameters for geometry, initial conditions, 

physical properties, and precipitation input can be set to constants or read from grids. Maximum 

depth, Zmax , water table depth, d , and initial or background infiltration rate, IZ , are specified in 

lines 7 and 8 (table 2-2); negative values indicate that these will be read from grids (lines 23–28, 

table 2-2) , rather than using spatially constant values. The TRIGRS program now uses property 

zones to assign physical properties, so that a single grid file (lines 21–22) can be used for all 

property values. Lists of properties for each zone are contained in the initialization file (table 2-2, 

lines 9–14). The parameters must be listed in the following order: c ,  ,  s , D0 , KZ ,s , r , and 

 . Enter a positive value for each parameter in each property zone, the property zone grid 

determines how the properties vary spatially by mapping the parameters to their respective zones. 

If properties are uniform throughout the problem domain, then specify only one property zone. No 

property zone grid will be needed and the word “none” can be entered as a placeholder for the file 

name (table 2-2, lines 21–22). Use a positive value of alpha to use the unsaturated infiltration 
66

solution; a negative value will cause TRIGRS to use the saturated infiltration model for that 

property zone. 

Rainfall intensities.—The line of data immediately after the property zone data lists rainfall 

intensities corresponding to the number of periods (nper, table 2-2, line 6) of distinct rainfall 

intensity (table 2-2, lines15–16). The intensities should be listed in their correct temporal order. 

An array of spatially constant rainfall values, cri(:), (lines 15–16 of table 2-2) or grids of spatially 

variable rainfall must exist for each period (table 2-2, lines 29–31). To use spatially constant 

rainfall values, list nper positive values on line 16 in the initialization file and placeholders for file 

names, one per line, for a total of nper lines. To use spatially variable rainfall values, list nper 

negative values on line 16 in the initialization file and file names for the rainfall grids, one per line, 

for a total of nper lines. As with the background flux, entering a negative intensity, cri(), causes 

TRIGRS to read rainfall intensity from a grid file for the corresponding time period (table 2-2, lines 

29–31). 

Time increments.—The next line of numerical data in the initialization file (table 2-2, lines 

17–18) contains times that define the rainfall periods as defined in figure 4. These should be listed 

in their correct temporal order so that the numbers increase from left to right. Three times are listed 

in "tr_in.txt" for this example (table 2-2). The first time, 0 s, is at the beginning of the rainfall; the 

remaining times correspond to the end of each rainfall period and beginning of the next. Thus, 

period one, which has an intensity of 6 x 10-7 m/s (2.16 mm/h) ends after 172,800 s (48 h) and 

period two begins. After 216,000 s, (60 h) period two, which has an intensity of 2x10-5 m/s, ends. 

Thus period two lasts 43,200 s (216,000–172,800) or 12 h. Times and rainfall intensities used in 

this example were selected arbitrarily, but are within the observed ranges of rainfall duration and 

intensity. 

Input files.—The next group of lines lists the names of the input files in the order in which 

TRIGRS expects them (table 2-2, lines 19–39). All the files are in the tutorial data folder, so the 

folder name precedes the name of each file as described previously. Verify that the folder name 

agrees with the location of the input files. In addition to a file corresponding to each physical 

parameter listed on line 8 (maximum depth, initial water table depth, and steady pre-storm 

infiltration rate), the listing includes files for slope angle, rainfall intensity, and runoff routing data. 

Note that the number of files for rainfall intensity must equal the number of storm periods and that 

all are grouped under one heading. Two rainfall intensity files are listed in this example. Supply a 

placeholder, such as "none" for the names of any file names that do not exist. For example, if you 

assume that the initial water-table depth is constant over the entire area, you would not need to 

create a grid file for water-table depth. In that case, you would enter a positive constant for the 

water table depth (table 2-2, line 7–8) and type "none" below the heading that asks for the name of 

the water table grid file (table 2-2, lines 25–26). 

Output files.—The log file, "TrigrsLog.txt", is generated automatically, all others are 

optional. Table 2-2 shows the group of lines in the initialization file that control output (lines 40–

63). If desired, you may supply the name of a folder where the output files will be stored. In this 

case, TRIGRS will write the files in the folder named "data/tutorial/" (“data\tutorial” for the 

Windows version) within the TRIGRS folder (lines 40–41). Leave the line blank if you wish to put 

the files in the same folder where the TRIGRS executable program file resides (or present working 

directory on Mac OSX or UNIX). If desired, provide an identification code (as many as 8 

characters) for output files, otherwise leave a blank line (table 2-2, lines 42–43). For this example, 

the code "t=216000" corresponds the time since rainfall began, specified earlier. Be sure to change 

the code before running the next simulation. Generation of individual output grid files can be 

turned on or off by typing "T" or "F" on the corresponding lines of the initialization file (table 2-2, 

lines 44–55). Table 4 of the main text describes the various output files. An output flag in the 
67

initialization file controls output of a detailed listing of pressure head and factor of safety by depth 

for each grid cell (table 2-2, lines 56–57). For this tutorial, make sure that all of the grid files are 

turned on (T) and that the output flag is set to -2 to generate a full list, like that in table 9. When 

working with large grids, it is best to suppress output of the detailed listing by setting the flag to 0 

(zero) to save disk space.

When using the unsaturated infiltration model or multiple time steps for the saturated 

infiltration model (tx>1, table 2-2, lines 3–4), it is possible to save output grids at more than one 

elapsed time during the storm. Specify the number of times for saving the grids (table 2-2, lines 

58–59), the specific times for output (lines 60–61) and whether to skip the intervening time steps 

(lines 62–63). If the specified output times do not correspond to increments of the time step size, 

then TRIGRS will automatically adjust the output time to the nearest time step. 

Other user options.—A few other user options exist in the current version of TRIGRS (table 

2-2, lines 64–73. The option to use an analytic solution for fillable porosity (lines 64–65) allows 

the user to choose the efficient analytic solution which converges rapidly for intermediate and later 

times or a numerical integration for very early times. The option of estimating positive pressure 

head in the rising water table zone (lines 66–67) allows the user to save to disk either the positive 

pressure estimate below the rising water table (enter “T”) or the negative pore pressure 

computations above the initial water table (enter “F”). Negative pressure computations can be used 

to compare model results with published numerical solutions (Srivastava and Yeh, 1991). The next 

option (table 2-2, lines 68–69) allows TRIGRS to use the unsaturated computations either in the 

entire zone above the initial water table, 0 0 = , or only above the capillary fringe 0 −= 1 

(Savage and others, 2004). Using the second option (only above the capillary fringe) tends to 

produce more rapid water-table rise, because the unsaturated infiltration models operate on a 

thinner layer, but it may produce a better approximation of the soil-water characteristic curve. 

The option of logging mass-balance results (table 2-2, lines 70–71) allows the user to check 

behavior of the unsaturated zone solution for particular combinations of parameters. This option is 

mainly for use in the early stages of modeling on trial cells. Logging mass balance during 

simulations on large grids will slow down computations and greatly increase disk storage 

requirements.

The flow direction option (table 2-2, lines 72–73) allows the user to override the default 

slope of the line that determines the maximum allowable values of computed pore pressures given 

in equation 3A (see also Iverson, 2000). Table 3 of the main text describes the choices and their 

results under the program control variable “flowdir.”

Using TRIGRS to Analyze Slope Stability during a Simple Storm

After studying and editing the initialization file, close the file and double-click the TRIGRS 

application icon to perform an analysis. (If using Mac OSX, follow instructions in the Program 

Operation section of the main text.) The TRIGRS program will attempt to open an initialization 

file named “tr_in.txt” in the same folder as the TRIGRS application resides. If TRIGRS does not 

find the file, it will prompt the user to enter the name of the initialization file. After finding the 

initialization file, TRIGRS should run without any additional intervention by the user. After 

TRIGRS completes execution, you can open the output files in the folder “data/tutorial” to see the 

results. In this analysis, time was set to compute conditions at the end of the storm. To save these 

results before conducting another analysis, add the identification code to the name of the log file 

(change it to "TrigrsLogt=216000.txt") and move it to the folder with the output data so you have a 

record of what input parameters were used. Next, edit the initialization file to change the time, 

identification code, and any other parameters before running TRIGRS again. The TRIGRS 
68

program does not have any display capabilities. To observe your results graphically, import the 

ASCII grid output files into graphics or GIS software.

Time steps.—A typical analysis will involve running TRIGRS for several different times 

throughout the duration of a storm to observe how pore pressure and slope stability vary as the 

storm progresses. To observe pore pressures at an earlier time, open “tr_in.txt” and edit the time 

(table 2-2, lines 5–6). For example, change the elapsed time from 216000 s to 43200 s and run 

TRIGRS again to observe the pore pressure and factor of safety at the middle of the first storm 

period. Then repeat for 86400 s and 100800 s to compute pressure head at the end of the first and 

second storm periods. After each analysis, be sure to save the log file by renaming and moving it 

as described above. Alternately, one can leave the value of t=216000 and set the value of tx to 4 

(table 2-2 lines 3–4), to get output at even time increments. Using this option requires the user to 

specify output times (table 2-2, lines 58–63). 

Figure 2-3 shows graphical representations of the input data and the results of the analysis 

for the hypothetical study area that is represented by the tutorial data grids. 

Input data.—The digital elevation model (DEM, fig. 2-3) represents part of a small canyon 

that drains to the southwest. Geologic material that underlies the upper left part of the area has 

higher shear strength parameters (cohesion and phi, fig. 2-3), and lower hydraulic conductivity, KS , 

and diffusivity, D0 , than material in the remainder of the area. The depth to firm bedrock is 

assumed to be 2 m everywhere except on slopes that are steeper than 30°, where depth decreases 

linearly with increasing slope angle. The precipitation data represent a hypothetical 28-hour storm 

that consists of 24 hours of steady, light precipitation followed by four hours of heavy precipitation. 

We have assumed uniform precipitation over this small area. The steady infiltration rate (not 

shown) is also uniform across the area, except at a spring located in the channel bottom. A 

negative infiltration rate in the seventh column and eighth row of the grid file “rizero.asc,” which is 

found in the folder "TRIGRS/data/tutorial/" represents discharge from the spring.

Results. –Three series of plots in figure 2-3 show the runoff and progressive changes in 

pressure head, , and factor of safety, FS . Runoff occurred during both phases of the storm 

(RUNOFF 1 and 2 in fig. 2-3). The precipitation rate is higher than the hydraulic conductivity of 

material in the upper left so that water runs off this area and infiltrates into the more permeable 

material down slope, as reflected by the higher infiltration rates directly downslope from the 

contact between materials of higher and lower conductivity (INFILTRATION 1 and 2, Fig. 2-3). 

As a result, pressure head increased more rapidly downslope from the contact between the two 

materials. A spring exists at the channel bottom in the lower left corner of the area, which also 

produces runoff during both storm periods. Pressure head increased very little in the less 

permeable material, in the upper left corner of the area, but increased steadily throughout the storm 

in the more permeable material. Likewise, factor of safety changed gradually during the storm, and 

many of the steep cells become unstable by the end of the storm (Fig. 2-3). 

Starting a new problem.—Having completed the tutorial, you should now be ready to set up 

and solve a new problem. Use GIS software to create the input grids and save them in the required 

ASCII format. Create a new folder for all the input data (preferably in the TRIGRS “Data” folder) 

and copy all of the data files for the new problem into the new folder. Use the utility GridMatch to 

make sure that all the grid files are compatible before proceeding. If you plan to use the simple 

runoff routing features in TRIGRS, then use TopoIndex to generate the necessary files. 

Once all the grid and routing files are ready, make a copy of the TRIGRS initialization file 

and rename both the original file and the copy. If any of the initialization files in the TRIGRS 

folder are named “tr_in.txt,” be sure that it is the active (new) one. However, to avoid confusion, it 

might be best to give all the initialization files descriptive names and enter the name of the one you 
69

want to use at run time. After creating and naming a new initialization file, open it in a text editor 

and enter the correct folder and file names for your new problem into the new initialization file. 

Also, be sure to edit the project title and all the other input parameters to suit the new problem. 

Using the unsaturated infiltration model.—Once you are confident using the model to 

simulate saturated infiltration, you are ready to begin experimenting with the unsaturated 

infiltration model. Successful use of the infiltration model requires initial testing of parameters on 

single-cell grids to learn how the model will behave with particular combinations of soil properties 

and model or program control parameters. Results of these tests will guide selection of program 

control parameters, such as nmax, nzs, and tx, as well as ranges of physical properties to achieve 

accurate results with efficient computation. 
70

Figure 2-1. Screen image of the file "dem.asc" with grid dimensions highlighted.
71

Figure 2-2. Screen image of the file "TopoIndexLog.txt" with array length parameters 

highlighted.
65 103m

DEM SLOPE

0 35°

ZMAX

1.5 1.75 1.9 2.0m

COHESION

3.5 8.0kPa

Ks

1x10-7 1x10-4m/s 6x10-6 8x10-4m2/s

PHI Ds

35 31°

RUNOFF 1

0 2.8 4.4 6.0 8.0x10-7m/s 0 0.5 0.7 1.0 1.4x10-5m/s

RUNOFF 2

ψ t = 0.01s ψ t = 43200s ψ t = 86400s ψ t = 100800s

<0.3 0.5 1.0 2.0m <0.3 0.5 1.0 2.0m <0.3 0.5 1.0 2.0m <0.3 0.5 1.0 2.0m

FS t = 0.01s FS t = 43200s FS t = 86400s FS t = 100800s

>1.5 1.5 1.2 <1.0 >1.5 1.5 1.2 <1.0 >1.5 1.5 1.2 <1.0 >1.5 1.5 1.2 <1.0

INFILTRATION 1

0 0.3 0.5 0.7 1.1x10-6m/s

INFILTRATION 2

0 0.4 0.7 1.3 1.7x10-5m/s

Figure 2-3. Sample analysis of data from the folder “TRIGRS/data/tutorial”. The first two rows 

represent input data (elevation, slope angle, depth of soil, soil cohesion, internal angle of friction, 

hydraulic conductivity, and hydraulic diffusivity). The third row shows runoff and infiltration rates 

during the two storm periods. The fourth row shows pressure head at the base of the soil through 

time, and the bottom row shows the factor of safety through time. 

72
73

Table 2-1. Contents of the initialization file “tpx_in.txt” for tutorial problem, listed by line 

number. [File path names are limited to 255 characters; “Rows”, “Columns”, and “Number of iterations” are positive integers. 

“Exponent” is a real number as defined in table 1. The default folder for storing output is the folder where the program 

TopoIndex resides]

Line 

number

Text

1 Name of project (up to 255 characters)

2 Tutorial data for TopoIndex analysis

3 Rows, Columns, flow-direction numbering scheme (ESRI=1, TopoIndex=2)

4 10, 10, 2

5 Exponent, Number of iterations

6 25, 10

7 Name of elevation grid file 

8 data/tutorial/dem.txt

9 Name of direction grid 

10 data/tutorial/directions.txt

11 Save listing of D8 downslope neighbor cells? Enter T (.true.) or F (.false.) 

12 T 

13 Save grid of D8 downslope neighbor cells? Enter T (.true.) or F (.false.)

14 T 

15 Save cell index number grid? Enter T (.true.) or F (.false.)

16 T 

17 Save list of cell number and corresponding index number? Enter T (.true.) or F (.false.) 

18 T 

19 Save flow-direction grid remapped from ESRI to TopoIndex? Enter T (.true.) or F (.false.)

20 T 

21 Name of folder to store output? Leave blank to use default folder

22 data/tutorial/

23 Identification code for output files? (8 characters or less)

24 tutorial
74

Table 2-2. Contents of initialization file “tr_in.txt” for tutorial problem, listed by line 

number. 

Line 

number

Text

1 Name of project (up to 255 characters)

2 TRIGRS, version 2, Tutorial

3 imax, row, col, nwf, tx, nmax

4 100, 10, 10, 154, 1, 30

5 nzs, mmax, nper, zmin, uww, t, zones

6 10, -100, 2, 0.001, 9.8e3, 216000, 2

7 zmax, depth, rizero, Min_Slope_Angle (degrees)

8 -3.001, -2.4, -1.0e-9, 0.

9 zone, 1

10 cohesion,phi, uws, diffus, K-sat, Theta-sat,Theta-res,Alpha

11 3.5e+03, 35., 2.2e+04, 6.0e-06, 1.0e-07, 0.45, 0.05, -0.5

12 zone, 2

13 cohesion,phi, uws, diffus, K-sat, Theta-sat,Theta-res,Alpha

14 8.0e+03, 31., 2.2e+04, 8.0e-4, 1.0e-04, 0.45, 0.06, -8.

15 cri(1), cri(2), ..., cri(nper)

16 3.e-7, 9.e-5

17 capt(1), capt(2), ..., capt(n), capt(n+1)

18 0, 172800, 216000

19 File name of slope angle grid (slofil) 

20 Data/tutorial/slope.asc

21 File name of property zone grid (zonfil)

22 Data/tutorial/zones.asc

23 File name of depth grid (zfil) 

24 Data/tutorial/zmax.asc

25 File name of initial depth of water table grid (depfil)

26 Data/tutorial/depthwt.asc

27 File name of initial infiltration rate grid (rizerofil)

28 Data/tutorial/rizero.asc

29 List of file name(s) of rainfall intensity for each period, (rifil()) 

30 Data/tutorial/ri1.asc

31 Data/tutorial/ri2.asc

32 File name of grid of D8 runoff receptor cell numbers (nxtfil)

33 Data/tutorial/TIdscelGrid_tutorial.asc

34 File name of list of defining runoff computation order (ndxfil)

35 Data/tutorial/TIcelindxList_tutorial.txt

36 File name of list of all runoff receptor cells (dscfil)

37 Data/tutorial/TIdscelList_tutorial.txt

38 File name of list of runoff weighting factors (wffil)

39 Data/tutorial/TIwfactorList_tutorial.txt

40 Folder where output grid files will be stored (folder)

41 Data/tutorial/

42 Identification code to be added to names of output files (suffix)

43 itx1

44 Save grid files of runoff? Enter T (.true.) or F (.false.)

45 T

46 Save grid of minimum factor of safety? Enter T (.true.) or F (.false.)

47 T
75

Table 2-2. Contents of initialization file “tr_in.txt” for tutorial problem, listed by line 

number.—Continued.

Line 

number

Text

48 Save grid of depth of minimum factor of safety? Enter T (.true.) or F 

(.false.)

49 T

50 Save grid of pore pressure at depth of minimum factor of safety? Enter 

Enter T (.true.) or F (.false.)

51 T

52 Save grid files of actual infiltration rate? Enter T (.true.) or F 

(.false.)

53 T

54 Save grid files of unsaturated zone basal flux? Enter T (.true.) or F 

(.false.)

55 F

56 Save listing of pressure head and factor of safety ("flag")? (Enter -2 

detailed, -1 normal, 0 none)

57 -2

58 Number of times to save output grids

59 1

60 Times of output grids 

61 216000.

62 Skip other timesteps? Enter T (.true.) or F (.false.)

63 F

64 Use analytic solution for fillable porosity? Enter T (.true.) or F 

(.false.)

65 T

66 Estimate positive pressure head in rising water table zone (i.e. in 

lower part of unsat zone)? Enter T (.true.) or F (.false.)

67 T

68 Use psi0=-1/alpha? Enter T (.true.) or F (.false.) (False selects the 

default value, psi0=0)

69 F

70 Log mass balance results? Enter T (.true.) or F (.false.)

71 T

72 Flow direction (enter "gener", "slope", or "hydro")

73 slope

74 Add steady background flux to transient infiltration rate to prevent 

drying beyond the initial conditions during periods of zero 

infiltration? 

75 F