U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods 6–A18

User’s Guide to the Variably Saturated 

Flow (VSF) Process for MODFLOW

Constant-headcell

x

z

1.00

0.90

0.80

0.70

0.60

0.50

0.40

Soil saturation, dimensionless

User’s Guide to the Variably Saturated Flow 

(VSF) Process for MODFLOW

By R. Brad Thoms, Richard L. Johnson, and Richard W. Healy

U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods 6–A18
U.S. Department of the Interior

Dirk Kempthorne, Secretary

U.S. Geological Survey

Mark D. Myers, Director

U.S. Geological Survey, Reston, Virginia: April 2006

Revised and reprinted: November 2006

For sale by U.S. Geological Survey, Information Services

Box 25286, Denver Federal Center

Denver, CO 80225

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

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

Any use of trade, product, or firm names in this publication 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 materials contained within this report.

Suggested citation:

Thoms, R.B., Johnson, R.L., and Healy, R.W., 2006, User’s guide to the Variably Saturated Flow (VSF) Process for MODFLOW: U.S. Geological Survey Techniques and Methods 6–A18, 58 p.

Prepared by the U.S. Geological Survey in Lawrence, Kansas (http://ks.water.usgs.gov)
 iii

Preface

A new Variably Saturated Flow (VSF) Process was written for use with the U.S. Geological Survey 

(USGS) MODFLOW–2000 ground-water model. The VSF Process is designed to simulate threedimensional flow through variably saturated soil media within the vadose zone. The performance 

of this computer program has been tested in models of hypothetical ground-water flow systems; 

however, future applications of the programs may reveal errors that were not detected in the test 

simulations. Users are requested to notify the authors if errors are found in the documentation 

report or in the computer program. Correspondence regarding the report or program should be sent 

to: 

Center for Groundwater Research

OGI School of Science & Engineering

Oregon Health & Science University

20000 NW Walker Road

Beaverton, Oregon 97006–8921 USA

Attention: R. Brad Thoms or Richard L. Johnson

 Although the computer program has been written and used by the USGS, no warranty, expressed 

or implied, is made by the USGS or the United States Government as to the accuracy and functionality of the program and related program material. Nor shall the fact of distribution constitute any 

such warranty, and no responsibility is assumed by the USGS in connection therewith. 

MODFLOW–2000, the VSF Process, and other ground-water programs are available from the 

USGS at the following World Wide Web (WWW) addresses:

http://water.usgs.gov/software/

http://water.usgs.gov/nrp/gwsoftware/

or by anonymous ftp file transfer from directory /pub/software/ground_water/modflow at ftp 

address: water.usgs.gov
iv 
 v

Contents

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

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

Authors' Note . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

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

Overview, Compatibility, and Applicability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

Richards' Equation Flow Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Governing Equation for Variably Saturated Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Soil Characteristic Functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Seepage Face Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Surface Ponding Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Surface Evaporation Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Root Zone Evapotranspiration Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Adaptive Time-Stepping Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Conceptualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Program Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Input Instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Output Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Simulation Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Case 1—One-Dimensional Analytical Solution for Unsaturated Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Case 2—One-Dimensional Field Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Case 3—Two-Dimensional Variably Saturated Transient Infiltration Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Case 4—Two-Dimensional Variably Saturated Transient Drainage Study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
vi 

Case 5—Two-Dimensional Variably Saturated Model Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Case 6—Three-Dimensional Variably Saturated Hypothetical Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Summary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

References Cited. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Appendix 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Modified Picard iteration technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Appendix 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

REF1 Package example input and output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

SPF1 Package example input and output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

PND1 Package example input and output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

SEV1 Package example input and output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

RZE1 Package example input and output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

ATS1 Package example input and output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Appendix 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Custom soil characteristic functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figures 

1. Flowchart representing the structure of the Variably Saturated Flow (VSF) Process and its 

interactionswith MODFLOW–2000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3

2. Graphs showing soil characteristic functions for hydraulicconductivity, K(ψ ), and soil saturation, 

Θ(ψ ),as a function of pore-water pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .8

3. Flowchart of subroutines used by the Richards' Equation Flow Package . . . . . . . . . . . . . . . . . . . . . . . . . . . .9

4. Flowchart of subroutines used by the Seepage Face Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5. Graph showing surface ponding of a medium-coarse soil as a function of relative rainfall rate for two 

initial pressure conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

6–10. Flowcharts of:

6. Subroutines used by the Surface Ponding Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

7. Subroutines used by the Surface Evaporation Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

8. Subroutines used by the Root Zone Evapotranspiration Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

9. Algorithm for the Adaptive Time-Stepping Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

10. Procedures used by the Adaptive Time-Stepping Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

11–14. Graphs showing:

11. Comparison between simulated water saturation (Θ) profiles and Phillip's analytical solution for

one-dimensional, constant-head infiltration in Yolo light clay and sand . . . . . . . . . . . . . . . . . . . . . . 29

12. Model results from the one-dimensional simulation plotted in relation to experimental data and 

the results of a comparable model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

13. Modeled water-table positions for two-dimensional infiltration study plotted in relation to the 

experimental data of Vauclin and others (1979). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

14. Modeled water-table positions for two-dimensional drainage study plotted in relation to experimental data of Vauclin and others (1975). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

15. Diagram of conceptual model for a catchment flow path that feeds a steady irrigation 

canal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

16–19. Graphs showing:

16. Daily precipitation and potential evapotranspiration totals used forcatchment flow-path 

simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
 vii

17. Comparison of VSF and VS2DT model results for depth to water at locations 180 and 

355meters upslope from irrigation canal plotted with daily precipitation. . . . . . . . . . . . . . . . . . . . . 35

18. Comparison of VSF and VS2DT model results for soil saturation at locations 180 and 

355meters upslope from irrigation canal plotted with daily precipitation. . . . . . . . . . . . . . . . . . . . . 36

19. Comparison of VSF and VS2DT model results forcumulative surface soil evaporation. . . . . . . 36

20. Diagram showing conceptual model for hypothetical agricultural drainage that feeds a 

steady irrigation canal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

21. Graph showing simulated depth to water under drainage swale for locations 180 and 

355meters upslopefrom irrigation canal plotted with daily precipitation. . . . . . . . . . . . . . . . . . . . . . . . . . . 38

22. Graph showing simulated ponded-surface cell totals plotted with daily actual 

evapotranspiration computed by RZE1 Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

23. Longitudinal cross-section view showing simulated soil saturations and flow vectors along grid 

row1at start of day 299. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

24. Plan view of simulated soil saturations and ponded cells for day 299 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

25. Example input for Richards' Equation Flow (REF1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

26. Example output for Richards' Equation Flow (REF1) Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

27. Example model budget for Richards' Equation Flow (REF1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

28. Example input forSeepage Face (SPF1) Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

29. Example output forSeepage Face (SPF1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

30. Example input forSurface Ponding (PND1) Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

31. Example output forSurface Ponding (PND1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

32. Example input forSurface Evaporation (SEV1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

33. Example output forSurface Evaporation (SEV1) Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

34. Example model budget forSurface Evaporation (SEV1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

35. Example input for Root Zone Evapotranspiration (RZE1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

36. Example output for Root Zone Evapotranspiration (RZE1) Package. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

37. Example model budget for Root Zone Evapotranspiration (RZE1) Package . . . . . . . . . . . . . . . . . . . . . . . . 55

38. Example input for Adaptive Time-Stepping (ATS1) Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

39. Code listing of FORTRAN FUNCTIONS defining soil characteristic relations . . . . . . . . . . . . . . . . . . . . . . . 57

Tables

1. Adaptive time-stepping results for day 333 of simulation case 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2. Model parameters forconstant-head infiltration in a soil column. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3. Model parameters for Warrick and others (1971) field study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4. Model parameters for Vauclin and others (1979) infiltration study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5. Model parameters for Vauclin and others (1975) drainage study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6. Model parameters for two-dimensional hypothetical comparison case. . . . . . . . . . . . . . . . . . . . . . . . . . . 34

7. Model parameters for three-dimensional hypothetical case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
viii 

Conversion Factors and Abbreviations 

Multiply By To obtain

centimeter (cm) 0.3937 inch (in.)

centimeter per day (cm/d) 0.3937 inch per day (in/d)

centimeter per hour (cm/h) 0.3937 inch per hour (in/h)

cubic centimeter (cm3) 0.06102 cubic inch (in3) 

cubic meter (m3) 35.31 cubic foot (ft3) 

cubic meter per day (m3/d) 35.31 cubic foot per day (ft3/d)

meter (m) 3.281 foot (ft)

meter per day (m/d) 3.281 foot per day (ft/d)

meter per hour (m/h) 3.281 foot per hour (ft/h)

square kilometer (km2) 0.3861 square mile (mi2)

square meter (m2) 0.0929 square foot (ft2)
 ix

Symbols and Dimensions

Number in parentheses refers to the equation where the symbol first appears. Symbols also are 

defined in the text.

Symbol Dimensions Description

C L-1 Specific moisture capacity (2).

C(ψ) L-1 Specific moisture capacity as a function of pore-water pressure head (3).

C(Θe) L-1 Specific moisture capacity as a function of effective soil saturation (7).

CR,CC,CV L2T-1 Hydraulic conductances between node (i,j,k) and an adjacent node (17).

DELRj L Cell width of column j in all rows (17).

DELCi  L Cell width of row i (17).

EV LT-1 Evaporation rate across surface boundary cell (10).

Hn,m L Approximate value for hydraulic head, h, at nth discrete time level and mth

iteration (15).

HCOFi,j,k L2T-1 Sum of coefficients for all source/sink terms and the negative part of the storage term that includes the head in the current time step, n for cell (i,j,k)

(18).

Ki(ψ) LT-1 Hydraulic conductivity as a function of pore-water pressure head in the i

direction (1).

Kc LT-1 Saturated hydraulic conductivity of surface crust (11).

Ksat LT-1 Saturated hydraulic conductivity of the soil (4).

Ktop LT-1 Saturated hydraulic conductivity of top boundary cell (11).

Pi,j,k L2T-1 Sum of head coefficients from source and sink terms for cell (i,j,k) (17).

QET L3T-1 Evapotranspiration flux (12).

Qi,j,k L3T-1 Sum of constants from source and sink terms for cell (i,j,k) (17).

RHSi,j,k L3T-1 Sum of negative constants from source/sink terms and the part of storage 

term that is multiplied by the head from time step m-1 for cell (i,j,k) (18). 

Ss L-1 Specific storage (3).

SSi,j,k L-1 Specific storage for cell (i,j,k) (17).

SRES L-1 Surface resistance factor (10).

THCK L Thickness of surface boundary cell (11).

THICKi,j,k L Vertical thickness of cell (i,j,k) (17).

VOLi,j,k L3 Volume of cell (i,j,k) (12).

W T-1 Volumetric flux per unit volume representing sources and/or sinks (3).

di,j,k L Depth of cell (i,j,k) (12).

droot L Depth of root zone at horizontal location (i,j) (12).

g L/T2 Gravitational acceleration (1).

h L Total hydraulic head, equal to the sum of the pore-water pressure head and 

the elevation head (1).

i,j,k dimensionless Coordinate indices in model grid corresponding to row, column, and layer, 

respectively (1).

ki L2 Intrinsic permeability of the soil medium (1).

kr dimensionless Relative permeability of a soil medium (1).

kr(ψ) dimensionless Relative permeability of a soil medium as a function of pore-water pressure 

head (1).

kr (Θe) dimensionless Relative permeability of a soil medium as a function of effective soil saturation (8).

n dimensionless Time step level (14).
x 

ne dimensionless Porosity of the medium, equal to void volume divided by total bulk volume 

(1).

nv dimensionless Soil parameter (van Genuchten, 1980) representing the degree of pore-size 

uniformity (as nv increases, uniformity increases) (6).

m dimensionless Iteration level (14).

r(z,t) L-2 Root activity as a function of depth and time (12).

t T Time (2).

t

n T Simulation time at time level n (14).

v LT-1 Average pore velocity of the fluid (1).

xi L Spatial coordinate (1).

z L Elevation head (1).

z(θ,t) L Depth as a function of soil-water content and time, t (13).

Θ dimensionless Soil saturation, equal to water volume divided by drainable void volume (2).

Θ(ψ) dimensionless Soil saturation as a function of pore-water pressure (3).

Θe(ψ) dimensionless Effective soil saturation as a function of pore-water pressure (6).

Θr dimensionless Residual soil saturation (6).

Θ0 dimensionless Initial soil saturation (13).

α L-1 Soil parameter (van Genuchten, 1980) representing the inverse characteristic 

length of the soil pores (6).

αc LT2M-1 Matrix compressibility (2).

ßc LT2M-1 Fluid compressibility (2).

θ dimensionless Volumetric soil water content, equal to water volume divided by total bulk 

volume (2).

μ ML-1T-1 Dynamic viscosity of water (1).

ρ ML-3 Density of the active fluid (that is, water) (1).

φn dimensionless nth order function of soil moisture, θ (13).

ψ L Pore-water pressure head (1).

ψA, ψB L Interval between which pore-water pressure values are stored in internal 

table (9).

ψa L Atmospheric pressure potential (10).

ψi , ψi+1 L Specified pore-water pressure values logarithmically distributed between 

interval ψA, ψB within internal table (9).

ψi,j,k L Pore-water pressure head within cell (i,j,k) (12).

ψroot L Pressure head equal to the permanent wilting point of the plant roots (12).

ψs L Pore-water pressure head of surface boundary cell (10).

Symbol Dimensions Description
Abstract 1

User’s Guide to the Variably Saturated Flow (VSF) Process for

MODFLOW

By R. Brad Thoms1, Richard L. Johnson1, and Richard W. Healy2

Abstract

A new process for simulating three-dimensional (3-D) variably saturated flow (VSF) using Richards’ equation has been added 

to the 3-D modular finite-difference ground-water model MODFLOW. Five new packages are presented here as part of the VSF 

Process—the Richards' Equation Flow (REF1) Package, the Seepage Face (SPF1) Package, the Surface Ponding (PND1) Package, 

the Surface Evaporation (SEV1) Package, and the Root Zone Evapotranspiration (RZE1) Package. Additionally, a new Adaptive 

Time-Stepping (ATS1) Package is presented for use by both the Ground-Water Flow (GWF) Process and VSF. The VSF Process 

allows simulation of flow in unsaturated media above the ground-water zone and facilitates modeling of ground-water/surfacewater interactions. 

Model performance is evaluated by comparison to an analytical solution for one-dimensional (1-D) constant-head infiltration 

(Dirichlet boundary condition), field experimental data for a 1-D constant-head infiltration, laboratory experimental data for twodimensional (2-D) constant-flux infiltration (Neumann boundary condition), laboratory experimental data for 2-D transient drainage through a seepage face, and numerical model results (VS2DT) of a 2-D flow-path simulation using realistic surface boundary 

conditions. A hypothetical 3-D example case also is presented to demonstrate the new capability using periodic boundary conditions (for example, daily precipitation) and varied surface topography over a larger spatial scale (0.133 square kilometer). The 

new model capabilities retain the modular structure of the MODFLOW code and preserve MODFLOW’s existing capabilities as 

well as compatibility with commercial pre-/post-processors. The overall success of the VSF Process in simulating mixed boundary 

conditions and variable soil types demonstrates its utility for future hydrologic investigations.

This report presents a new flow package implementing the governing equations for variably saturated ground-water flow, four 

new boundary condition packages unique to unsaturated flow, the Adaptive Time-Stepping Package for use with both the GWF 

Process and the new VSF Process, detailed descriptions of the input and output files for each package, and six simulation examples 

verifying model performance.

Introduction

The increasing focus on flow and transport interactions between ground and surface water necessitates a more sophisticated 

method for simulating variably saturated flow using the three-dimensional (3-D), modular finite-difference ground-water flow 

model MODFLOW (Harbaugh and others, 2000). Specifically, there is an increasing demand for a MODFLOW-compatible tool 

that allows practitioners to simulate vadose zone flow within a small-scale area of concern (for example, an agricultural field) in 

a way that facilitates interaction with existing regional ground-water MODFLOW models. The Variably Saturated Flow (VSF) 

Process is introduced here as an approach separate from the existing MODFLOW Ground-Water Flow (GWF) Process (Harbaugh 

and others, 2000). In a new flow package within the VSF Process, the governing equation for saturated ground-water flow is modified into the “mixed form” of Richards’ equation (Huyakorn and Pinder, 1983) for variably saturated ground-water flow. In addition, four new boundary conditions unique to unsaturated flow are presented here as packages within VSF. The new Adaptive 

Time-Stepping Package also is introduced to increase MODFLOW's efficiency during lengthy transient simulations and is implemented for use with both the existing GWF Process and the new VSF Process. 

The VSF Process expands MODFLOW's spatial domain to include the vadose zone. This requires a more general governing 

equation and several new boundary conditions. The Richards’ Equation Flow (REF1) Package uses soil characteristic functions to 

update the unsaturated soil’s drainage and storage properties as the soil pore-water pressure changes during a simulation. Interactions at the ground surface are represented by four new boundary condition packages—evaporation from the surface soil layer 

1Oregon Health & Science University, Portland, Oregon.

2U.S. Geological Survey, Lakewood, Colorado.
2 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

(SEV1), evapotranspiration from the root zone (RZE1), seepage faces (SPF1), and surface ponding (PND1). Expansion of the 

MODFLOW domain to include the vadose zone gives the model the physical basis for rainfall partitioning at the surface boundary 

by representing infiltration through the vadose zone to the aquifer from precipitation data (for example, Downer and Ogden, 2004). 

This enables the user to specify precipitation rates directly in the existing Recharge (RCH6) Package. Constant-head boundaries 

also can be specified using the existing MODFLOW Basic (BAS6) Package. 

Due to the nonlinear nature of Richards’ equation, the VSF Process requires input specifications that are different than typical 

saturated ground-water simulations. It requires finer spatial discretization within the vadose zone and temporal scales that are on 

the order of hours or days. Model cells above the water table are prevented from drying by REF1, and the specific yield for those 

cells is replaced by the specific moisture capacity in the VSF governing equation. The conductances for these cells are computed 

by multiplying the area of the cell face with the saturated hydraulic conductivity and the relative permeability (a function of the 

pore-water pressure with values between 0 and 1). Both the relative permeability and the specific moisture capacity are updated 

each iteration from the latest total head value for that cell. Additionally, to alleviate convergence problems in steady-state simulations, an option is included that allows the user to specify an initial condition of hydrostatic equilibrium. 

VSF is implemented in a modular fashion that allows the use of existing GWF Packages within model cells that remain fully 

saturated during a given simulation. Additionally, existing head-dependant boundary conditions (for example, river, stream, lake, 

drain) specified within cells that become unsaturated also will function without code modifications as long as the resultant leakage 

rate (unit downward hydraulic gradient) remains less than the cell's saturated hydraulic conductivity. Convergence problems will 

occur for large specified source/sink fluxes (for example, using the well package) within unsaturated cells. The VSF Process is 

currently (2005) incompatible with the Ground-Water Transport (GWT), Observation (OBS), Sensitivity (SEN), and ParameterEstimation (PES) Processes as well as current transport models such as MT3DMS (Zheng and Wang, 1999) or RT3D (Clement, 

1997). However, the Richards’ Equation Flow (REF1) Package is designed to output the water-filled porosity distribution that 

would be needed for simulation of contaminant transport in the unsaturated zone, facilitating future transport model development 

(Thoms, 2003). Additionally, the modular design of the VSF Process will facilitate updates as future versions of MODFLOW 

are released.

The purpose of this report is to present a new process for simulating 3-D variably saturated flow (VSF) using Richards’ equation. The new VSF Process can be added to the 3-D modular finite-difference ground-water model MODFLOW. The report also 

includes four new boundary condition packages, the Adaptive Time-Stepping Package for use with both the GWF Process and the 

new VSF Process, detailed descriptions of the input and output files for each package, and six simulation examples verifying model 

performance.

Authors’ Note

The VSF approach is intended primarily for applications where unsaturated zone thicknesses are small (less than 50 m), seasonal water-table fluctuations are significant, and vadose zone flow is highly variable. VSF is an effective tool for examining nearsurface flow behavior such as subsurface stormflow, perched water tables, root zone evapotranspiration losses, and variable source 

areas. However, the high computational demands of the VSF Process preclude its reasonable application to the regional scale (for 

example, greater than 100 km2

). 

Readers should be aware of another MODFLOW package. The Unsaturated Zone Flow (UZF1) Package (Niswonger and 

others, 2006) is designed for coupling watershed-runoff models to MODFLOW–2000, but it can be used without that coupling. It 

uses a kinematic wave approach to simulate one-dimensional (1-D) vertical flow through the unsaturated zone. UZF1 provides an 

efficient means of simulating recharge in MODFLOW–2000 that accounts for evapotranspiration and storage in the unsaturated 

zone. UZF1 offers a less rigorous, and less computationally demanding, treatment of flow through the unsaturated zone than that

provided by VSF.

Acknowledgments

This work was partially supported by the National Water-Quality Assessment Program of the U.S. Geological Survey. Additional support was provided by the Department of Environmental & Biomolecular Systems at the OGI School of Science & Engineering at Oregon Health & Science University (Portland, Oregon). 
Overview, Compatibility, and Applicability 3

Overview, Compatibility, and Applicability

Vadose zone flow traditionally is incorporated within MODFLOW simulations through recharge estimation methods and the 

use of specific yield for unconfined cells. The approach presented here allows MODFLOW to simulate 3-D, variably saturated 

flow in a physically rigorous way using Richards’ equation. An overview of the VSF Process operation within MODFLOW is 

presented in figure 1. VSF is structured in the same fashion as GWF; packages perform specific hydrologic capabilities (for example, PND1 simulates surface ponding, SEV1 simulates surface evaporation), and each package contains procedures that perform 

similar functions (for example, allocate reserves array space).

The ability to simulate fully saturated cells is retained using modified versions of the Block Centered Flow (BCF6) and Layer 

Property Flow (LPF1) Packages. The new versions of these packages operate within VSF to formulate each cell’s conductance 

and storage properties for partially saturated cells while their original functionality remains unchanged for fully saturated cells. 

For unsaturated cells, LPF1 and BCF6 neglect the specific storage and specific yield terms and set the saturated thickness equal to 

the cell thickness. The new Richards’ Equation Flow (REF1) Package scales the unsaturated cell’s saturated conductance and leakance values by the relative permeability, scales the specific storage by soil saturation, and adds the specific moisture capacity to 

the storage term. The user must specify model layers containing the anticipated vadose zone in the REF1 input file, and those 

layers must be specified convertible or unconfined in the LPF1 or BCF6 Packages, respectively. 

PND? Y

GLOBAL: Define

VSF: Allocate

REF1, PND1, SPF1, SEV1, RZE1

VSF: Formulate

LPF1* / BCF6*

VSF: Formulate

REF1, SEV1, RZE1, SPF1

VSF: Budget

REF1, PND1, SPF1, SEV1, RZE1

VSF: Budget

LPF1* / BCF6*

MODFLOW–2000 Main

VSF:Read and Prepare

REF1, SEV1, RZE1

PND1 / SPF1 – Reset flags, IBOUND = 1

Stress-period

loop

Time-step

loop

Iteration

loop

VSF: Restart

Surface ponding? IBOUND = -1, 

Restart time step.

Y

New Packages

REF1: Richards’ Equation Flow

PND1: Surface Ponding

SPF1: Seepage Face

SEV1: Soil Evaporation

RZE1: Root Zone Evapotranspiration

LPF1*: Modified Layer Property Flow

BCF6*: Modified Block Centered Flow

Bold Text denotes Processes

Italic Text denotes Procedures

* Modified for use with VSF

Figure 1. Structure of the Variably Saturated Flow (VSF) Process and its interactions with MODFLOW–2000. 

Simulations are controlled by MODFLOW's main program (left) and proceed from top to bottom, cycling through 

the nested iteration, time-step, and stress-period loops.
4 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

The Observation, Sensitivity, and Parameter Estimation Processes (Hill and others, 2000) as well as the Hydrologic Unit Flow 

Package (Anderman and Hill, 2000) are not supported by the VSF Process. VSF is structured so that, with some user discretion, 

existing GWF packages can be used within cells that remain fully saturated during a given simulation. The increased spatial discretization required by VSF will cause problems with surface-interaction packages such as the River (RIV6), Lake (LAK3), and 

Stream (STR6) Packages that specify leakance terms along a linear sequence of cells. For example, surface-water bodies that are 

wider than the (reduced) cell width required to maintain a reasonable width-to-depth aspect ratio may not be effectively modeled 

with the VSF capability active. Additionally, the seepage face and surface ponding boundary conditions use the IBOUND array 

to make adjustments to the surface cells on the basis of simulation conditions. These adjustments will override all other active 

boundary conditions in shared cells specified by other packages (for example, the Recharge, General Head Boundary, or TimeVariant Specified Head Packages). Careful definition of boundary conditions is required for complex simulations with multiple 

boundary conditions to avoid unwanted interference. However, appropriate use of the Seepage Face (SPF1) and Surface Ponding 

(PND1) Packages will enable MODFLOW to directly specify precipitation data as the ground-surface flux using the existing 

Recharge (RCH6) Package. Constant-head boundaries also can be specified using the existing Basic (BAS6) Package. 

Richard’s equation typically is used to simulate vadose zone fluxes on spatial scales ranging from tens of meters (for example,

Scanlon and others, 2002) to hundreds of meters (for example, Downer and Ogden, 2004). In conjunction with existing 

MODFLOW tools such as local grid refinement (Mehl and Hill, 2002), VSF can be a powerful tool for the user to examine smaller 

scale surface-vadose zone interactions within the context of a larger regional ground-water model. Quality data sets containing 

high collection frequency (hourly or daily data) will allow the modeler to examine in detail transient processes such as lateral subsurface stormflow that GWF is typically not able to represent. However, there are two major limiting factors in applying Richards’ 

equation to large-scale systems: (1) realistic parameter estimation of capillary saturation characteristics for the domain of interest, 

and (2) computational efficiency of large 3-D simulations. Depending on the specific goals of a particular study, the first limitation 

may be overcome by simplifying assumptions or strategic field data collection. The second limitation is overcome only through 

increased technological capability (that is, a fast computer). 

There are a few constraints limiting VSF effectiveness under certain conditions. Large domains (greater than 100,000 nodes) 

will require lengthy simulation times and a large file-storage capacity. The nonlinearity introduced by the soil characteristic functions to the hydraulic-conductivity and storage terms requires fine spatial discretization (Δz less than 30 cm) to produce accurate 

results, especially near boundaries with strong hydraulic gradients. Extremely nonlinear soils (for example, coarse sand) and/or 

boundary conditions (for example, large rainfall on initially dry soil) as well as complex heterogeneous systems will cause convergence problems. These difficulties often can be overcome by trying different solvers or adjusting the solver’s damping parameters and iteration settings. General troubleshooting tips include checking the model’s mass balance, adjusting the adaptive timestepping parameters (for example, decrease the minimum time step), refining the spatial grid, avoiding assigned fluxes greater than 

the soil’s saturated hydraulic conductivity, increasing the solver’s maximum iterations, and setting a realistic initial head profile 

(for example, using the solution from a previous simulation run). Simulation efficiency can be maximized by increasing the spatial 

grid discretization gradually and comparing simulation results and mass balance as the grid size increases. Additionally, nonlinear 

initial conditions (for example, large head gradients) specified for steady-state simulations will cause convergence problems. For 

most steady-state applications, the user is advised to employ the hydrostatic equilibrium option by setting the IHYDST flag in the 

REF1 Package equal to 1, which specifies the total head in each unsaturated cell equal to the water-table location defined by the 

BAS6 Package.

Richards’ Equation Flow (REF1) Package

Conceptualization

The VSF Process includes several packages that convert MODFLOW into a more general 3-D variably saturated flow model. 

The Richards’ Equation Flow (REF1) Package converts MODFLOW’s governing equation from a special (saturated) case of Richards’ equation to the full equation for variably saturated flow. REF1 works concurrently with a modified version of the active flow 

package (BCF6 or LPF1) to incorporate pressure-head dependence within the hydraulic-conductivity and storage terms for unsaturated cells (that is, total head below nodal elevation). For these cells, the saturated thickness within BCF6 or LPF1 is set equal 

to the cell thickness, and the saturated storage terms are ignored. Within REF1, the saturated conductance and leakance values are 

scaled by the relative permeability, kr, from the current head in that cell, and the unsaturated storage terms are computed. The 

unsaturated soil’s storage and relative permeability properties are defined by soil characteristic functions that describe the unsaturated behavior as the pore-water pressure changes. Under saturated conditions, these characteristic functions converge to 1, 

allowing the modified governing equation to converge back to the original MODFLOW formulation for saturated cells. The formulation of the new governing equation and examples of the soil characteristic functions are described in this section.
Richards’ Equation Flow (REF1) Package 5

Governing Equation for Variably Saturated Flow 

Three-dimensional, variably saturated isothermal fluid flow in a heterogeneous porous medium is described by an equation 

that is derived by incorporation of Darcy’s law into the conservation of mass equation for fluid flow (Lappala and others, 1987). 

Darcy’s law for variably saturated media can be expressed as: 

, (1)

where

i = coordinate indices (i = 1, 2, 3);

v = average pore velocity of the fluid [L/T];

kr(ψ) = relative permeability of a particular medium (soil type) as a function of pressure head, ψ [dimensionless];

ψ = pore-water pressure head (ψ = h-z) [L];

z = elevation head [L];

ki = intrinsic permeability of the medium [L2];

ρ = density of the active fluid (that is, water) [M/L3];

g = gravitational acceleration [L/T2] ;

h = total hydraulic head (equal to the sum of the average pore-water pressure head, ψ, and the elevation head, z [L];

ne = porosity of the medium (void volume/total bulk volume) [dimensionless];

μ = dynamic viscosity of water [M/LT];

xi = spatial coordinate [L]; and

Ki(ψ) = hydraulic conductivity as a function of pore-water pressure head [L/T].

The continuity equation represents the conservation of mass during fluid flow through an elemental volume of the porous 

medium. This equation takes the form (Freeze, 1978): 

, (2)

where

Θ = saturation of the soil (water volume/void volume) [dimensionless]; and

t = time [T].

The bracketed terms on the right-hand side of equation 2 describe the storage properties of the medium. These account for 

changes in liquid stored in the elementary volume due to changes in soil saturation (Θ), changes in the soil pore space (ne), and 

compression or expansion of the liquid (changes in the fluid density, ρ). For a given soil type or geologic formation, the porosity, 

density, and soil-moisture content are assumed to be solely dependent on ψ. Several terms can be defined to simplify the righthand side of equation 2 (Lappala and others, 1987):

Specific moisture capacity [L-1]: ;

Matrix compressibility [LT2/M]: ;

Fluid compressibility [LT2/M]: ; and

Specific storage [L-1]: ;

where

θ = volumetric soil water content (Θ . ne) (water volume / total bulk volume) [dimensionless].

vi kr( ) ψ

ki

ne

----ρg

μ

------ ∂h

∂xi

------ Ki

( ) ψ

ne

-------------- ∂h

∂xi

= = ------

∂ neρvx ( )

∂x --------------------- ∂ neρvy ( )∂y --------------------- ∂ neρvz ( )∂z + + --------------------- ∂h

∂t

----- ρne

∂Θ

∂h

------- ρΘ

∂ne

∂h -------- neΘ∂ρ∂h = + + ------

C ∂θ

∂ψ = -------

αc

∂ne

∂ψ = -------

βc

1

ρ

-- ∂ρ

∂ψ = -------

Ss ρg neβc + αc = ( )
6 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

For simplification, it frequently is assumed that the fluid and soil matrix are only slightly compressible and the specific storage 

term, Ss, is constant for a given soil type or geologic formation. Ss represents the volume of water per unit aquifer area that a unit 

volume of soil matrix releases (or takes up) in response to a unit decrease (or increase) in hydraulic head (Dingman, 2002). This 

assumption, combined with equations 1 and 2, results in what is known as Richards’ equation (Huyakorn and Pinder, 1983):

, (3)

where W = volumetric flux per unit volume representing sources and (or) sinks [T-1].

This parabolic partial differential equation is nonlinear due to the nature of the hydraulic conductivity K(ψ) and specific moisture capacity C(ψ) functions. To alleviate convergence difficulties arising from steep gradients in soil moisture storage, the specific moisture capacity term is linearized using the modified Picard iteration technique (Celia and others, 1990). A detailed discussion on this technique is included in Appendix 1. The other nonlinear function, K(ψ), is defined here as the product of the 

relative permeability, kr(ψ), and the saturated hydraulic conductivity (Ksat):

, (4)

where kr(ψ) has values between 0 and 1, depending upon the pore-water pressure, ψ, of the soil. 

As ψ becomes positive and the pores saturate with water, the specific moisture capacity C(ψ) converges to 0, the relative permeability kr(ψ) converges to 1, and the soil saturation Θ converges to 1. Thus, equation 3 converges to the general flow equation 

used by MODFLOW (equation 1, Harbaugh and others, 2000) for saturated subsurface flow:

. (5)

The new formulation (equation 3) allows MODFLOW to simulate a subsurface domain from the ground surface through the 

unsaturated (vadose) zone past the water table into the saturated aquifer. In effect, VSF is a more sophisticated replacement for 

unconfined cells used in MODFLOW. However, as discussed below, it should be noted that VSF assigns cell properties on the 

basis of total head at nodal elevations, whereas the GWF Process scales the cell conductance and storage properties by the cell’s 

saturated thickness. For GWF this is appropriate because cell thicknesses at the water table often are quite large, whereas for VSF 

cell thicknesses should be small (less than or equal to 30 cm).

Soil Characteristic Functions

For MODFLOW to simulate variably saturated porous media flow, several relationships describing the unsaturated flow and 

storage properties intrinsic to a modeled soil must be defined. The soil characteristic functions K(ψ), Θ(ψ), and C(ψ) have been 

represented in the literature by several different empirical and theoretical methods (Millington and Quirk, 1961; Brooks and Corey, 

1966; Mualem, 1976; Haverkamp and others, 1977; van Genuchten, 1980; Broadbridge and White, 1988). The REF1 Package 

uses the popular nonhysteretic (singular) van Genuchten-Mualem soil characteristic functions. These functions take the form: 

Effective saturation as a function of pore-water pressure head (van Genuchten, 1980):

; (6)

∂x

∂ Kx( ) ψ

∂h

∂x

----- ∂

∂y

----- Ky( ) ψ

∂h

∂y

----- ∂

∂z

----- Kz( ) ψ

∂h

∂z +++ ----- W ∂h∂t

----- Θ( ) ψ S = [ ] s + C( ) ψ

K( ) ψ k = r( ) ψ Ksat

∂

∂x

----- Ksat

x ∂h

∂x

----- ∂

∂y

----- Ksat

y ∂h

∂y

----- ∂

∂z

----- Ksat

z ∂h

∂z +++ ----- W Ss

∂h

∂t = -----

Θe( ) ψ

Θ Θr –

1 Θr – ----------------

1

1 ( ) αψ nv [ ] + m --------------------------------------- ψ < 0

⎩ 1 ψ ≥ 0

⎪

⎨

⎪

⎧

= =
Richards’ Equation Flow (REF1) Package 7

Specific moisture capacity as a function of effective saturation:

; (7)

Relative permeability as a function of effective saturation (Mualem, 1976):

; (8)

where

Θr = residual soil saturation [dimensionless];

nv, m = soil parameters (m = 1-1/nv) representing the degree of pore-size uniformity (as nv increases, uniformity increases) 

[dimensionless]; and

α = parameter representing the inverse characteristic length of the soil pores [L-1].

These functions are used to evaluate the coefficients in equation 3 iteratively as the model solves for hydraulic head. Examples of 

unsaturated hydraulic conductivity and soil saturation behavior are shown in figures 2A and 2B, respectively. The figures show 

that pore-water pressure (suction) has a strong, nonlinear effect on both hydraulic conductivity and saturation of the soil.

The soil saturation and specific moisture capacity functions are used to iteratively evaluate the unsaturated storage properties 

for the current head in each cell (equations 6 and 7). It is noted that the soil pressure is evaluated at the nodal elevation (that is, 

cell center). This is the same approach used, for example, by VS2D (Lappala and others, 1987). Thus within VSF, the cell is

treated as fully saturated when the total head is above the cell center elevation and unsaturated when the head is below it. In a 

physical context, the capillary rise for most soil media is 10 cm or more (van Genuchten, 1978), which is on the order of the halfthickness of a typical VSF cell. Therefore, if the simulated water-table elevation rises to the model cell’s nodal elevation, it can 

be assumed that capillary rise will cause the cell to become saturated. Intercell averages of the relative permeability terms are 

computed for the current head values and multiplied by the saturated transmissivity and leakance arrays, giving the unsaturated

hydraulic-conductivity values in equation 3. There are numerous methods suggested in the literature for computing the intercell 

average of relative permeability values (Brutsaert, 1971; Warrick and others, 1991; Zaidel and Russo, 1992; Baker, 1995; Desbarats, 1995) for block-centered, finite-difference models. In general, the appropriate method depends on the hydraulic properties of 

the system and how well the system is characterized. In VSF, the user can specify the averaging method in the input file as arithmetic mean, geometric mean, or upstream weighting (the relative permeability of the cell with the greater total head is used in the 

intercell computation).

At the start of each simulation an internal table of water contents, hydraulic conductivities, and specific moisture capacities

is generated for each soil type to minimize the computationally intensive direct evaluation of the soil hydraulic functions (Simunek 

and others, 1995). The user specifies the interval of pore-water pressure (suction) head expected within the simulation (ψA, ψB), 

and the values of Θ, kr, and C are evaluated at prescribed pore-water pressure heads (ψi

) within the interval. These prescribed 

pressure heads are distributed logarithmically within the interval such that: 

constant , (9)

where

ψA > ψi+1 > ψB,

ψA > ψi > ψB.

During the simulation, the values for Θ, kr , and C within a model cell are computed for a given head value by linearly interpolating 

between the corresponding values in the table. If the given head value is outside of the table’s specified interval (ψA, ψB), the 

hydraulic characteristics are evaluated directly from the soil characteristic functions.

Note to advanced MODFLOW users: The majority of VSF simulations can be performed using the default functions described 

above and an interpolation utility such as RETC (van Genuchten and others, 1997) to determine the best-fit van Genuchten-Mualem 

soil characteristics for the soil(s). However, several options for customizing the soil characteristic functions are included to 

(1) allow users to define a separate set of parameters (nv , m) for the relative permeability, kr(ψ), and soil saturation, Θ(ψ), relationships, respectively (the same parameter set is used to define both the soil saturation and the specific moisture capacity relationships); (2) allow the model to use soil characteristic functions other than the default van Genuchten-Mualem relationships 

previously defined (equations 6–8) (for example, Vauclin and others, 1975, or Haverkamp and others, 1977), and (3) provide a 

C Θe ( ) ∂θ

∂ψ------- nemα 1 Θr ( ) –

1 – m ----------------------------------Θe

1/m 1 Θe1/m ( ) – m

= = –

kr Θe () Θe

1/2 1 1 Θe1/m ( ) – m [ ] – 2

=

ψi + 1/ψi=
8 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

0.001

0.01

0.10

10

100

Saturated hydraulic conductivity, Ksat

(A)

1.0

0.3

0.4

0.5

0.6

0.7

0.8

0.9

-150 -200 -250 -300

Pore pressure, , in centimeters

0 -50 -100

Residual soil saturation, r

(B) 1.0

Soil saturation, , dimensionless Θ Hydraulic conductivity, 

K, in centimeters per day

Figure 2. Soil characteristic functions for (A) hydraulic conductivity, K(ψ), and (B) soil saturation, Θ(ψ), as 

a function of pore-water pressure. This soil type is used in simulation example case 2.
Richards’ Equation Flow (REF1) Package 9

logical structure within the FORTRAN code for creating user-defined soil characteristic functions. The soil characteristic 

functions consist of separate subroutines that can be customized with a working knowledge of FORTRAN and a compiler. 

Several of the example cases presented in this document (cases 1 and 4) use functions different from the default van GenuchtenMualem relationships. A full explanation of the options available for custom soil characteristic functions and the code listing are 

presented in Appendix 3.

Program Description

The REF1 Package is built in the modular format characteristic of all MODFLOW packages and operates in conjunction with 

either the BCF6 Package or the LPF1 Package. The subroutines within REF1 are constructed to fit the standard MODFLOW process-package-procedure framework (fig. 3); VSF1REF1AL allocates space for the saturation or water-filled porosity output control variables, the saturated leakance, the soil characteristic function parameters, and the associated internal tables; VSF1REF1RP 

reads in the user settings and parameters for unsaturated model cells, generates starting head distribution from initial saturations 

(if requested) using SVSF1REF1PS, and creates the internal tables relating soil pore-water pressure to relative permeability using 

SVSF1REF1KP, to soil saturation using SVSF1REF1SP, and to specific moisture capacity using SVSF1REF1CP; VSF1REF1FM 

calculates the storage of unsaturated cells using soil saturation and specific moisture capacity values from the internal table and 

uses the relative permeability values from the internal table to adjust the conductance and leakance terms computed by 

VSF1LPF1SP / VSF1LPF1FM or VSF1BCF6FM; and VSF1REF1BD outputs soil saturation or water-filled porosity (if 

requested), computes the unsaturated storage, and stores the total storage in the flow budget array, VBVL. Selected subroutines 

of the BCF6 and LPF1 Packages are modified to include an integer flag array of unsaturated cell locations for use in the VSF Process. The VSF1LPF1FM and VSF1BCF6FM subroutines are modified to ignore storage calculations and leakage corrections for 

the user-specified unsaturated cells. The VSF1LPF1SP and SVSF1LPF1N subroutines are modified to pass the unsaturated cell 

locations flag array to the SVSF1LPF1VCOND and SVSF1LPF1HCOND subroutines. The SVSF1LPF1VCOND, 

SVSF1LPF1HCOND, and SVSF1BCF6H subroutines set the saturated thickness of unsaturated cells equal to the cell thickness. 

The SVSF1LPF1B, SVSF1BCF6B, SVSF1LPF1F, and SVSF1BCF6F subroutines compute flow using the entire cell face for 

unsaturated cells. The SVSF1LPF1S and SVSF1BCF6S subroutines ignore the storage calculations for unsaturated cells. The 

REF1 subroutines are written in standard FORTRAN90 and should be compatible with any FORTRAN90 compiler.

MODFLOW–2000 Main

VSF1REF1BD

SVSF1LPF1N*

SVSF1LPF1VCOND*

SVSF1LPF1HCOND* Modified subroutines*

VSF1REF1AL VSF1REF1RP VSF1REF1FM

SVSF1BCF6H*

VSF1LPF1SP*

VSF1LPF1FM*

VSF1BCF6FM*

LPF or BCF? LPF or BCF?

SVSF1LPF1S*

SVSF1LPF1F*

SVSF1LPF1B*

SVSF1BCF6S*

SVSF1BCF6F*

SVSF1BCF6B*

SVSF1REF1KP

SVSF1REF1SP

SVSF1REF1CP

SVSF1REF1PS

SVSF1REF1KP

SVSF1REF1SP

SVSF1REF1CP

SVSF1REF1KP

SVSF1REF1SP

SVSF1REF1CP

Figure 3. Subroutines used by the Richards’ Equation Flow (REF1) Package. A typical simulation proceeds from left to 

right. Existing subroutines that have been modified are contained in a dashed box.
10 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Input Instructions

Input for the REF1 Package is read from a file specified in the MODFLOW name file, identified as the ‘REF’ file type. Free 

format is used for reading all variables. The REF1 Package is used in conjunction with either the LPF1 Package or BCF6 Package

only; the HUF2 Package currently (2005) is not supported by the VSF Process. Packages from the Ground-Water Flow (GWF) 

Process can only be used within fully saturated cells, which must be designated as such in the integer flag array, STYPE (STYPE=

0).

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. IREFCB

Record 2 is only read if 

2. IREFOC ISATOC

3. SLNUM NTAB

4. IAVG ISATFL IHYDST

Record 5 is only read if IREFOC > 0

5. REFOC (IREFOC, 2)

6. PTABA PTABB

Read record 7 for each soil type (SLNUM types)

7. ISC

ISC defines the parameter set for the soil characteristic function of each soil type as either singular or dual.

Read record 8 once for each soil type (SLNUM types) if the parameter set is singular (|ISC|=1) or twice if the parameter set is dual 

(|ISC|>1). 

8. ALPHA VGN

Read record 9 once for each soil type (SLNUM types)

9. RSAT EFFP

Record 10 is only read if ISATFL > 0

10. SAT (NCOL, NROW, NLAY) - U2DREL.

11. STYPE(NCOL, NROW, NLAY) - U2DINT.

Explanation of the variables read by the REF1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1. 

1. IREFCB—The flag containing the file unit for saving cell-by-cell soil saturations.

> 0, requested output will be saved to an ASCII file with the ‘.SAT’ extension located at unit IREFCB.

= 0, no output.

< 0, requested output will be saved to a BINARY file with the ‘.SAT’ extension located at unit IREFCB. 

2. IREFOC—The flag for REF1 cell-by-cell output control.

> 0, output recorded at times specified in the REFOC array.

< 0, output recorded at all times.

ISATOC—The flag indicating output of either soil saturations or water-filled porosities.

> 1, requested output recorded as soil saturations.

< 1, requested output recorded as water-filled porosities.

Read only if .

3. SLNUM—The flag indicating the number of soil types (sets of unsaturated hydraulic parameters) defined in the model 

domain.

NTAB—The total number of specified soil pore-water pressure values within the internal soil characteristics table.

IREFCB 0 ≠

IREFCB 0 ≠
Richards’ Equation Flow (REF1) Package 11

4. IAVG—The flag specifying the intercell averaging of the relative permeability.

= 1, Geometric mean.

= 2, Arithmetic mean.

= 3, Upstream weighted (the relative permeability of the cell with the greater total head is used).

ISATFL—The flag indicating if initial soil saturations are to be read to set the initial head distribution.

< 0, initial heads read in the BAS file are used.

> 0, initial saturations are read as a real array (SAT) for all model cells.

IHYDST—The flag indicating that the initial head distribution for unsaturated cells will be set to hydrostatic equilibrium 

(total head = water-table elevation specified in BAS input file).

< 0, initial heads in unsaturated cells remain unchanged.

> 0, initial heads in unsaturated cells are set equal to the water-table elevation.

5. REFOC (IREFOC, 2) —The two-dimensional (2-D) output control array containing real values corresponding to the stress 

period and time step or total simulation time. There are IREFOC entries corresponding to the total number of requested 

output times. The first value in each entry contains the requested stress period. The second value contains either a negative

integer corresponding to the time step number or a positive real number corresponding to the total simulation time. Read 

only if IREFOC > 0.

Note: If the adaptive time-stepping (ATS1) is active, the requested output must correspond to output times requested by the 

Output Control (OC) Option.

6. PTABA—The lower bounds (near saturation) of the specified pore-water pressure interval for the internal soil 

characteristics table.

PTABB—The upper bounds (high suction) of the specified pore-water pressure interval for the internal soil characteristics 

table.

7. ISC (SLNUM) —The integer flag array that determines whether the soil characteristic function for each soil type is defined 

by a singular or dual parameter set (ALPHA, VGN) and determines the form of the function used.

ISC = 1, singular set of parameters are used in the van Genuchten-Mualem functions.

ISC = 2, dual set of parameters are used in the van Genuchten-Mualem functions.

ISC < 0, singular/dual set of parameters are used in custom functions (see Appendix 3).

8. ALPHA (SLNUM, 2) —The real array containing van Genuchten soil parameter, α.

VGN (SLNUM, 2) —The real array containing van Genuchten soil parameter, nv.

Note: For dual parameter sets (|ISC|>1), the first set of entries defines the relative permeability relationship, and the second 

set of entries defines the soil saturation/specific moisture capacity relationships. Only one set of entries is read if |ISC|=1. 

See Appendix 3 for more information.

9. RSAT (SLNUM)—The real array containing the residual soil saturation, Θr.

EFFP (SLNUM) —The real array containing the soil porosity, ne. 

10. SAT (NCOL, NROW, NLAY) —The real array containing initial soil saturations. 

Only read if ISATFL > 0.

11. STYPE (NCOL, NROW, NLAY) —The integer array containing the soil type corresponding to the set of hydraulic 

parameters (STYPE = [1...SLNUM]). 

Note: Setting STYPE (i,j,k) equal to zero deactivates the REF1 functionality for cell (i,j,k), and the cell is treated as fully 

saturated by the active flow package (LPF1 or BCF6).

Output Data

Simulation results for saturation or water-filled porosity are recorded to a file with a ‘.SAT’ extension for user-requested 

times. The output is recorded by the ULASAV (ASCII format) or ULASV2 (BINARY format) subroutines. The format of these 

files is identical to the head (HED) and drawdown (DRW) output files, allowing existing post-processors to read and process the

results efficiently. Example output can be found in Appendix 2.
12 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Seepage Face (SPF1) Package

Conceptualization

Expansion of the MODFLOW simulation domain into the vadose zone introduces several boundary conditions at the ground 

surface. The Seepage Face (SPF1) Package enables MODFLOW to accurately simulate flux out of the subsurface domain across 

the ground-surface interface via a “seepage face” (for example, areas with steep surface topography that intersect with the 

water table).

The seepage face boundary condition is applied to cells within the model domain defined by the user as an interface where 

ground water could potentially exfiltrate during the simulation period. The user specifies in the input file the locations of all the 

cells that might become saturated during the course of the simulation. The specified cells do not allow flow across the surface 

boundary as long as the water table remains below the respective nodal elevation. Once the water table rises to saturate the cell, 

the cell is converted into a constant-head cell with the head set to the nodal elevation value. This is equivalent to setting the cell’s 

pore-water pressure equal to atmospheric pressure (that is, zero gauge pressure). Once a cell is set as a constant-head cell, the total 

flow through all the cell faces is evaluated. The cell is converted back to a variable-flow cell if the total flow is positive (into the 

model domain). The boundary condition is evaluated within the solver iteration loop and constrains the exiting flow to only those 

seepage-face cells that are saturated. This seepage-face flow is assumed to exit the domain as surface runoff and does not re-enter 

the system. Simulation example case 4 demonstrates this new boundary condition during transient drainage of a fine sandy soil. 

Program Description

SPF1 subroutines are constructed to fit the standard MODFLOW process-package-procedure framework (fig. 4);

VSF1SPF1AL allocates space for the seepage-face location array; VSF1SPF1RP reads the potential seepage-face locations for the 

entire simulation and the output control options; VSF1SPF1FL resets the IBOUND values for seepage-face cells back to variable 

head (IBOUND=1) at the start of each stress period; VSF1SPF1FM activates and deactivates the seepage-face cells on the basis 

of the most recent total head values using the algorithm described in the previous section; VSF1SPF1BD records the seepage-face

VSF1SPF1AL

VSF1SPF1RP

VSF1SPF1FL

VSF1SPF1FM

VSF1SPF1BD

Stress-period

loop

Time-step

loop

Iteration

loop

MODFLOW–2000 Main

Figure 4. Subroutines used by the Seepage Face (SPF1) 

Package. A typical simulation proceeds from top to 

bottom. 
Seepage Face (SPF1) Package 13

locations and fluxes if requested. Flow from active seepage-face cells (IBOUND=-1) is treated as constant-head flow and 

accounted for by the Budget Procedure within the LPF1 or BCF6 Packages. The SPF1 subroutines are written in standard 

FORTRAN90 and should be compatible with any FORTRAN90 compiler.

Input Instructions

Input for the SPF1 Package is read from a file specified in the MODFLOW name file, identified as the ‘SPF’ file type. Free 

format is used for reading all variables. All boundary condition packages (for example, RCH6, CHD6) sharing the same cell locations will be overridden by the IBOUND manipulation within SPF1.

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. ISPFCB

Record 2 is only read if 

2. ISPFOC

Record 3 is only read if ISPFOC > 0

3. SPFOC (ISPFOC, 2)

4. SEEP (NCOL, NROW, NLAY) - U2DINT.

Explanation of the variables read by the SPF1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1.

1. ISPFCB—The flag containing the file unit for saving seepage-face locations and flows.

> 0, requested output will be saved to an ASCII file with the ‘.SFO’ extension located at unit ISPFCB.

= 0, no output.

< 0, requested output will be saved to a BINARY file with the ‘.SFO’ extension located at unit ISPFCB. 

2. ISPFOC—The flag for SPF1 output control. 

> 0, output recorded at times specified in the SPFOC array.

< 0, output recorded at all times.

3. SPFOC (ISPFOC, 2) —The 2-D output control array containing real values corresponding to the stress period and time step 

or total simulation time. There are ISPFOC entries corresponding to the total number of requested output times. The first 

value in each entry contains the requested stress period. The second value contains either a negative integer corresponding 

to the time-step number or a positive real number corresponding to the total simulation time. Read only if ISPFOC > 0.

Note: If the adaptive time stepping (ATS1) is active, the requested output must correspond to output times requested by the 

Output Control (OC) Option.

4. SEEP (NCOL, NROW, NLAY)—The integer array containing the locations of all potential seepage-face cells.

> 0 denotes a potential seepage-face cell.

< 0 denotes a normal MODFLOW cell.

Output Data

Simulation results for active (that is, exfiltrating) seepage-face cells are recorded to a ‘.SFO’ file for user-requested times. If 

ASCII output is requested, cell index locations (column, row, layer) and total flow within the cell are recorded under a header specifying the stress period, time step, elapsed simulation time, and total number of active seepage-face cells. For BINARY output, 

the MODFLOW utility UBUDSV is used to record a 3-D array of seepage-face fluxes for all grid cells at the requested simulation 

time. Flow from active seepage cells is accounted in the model budget as constant-head flow by the SVSF1BCF6F and 

SVSF1LPF1F subroutines and stored in the model’s cell-by-cell flow output file as requested in the OC Option. Example output 

can be found in Appendix 2.

IREFCB 0 ≠
14 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Surface Ponding (PND1) Package

Conceptualization

The Surface Ponding (PND1) Package enables MODFLOW to control surface ponding appropriately by restricting unrealistic 

hydraulic-head buildup in surface grid cells caused by exfiltration or large specified precipitation rates on dry soil. This approach 

is identical to the algorithm described by Lappala and others (1987). Simulation example cases 5 and 6 demonstrate this new 

boundary condition package for 2- and 3-D domains, respectively, with variable specified surface fluxes.

The surface ponding boundary condition interacts with the existing Recharge (RCH6) Package to switch the constant-flux 

boundary condition to a constant head (IBOUND=-1) if the hydraulic head in a surface cell exceeds a user-specified ponding 

elevation. The ponding elevation is determined by the sum of the cell’s top elevation and a ponding depth (POND) prescribed in 

the VSF input file. If the total flux from a ponded cell into the domain exceeds the specified surface flux by more than 1 percent, 

the boundary condition is switched back to constant flux. Ponded cells are accounted for in the model flow budget as constanthead flow; rainfall specified for those cells is neglected. In order to examine the behavior of the ponding algorithm under a range 

of hypothetical infiltration conditions, a series of simulations using a range of specified infiltration rates was performed for a 

medium-coarse soil for two initial conditions (fig. 5). The time required to produce ponding on the soil surface (POND = 0.0 cm) 

is plotted as a function of the rainfall rate scaled by the saturated hydraulic conductivity. It is shown that ponding time decreases 

rapidly with increasing rainfall rate and wetter initial conditions (ψ = -80 cm). 

Program Description

The subroutines are constructed to fit the standard MODFLOW process-package-procedure framework (fig. 6); 

VSF1PND1AL allocates space for the ponding depth array; VSF1PND1RP reads the ponding depths and the output control 

options; VSF1PND1FL resets the IBOUND values for ponded cells back to variable head (IBOUND=1) at the start of each stress 

period; VSF1PND1RST checks for ponded cells using the algorithm described in the previous section and restarts the iteration 

loop if ponding occurs; and VSF1PND1BD records the ponded cell locations and fluxes if requested. Flow from ponded cells 

Relative rainfall rate, dimensionless

 Ponding time, in hours

0

10

20

30

40

50

60

70

80

90

1.0 1.5 2.0 2.5 3.0

Initial pressure = -200 centimeters

Initial pressure = -80 centimeters

Figure 5. Surface ponding of a medium-coarse soil as a function of relative rainfall rate for two initial pressure 

conditions.
Surface Ponding (PND1) Package 15

(IBOUND=-1) is treated as constant-head flow and accounted for by the Budget Procedure within the LPF1 or BCF6 Packages. 

The PND1 subroutines are written in standard FORTRAN90 and should be compatible with any FORTRAN90 compiler.

Input Instructions

Input for the PND1 Package is read from a file specified in the MODFLOW name file, identified as the ‘PND’ file type. Free 

format is used for reading all variables. All boundary condition packages (for example, RCH6, CHD6) sharing designated surfacecell locations will be overridden by the IBOUND manipulation within PND1.

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. IPNDCB

Record 2 is only read if 

2. IPNDOC

Record 3 is only read if IPNDOC > 0

3. PNDOC (IPNDOC, 2)

4. POND (NCOL, NROW) - U2DREL.

Explanation of the variables read by the PND1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1.

1. IPNDCB—The flag containing the file unit for saving ponded cell locations and flows.

> 0, requested output will be saved to an ASCII file with the ‘.PDO’ extension located at unit IPNDCB.

VSF1PND1AL

VSF1PND1RP

VSF1PND1FL

VSF1PND1BD

Stress-period

loop

Time-step

loop

Iteration

loop

MODFLOW–2000 Main

VSF1PND1RST:

Surface ponding? IBOUND = -1, 

Restart time step.

Y

Figure 6. Subroutines used by the Surface Ponding (PND1) 

Package. A typical simulation proceeds from top to bottom.

IPNDCB 0 ≠
16 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

= 0, no output.

< 0, requested output will be saved to a BINARY file with the ‘.PDO’ extension located at unit IPNDCB. 

2. IPNDOC—The flag for PND1 output control.

> 0, output recorded at times specified in the PNDOC array.

< 0, output recorded at all times.

3. PNDOC (IPNDOC, 2) —The 2-D output control array containing real values corresponding to the stress period and time 

step or total simulation time. There are IPNDOC entries corresponding to the total number of requested output times. The 

first value in each entry contains the requested stress period. The second value contains either a negative integer 

corresponding to either the time-step number or a positive real number corresponding to the total simulation time. Read 

only if IPNDOC > 0.

Note: If the adaptive time stepping (ATS) is active, the requested output must correspond to output times requested by the 

Output Control Option.

4. POND (NCOL, NROW)—The real array containing the ponding depths of the surface cells (layer 1).

Output Data

Simulation results for ponded cells are recorded to a file with the ‘.PDO’ extension for user-requested times. If ASCII output

is requested, the output is recorded as list of cell index locations (column, row) and total flow within the cell under a header specifying the stress period, time step, elapsed simulation time, and total number of ponded cells. For BINARY output, the 

MODFLOW utility subroutine UBDSV3 is used to record the flow for each ponded cell as a 2-D array for the surface layer. Flow 

from ponded cells is accounted for in the model flow budget as constant-head flow by the SVSF1BCF6F and SVSF1LPF1F subroutines and stored as constant-head flow in the model’s cell-by-cell flow output file as requested in the Output Control Option. 

Example output can be found in Appendix 2.

Surface Evaporation (SEV1) Package

Conceptualization

A two-stage relationship is used to represent evaporation from bare soil at the ground surface (that is, the top model layer) 

within the Surface Evaporation (SEV1) Package presented below. Cyclic variations (that is, hourly or daily) in soil evaporation 

can be simulated through the specification of multiple stress periods representing each stage of the cycle (for example, night and 

day). The evaporation rate is driven by the pressure potential gradient between the soil pores and the atmosphere, described by the 

equation (Lapalla and others, 1987):

, (10)

where 

EV = evaporation rate across surface boundary cell (evaporation out of cell is positive) [L/T];

K(ψ) = hydraulic conductivity of surface cell as a function of pore-water pressure head [L/T];

SRES = surface resistance factor (that is, crusting) [L-1];

ψa = atmospheric pressure potential [L]; and

ψs = pore-water pressure head of surface boundary cell [L].

The actual flux extracted from the boundary cell is determined by the value of EV; if EV < PEV (the potential evaporation for a 

given stress period), the flux is computed as EV times the surface area of the cell’s top face; if EV > PEV, the flux is set equal to 

the potential evaporation times the surface area of the cell’s top face.

The SEV1 Package requires three variables to be specified for each stress period—the potential soil evaporation rate (PEV), 

the approximate atmospheric pressure potential (ψa), and the surface resistance factor (SRES). The potential evaporation rate 

EV K( ) ψ SRES ψs ψa = ⋅ ⋅ ( ) –
Surface Evaporation (SEV1) Package 17

(PEV) may be estimated from wind and temperature conditions using established methods such as the Penman equation (Campbell, 

1977). The atmospheric pressure potential (ψa) may be calculated with the Kelvin equation but typically is estimated as -103 m 

or greater to alleviate numerical instability (Lapalla and others, 1987). The surface resistance factor (SRES) is the reciprocal of 

the flow distance between the cell node and the ground surface scaled by the ratio of the soil crust hydraulic conductivity to the 

boundary cell’s hydraulic conductivity, as described by the following equation (Lapalla and others, 1987):

, (11)

where

THCK = thickness of surface boundary cell [L];

Kc = saturated hydraulic conductivity of surface crust [L/T]; and

Ktop = saturated hydraulic conductivity of top boundary cell [L/T].

Simulation example case 5 (Run 2) demonstrates the performance of the SEV1 Package by comparing the cumulative surface evaporation volume with results from VS2DT (Healy, 1990).

Program Description

The subroutines are constructed to fit the standard MODFLOW process-package-procedure framework (fig. 7); 

VSF1SEV1AL allocates space for the potential evaporation and surface resistance factor arrays; VSF1SEV1RP reads the potential 

evaporation, atmospheric pressure potential, surface resistance factor, and output control options for each stress period; 

VSF1SEV1FM computes the evaporative flux for each surface cell from the current soil-pore pressure distribution at the surface 

boundary; and VSF1SEV1BD records the soil evaporation from the surface boundary, if requested, and stores the total in the flow

budget array, VBVL. The SEV1 subroutines are written in standard FORTRAN90 and should be compatible with any 

FORTRAN90 compiler.

SRES 2

THCK

---------------- Kc

Ktop

= ⋅ ---------

VSF1SEV1AL

VSF1SEV1RP

VSF1SEV1FM

VSF1SEV1BD

Stress-period

loop

Time-step

loop

Iteration

loop

MODFLOW−2000 Main

Figure 7. Subroutines used by the Surface Evaporation 

(SEV1) Package. A typical simulation proceeds from top 

to bottom.
18 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Input Instructions

Input for the SEV1 Package is read from a file specified in the MODFLOW name file, identified as the ‘SEV’ file type. Free 

format is used for reading all variables.

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. ISEVCB

Record 2 is only read if 

2. ISEVOC ISEVCL

Record 3 is only read if SEVOC > 0

3. SEVOC (ISEVOC, 2)

Record 4 is only read if ISEVCL > 0

4. SEVCL (ISEVCL, 2)

FOR EACH STRESS PERIOD

5. INPEV INHA INSRES INSEL

Record 6 is only read if INPEV > 0; otherwise, the previous stress period's values are used.

6. PEV (NCOL, NROW) U2DREL

Record 7 is only read if INHA > 0; otherwise, the previous stress period's values are used.

7. HA

Record 8 is only read if INSRES > 0; otherwise, the previous stress period's values are used.

8. SRES (NCOL, NROW) U2DREL

Record 9 is only read if INSEL > 0; otherwise, the previous stress period’s values are used.

9. SEL(NCOL, NROW) U2DINT

Explanation of the variables read by the SEV1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1.

1. ISEVCB—The flag containing the file unit for saving soil evaporation flow from the surface cells.

> 0, requested output will be saved to an ASCII file with the ‘.SEO’ extension located at 

unit ISEVCB.

= 0, no output.

< 0, requested output will be saved to a BINARY file with the ‘SEO’ extension located at unit ISEVCB. 

2. ISEVOC—The flag for SEV1 output control.

> 0, output recorded at times specified in the SEVOC array.

< 0, output recorded at all times.

ISEVCL – The flag for SEV1 requested cell locations for output to the ‘.SEO’ file.

> 0, output recorded at locations specified in the SEVCL array for requested times.

= 0, only header containing total evapotranspiration flow rate is recorded for requested times.

< 0, output recorded at all locations for requested times.

3. SEVOC (ISEVOC, 2) —The 2-D output control array containing real values corresponding to the stress period and time step 

or total simulation time. There are ISEVOC entries corresponding to the total number of requested output times. The first 

value in each entry contains the requested stress period. The second value contains either a negative integer corresponding 

to the time-step number or a positive real number corresponding to the total simulation time. Read only if ISEVOC > 0.

Note: If the adaptive time stepping (ATS1) is active, the requested output must correspond to output times requested by the 

Output Control Option.

ISEVCB 0 ≠
Root Zone Evapotranspiration (RZE1) Package 19

4. SEVCL (ISEVCL, 2)—The 2-D cell location array containing integer values corresponding to the cell column and row 

numbers. There are ISEVCL entries corresponding to the total number of requested locations. The first value in each entry 

contains the requested cell’s column number, j. The second value contains the requested cell’s row number, i. Read only if 

ISEVCL > 0.

Note: This setting only applies to the ASCII formatted output file. For BINARY output (ISEVCB<0), results for the entire 

cell grid are recorded by the UBUDSV utility subroutine.

5. INPEV—The flag indicating whether to re-use the previous stress period’s value for potential evaporation (PEV).

< 0, the previous stress period’s value for PEV is used.

> 0, PEV is read from the SEV1 input file.

INHA—The flag indicating whether to re-use the previous stress period’s value for atmospheric pressure potential (HA).

< 0, the previous stress period’s value for HA is used.

> 0, HA is read from the SEV1 input file.

INSRES—The flag indicating whether to re-use the previous stress period’s value for surface resistance

factor (SRES).

< 0, the previous stress period’s value for SRES is used.

> 0, SRES is read from the SEV1 input file.

INSEL—The flag indicating whether to re-use the previous stress period’s value specifying the model layer from which 

soil evaporation is extracted (SEL).

< 0, the previous stress period’s value for SEL is used.

> 0, SEL is read from the SEV1 input file.

6. PEV (NCOL, NROW) —The real array containing the potential soil evaporation rate for the top boundary cells. Only read 

if INPEV > 0.

7. HA—The atmospheric pressure potential for the current stress period. Only read if INHA > 0.

8. SRES (NCOL, NROW) —The real array containing the surface resistance factors for the top boundary cells. Only read if 

INSRES > 0.

9. SEL (NCOL, NROW) —The integer array indicating the model layer from which soil evaporation is extracted. This array is 

typically set as layer 1; however, to avoid interference from Surface Ponding (PND1) or other packages sharing the top 

layer, the user may want to specify another model layer. Only read if INSEL > 0.

Output Data

Simulation results for evaporation flux are recorded to a file with the ‘.SEO’ extension for user-requested times. If ASCII 

output is requested, cell index locations (column, row, layer) and evaporation flux from each surface cell are recorded under a

header specifying the stress period, time step, and total simulation time. For BINARY output, the MODFLOW utility UBDSV3 

is used to record the flow as a 2-D array for the surface layer. The soil evaporation totals also are computed at the end of each time 

step and included in MODFLOW’s model budget. Example output can be found in Appendix 2.

Root Zone Evapotranspiration (RZE1) Package

Conceptualization

Evapotranspiration loss caused by plant-root extraction of surface vegetation in a given near-surface model cell can be 

described by the following equation (Lapalla and others, 1987):

, (12)

where

QET = evapotranspiration flux in cell (i,j,k), positive for evapotranspiration out of cell [L3/T];

QET

K( ) ψ rzt ( ) , ψijk , , ψroot ( ) – VOLijk , , ⋅ ⋅ ⋅ for ψijk , , ψroot and dijk , , d > < root

0 for ψijk , , ψroot or dijk , , d ⎩ ≤ ≥ root

⎨

⎧ =
20 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

r(z,t) = root activity function of depth and time [L-2];

ψroot = pressure head equal to the permanent wilting point of the plant roots [L];

ψi,j,k = pore-water pressure head in cell (i,j,k) [L];

VOLi,j,k = volume of cell (i,j,k) [L3];

di,j,k = depth of cell (i,j,k) [L]; and

droot = depth of root zone at horizontal location (i,j) [L].

The Root Zone Evapotranspiration (RZE1) Package computes fluxes within a given row i and column j location by totaling the 

evapotranspiration flux for each cell down the layers (k=1 to n) and comparing the total to the volumetric potential evaporation for 

the current stress period. The volumetric potential evapotranspiration is equal to the potential evapotranspiration rate (PET) times 

the cell’s horizontal surface area (A). If the computed total flux within the horizontal location (i,j) is greater than the PET•A value, 

each of the fluxes for that location are multiplied by the ratio of PET/QET* to adjust the total to equal PET•A, where QET* is the 

sum of QET for all cells at that horizontal location. It is implicit in this approach that cells with nodal elevations below the root 

zone depth, droot, are not included in the total volumetric evapotranspiration (equation 12).

There are six variables that must be specified for each stress period. Similar to the Soil Evaporation (SEV1) Package, cyclic 

variations (that is, hourly or daily) in evapotranspiration can be simulated through the specification of multiple stress periods representing each stage of the cycle (for example, night and day). The variables required include the potential evapotranspiration 

(PET) generally computed by the user from weather data with a function such as the Penman equation (Jensen and others, 1990), 

the depth of the bottom of the root zone (RTDPTH), the permanent wilting point of the plant roots (HROOT) as a function of horizontal location (typically equal to -150 m), the rooting depth (RTDPTH) below which roots cannot extract water from the soil, 

and the root activity at the bottom of the roots (RTBOT) and at the land surface (RTTOP). The root activity generally is defined 

as the root length in a given soil volume divided by that soil volume and is linearly interpolated between RTBOT and RTTOP.

Program Description

The subroutines are constructed to fit the standard MODFLOW process-package-procedure framework (fig. 8); 

VSF1RZE1AL allocates space for the potential evapotranspiration, wilting point, root depth, and root activity arrays;

VSF1RZE1AL

VSF1RZE1RP

VSF1RZE1FM

VSF1RZE1BD

Stress-period

loop

Time-step

loop

Iteration

loop

MODFLOW–2000 Main

Figure 8. Subroutines used by the Root Zone Evapotranspiration (RZE1) Package. A typical simulation proceeds 

from top to bottom.
Root Zone Evapotranspiration (RZE1) Package 21

VSF1RZE1RP reads the potential evapotranspiration, root depth, wilting point, top root activity, bottom root activity, and output 

control options for each stress period; VSF1RZE1FM computes the evapotranspiration flux for the root zone at each horizontal 

location from the most recent computed soil-pore pressure distribution; and VSF1RZE1BD records the evapotranspiration from 

the root zone, if requested, and stores the total in the flow budget array, VBVL. The RZE1 subroutines are written in standard

FORTRAN90 and should be compatible with any FORTRAN90 compiler.

Input Instructions

Input for the RZE1 Package is read from a file specified in the MODFLOW name file, identified as the ‘RZE’ file type. Free 

format is used for reading all variables. 

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. IRZECB

Record 2 is only read if 

2. IRZEOC IRZECL

Record 3 is only read if IRZEOC > 0

3. RZEOC (IRZEOC, 2)

Record 4 is only read if IRZECL > 0

4. RZECL (IRZECL, 2)

FOR EACH STRESS PERIOD

5. INPET INRTDPTH INRTBOT INRTTOP INHROOT INRZL

Record 6 is only read if INPET > 0; otherwise, the previous stress period’s values are used.

6. PET (NCOL,NROW) U2DREL

Record 7 is only read if INRTDPTH > 0; otherwise, the previous stress period’s values are used.

7. RTDPTH (NCOL, NROW) U2DREL

Record 8 is only read if INRTBOT > 0; otherwise, the previous stress period’s values are used.

8. RTBOT (NCOL, NROW) U2DREL

Record 9 is only read if INRTTOP > 0; otherwise, the previous stress period’s values are used.

9. RTTOP (NCOL, NROW) U2DREL

Record 10 is only read if INHROOT > 0; otherwise, the previous stress period’s values are used.

10. HROOT (NCOL, NROW) U2DREL

Record 11 is only read if INRZL > 0; otherwise, the previous stress period’s values are used.

11. RZL (NCOL, NROW) U2INT

Explanation of the variables read by the RZE1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1.

1. IRZECB—The flag containing the file unit for saving evapotranspiration flows and cell locations.

IRZECB 0 ≠
22 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

> 0, requested output will be saved to an ASCII file with the ‘.RZO’ extension located at 

unit IRZECB.

= 0, no output.

< 0, requested output will be saved to a BINARY file with the ‘.RZO’ extension located at unit IRZECB. 

2. IRZEOC—The flag for RZE1 output control.

> 0, output recorded at times specified in the RZEOC array.

< 0, output recorded at all times.

IRZECL—The flag for RZE1 requested cell locations for output to the ‘.RZO’ file.

> 0, output recorded at locations specified in the RZECL array for requested times.

= 0, only header containing total evapotranspiration flow rate is recorded for requested times.

< 0, output recorded at all locations for requested times.

3. RZEOC (IRZEOC, 2) —The 2-D output control array containing real values corresponding to the stress period and time step 

or total simulation time. There are IRZEOC entries corresponding to the total number of requested output times. The first 

value in each entry contains the requested stress period. The second value contains either a negative integer corresponding 

to the time-step number or a positive real number corresponding to the total simulation time. Read only if IRZEOC > 0.

Note: If the adaptive time stepping (ATS1) is active, the requested output must correspond to output times requested by the 

Output Control Option.

4. RZECL (IRZECL, 2) —The 2-D output cell location array containing integer values corresponding to the cell column and 

row numbers (j,i). There are IRZECL entries corresponding to the total number of requested locations. The first value in 

each entry contains the requested cell’s column number, j. The second value contains the requested cell’s row number, i. 

Read only if IRZECL > 0.

Note: This setting only applies to the ASCII formatted output file. For BINARY output (IRZECB<0), results for the entire 

cell grid are recorded by the UBUDSV utility subroutine.

5. INPET—The flag indicating whether to re-use the previous stress period’s value for potential evapotranspiration (PET).

< 0, the previous stress period’s value for PET is used.

> 0, PET is read from the RZE1 input file.

INRTDPTH–The flag indicating whether to re-use the previous stress period’s value for root zone depth (RTDPTH).

< 0, the previous stress period’s value for RTDPTH is used.

> 0, RTDPTH is read from the RZE1 input file.

INRTBOT—The flag indicating whether to re-use the previous stress period’s value for root activity at the bottom of the 

root zone (RTBOT).

< 0, the previous stress period’s value for RTBOT is used.

> 0, RTBOT is read from the RZE1 input file.

INRTTOP—The flag indicating whether to re-use the previous stress period’s value for root activity at the top of the root 

zone (RTTOP).

< 0, the previous stress period’s value for RTTOP is used.

> 0, RTTOP is read from the RZE1 input file.

INHROOT—The flag indicating whether to re-use the previous stress period’s value for the plants’ permanent wilting point 

(HROOT).

< 0, the previous stress period’s value for HROOT is used.

> 0, HROOT is read from the RZE1 input file.

INRZL—The flag indicating whether to re-use the previous stress period’s value specifying the topmost model layer from 

which to extract root zone evapotranspiration (RZL).

< 0, the previous stress period’s value for RZL is used.

> 0, RZL is read from the RZE1 input file.

6. PET (NCOL, NROW) —The array containing the potential soil evapotranspiration rate for the cells within the root zone. 

Only read if INPET > 0.

7. RTDPTH (NCOL, NROW) —The array containing the depth of the root zone from the ground surface. Only read if 

INRTDPTH > 0.

8. RTBOT (NCOL, NROW) —The array containing the root activity at the bottom of the root zone. Only read if 

INRTBOT > 0.

9. RTTOP (NCOL, NROW) —The array containing the root activity at the top of the root zone. Only read if INRTTOP > 0.
Adaptive Time-Stepping (ATS1) Package 23

10. HROOT (NCOL, NROW) —The array containing the pressure equal to the permanent wilting point of the plants. Only read 

if INHROOT > 0.

11. RZL (NCOL, NROW) —The integer array indicating the top model layer of the stack of cells within the root zone from 

which evapotranspiration is extracted. This array is typically set as layer 1; however, to avoid interference from the Surface

Ponding (PND1) Package or other packages sharing the top layer, the user may want to specify another model layer. Only 

read if INRZL > 0.

Output Data

Simulation results for evapotranspiration flux are recorded to a file with the ‘.RZO’ extension for user-requested times. If 

ASCII output is requested, cell index locations (column, row, layer) and evapotranspiration flux from each cell are recorded under 

a header specifying the stress period, time step, and total simulation time. For BINARY output, the MODFLOW utility UBUDSV 

is used to record a 3-D array of evapotranspiration fluxes for all grid cells at the requested simulation time. The evapotranspiration 

totals also are computed at the end of each time step and included in MODFLOW’s model budget. Example output can be found 

in Appendix 2.

Adaptive Time-Stepping (ATS1) Package

Conceptualization

Transient flow simulations of nonlinear equations often have difficulty converging in MODFLOW (for both saturated and 

unsaturated simulations) due to a user-prescribed time step that is too long. Conversely, for long simulations of multiple stress 

periods with variable boundary conditions, prescribed time steps that are too short will cause the overall computation time to be 

unnecessarily long. MODFLOW currently (2005) does not have a time-stepping scheme that is able to adjust the time step during 

a simulation on the basis of computational conditions. The Adaptive Time-Stepping (ATS1) Package was developed to increase 

the model’s efficiency during transient simulations (fig. 9), similar to the algorithms implemented in the U.S. Salinity Laboratory 

code SWMS_3D (Simunek and others, 1995) and the commercial code MODFLOW-SURFACT (Hydrogeologic, 1997). 

The new scheme requires the user to input the maximum time-step size (TMX), the minimum time-step size (TMN), a timestep reduction factor (TSD), and a time-step multiplier (TSM), as well as a file number (ITIM) and array size (MXSTP) for recording the time-step history of a simulation in an output file (for interfacing with a transport model). The original MODFLOW input 

files are still used to set the number of time steps (NSTP) for each stress period and to determine the time-step size (DELT). During 

the simulation, after the computations for each time step are completed, the new scheme evaluates the number of iterations required 

for convergence. If convergence is achieved in more than 65 percent of the user-defined maximum number of iterations 

(MXITER), the time step is divided by TSD in anticipation of convergence difficulties in the next time step. If convergence is

achieved in less than 35 percent of MXITER, the time step is multiplied by TSM to speed up the simulation. If the solution fails 

to converge within a given time step, the time step is reduced by TSD, and the iterative computations are restarted with the reduced 

time step. If the adjusted time step exceeds the maximum size (TMX), the step is set equal to TMX. Conversely, if the time step 

is reduced to less than the minimum size (TMN), the time step is set equal to TMN. The counter variable that tracks elapsed time 

in the current stress period (PERTIM) is adjusted using the latest time step (DELT). The last time step of each stress period is 

adjusted as necessary to prevent exceeding the user-defined stress period length (PERLEN). Additionally, to maintain compatibility with the Output Control Option, the internal MODFLOW time-step counter (KSTP) is adjusted, and time-step adjustments 

are constrained as the solution proceeds to ensure that requested output corresponds correctly with the specified time steps and 

stress periods. 

As an example of ATS1 performance, the time-step adjustment from simulation day 333 of simulation example case 5 is listed 

in table 1. The maximum number of solver iterations was set to 10, which caused the adaptive time-stepping algorithm to reduce 

the time-step size if the iteration ratio was greater than 0.65 (iterations > 6) and increase the time-step size if the iteration ratio was 

less than 0.35 (iterations < 4).

Program Description

The ATS1 Package is constructed to fit within the modular framework of MODFLOW and is compatible with both the GWF 

Process as well as the new VSF Process. A flowchart of its operation and interactions with the GWF and VSF Processes is shown 

in figure 10. ATS1AL (Allocate Procedure) reads in the time-stepping parameters and allocates array space (if needed) for the 
24 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Table 1. Adaptive time-stepping results for day 333 of simulation case 5.

[--, not applicable]

Elapsed time (days) Solver iterations Iteration ratio Time-step size

0 -- -- 0.041666

.041667 7 0.7 .020833

.0625 4 .4 .020833

.083333 3 .3 .03125

.114583 4 .4 .03125

.145833 3 .3 .046875

.192708 4 .4 .046875

.239583 4 .4 .046875

.286458 3 .3 .070313

.356771 4 .4 .070313

.427083 4 .4 .070313

.497396 3 .3 .105469

.602865 4 .4 .105469

.708333 4 .4 .105469

.813802 4 .4 .105469

.919271 4 .4 .080729

1.0 -- -- --

End of time-step loop

Restart time step

No: 

Reduce 

time 

step

Solution 

convergence?

Too many 

iterations?

Yes: 

Continue

Continue to next 

time step

Yes: 

Reduce 

time step 

gradually

No: 

Continue

Too few 

iterations?

No: 

Continue

Yes: Increase time step 

gradually

Figure 9. Algorithm for the Adaptive Time-Stepping (ATS1) Package.
Adaptive Time-Stepping (ATS1) Package 25

time-step history; ATS1ADV (Advance Procedure) advances the time-stepping counters; ATS1RST (Restart Procedure) reduces 

the time-step size and restarts the time step if the solution did not converge; otherwise, the time-step size is stored in an array, and 

the iteration loop is exited; and ATS1ADJ (Adjust Procedure) adjusts the time-step size on the basis of the iteration history of the 

previous time step and if the end of the stress period has been reached, the time-step history is saved to a file (if requested).

Input Instructions

Input for the ATS1 Package is read from a file specified in the MODFLOW name file, identified as the ‘ATS’ file type. Free 

format is used for reading all variables. 

FOR EACH SIMULATION

0. [#Text]

Item 0 is optional—“#” must be in column 1. Item 0 can be repeated multiple times.

1. ITIM MXSTP 

2. TMX TMN TSM TSD

Explanation of the variables read by the ATS1 Package

Text—is a character variable (199 characters) that starts in column 2. Any characters can be included in Text. The “#” character 

must be in column 1. 

1. ITIM—The flag containing the file unit for saving time-step sizes in a separate file.

N

Iteration

loop

Time-step

loop

Stress-period

loop

VSF Process:

• BASIC

• Layer Property Flow *

or Block Centered Flow*

• Recharge

• Soil Evaporation

• Root Zone Evapotranspiration

• Surface Ponding 

• Seepage Face

*Modified for use with VSF Process

Other packages valid in

fully saturated cells 

Global Process

VSF?

Normal MODFLOW 

Operation:

• Ground

Flow 

Process

• Observation, 

Sensitivity, and 

Parameter 

Estimation 

Processes

Y

ATS Package Procedures

ATS: Allocate

ATS: Advance

ATS: Adjust

ATS: Restart

Figure 10. Procedures used by the Adaptive Time-Stepping (ATS1) Package. A typical simulation proceeds from top to 

bottom.
26 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

> 0, output will be saved to unit ITIM in a binary file with the ‘.TIM’ extension.

< 0, no output.

MXSTP—Maximum number of time steps (during a given stress period) expected for entire simulation.

Used to allocate the time-step array that is saved to the output (‘.TIM’) file.

2. TMX—Maximum allowable time-step size.

TMN—Minimum allowable time-step size.

TSM—Time-step multiplication factor.

TSD—Time-step reduction (division) factor.

Output Data

If requested, the time-step sizes for the entire simulation are recorded to a binary file with the ‘.TIM’ extension for post-processing. The output is recorded as a header containing the period number, period length, and number of actual time steps used 

(after adjustment). Under the header, an unformatted array of time-step sizes is recorded for every time step of the simulation. An 

example input file and description of the output file format are included in Appendix 2.

Simulation Examples

Six simulation examples are presented here to test the new packages described in the previous section. Model performance 

is evaluated by comparison to an analytical solution for 1-D constant-head infiltration (Dirichlet boundary condition) in case 1, 

field experimental data for a 1-D constant-head infiltration in case 2, laboratory experimental data for 2-D constant-flux infiltration 

(Neumann boundary condition) in case 3, laboratory experimental data for 2-D transient drainage through a seepage face in case 

4, and numerical model results (VS2DT) of a 2-D flowpath simulation using realistic surface boundary conditions in case 5. A 

hypothetical 3-D example case also is presented in case 6 to demonstrate the new capability using periodic boundary conditions 

(for example, daily precipitation) and varied surface topography over a larger spatial scale (0.133 km2).

Case 1—One-Dimensional Analytical Solution for Unsaturated Flow

The strong nonlinearity of Richards’ equation precludes a closed-form analytical solution of the flow equation except under 

very restrictive initial and boundary conditions. Phillip (1969) presented a quasi-analytical solution for 1-D infiltration in a semiinfinite column under a constant-head boundary condition. It is used here to verify the accuracy of this model’s formulation.

The Philip solution takes the form of an infinite series (valid for finite values of t):

, (13)

where

z(θ,t) = depth [L];

θ = volumetric soil moisture content (dimensionless);

t = time [T]; and

φn = nth order function of θ.

This solution is valid for the following initial and boundary conditions:

To evaluate the performance of the unsaturated flow component of the VSF Process, two 1-D simulations were performed 

using the initial and boundary conditions just listed. These simulations represent infiltration in a sand column and clay column 

and are based on the work of Haverkamp and others (1977). The initial and boundary conditions as well as the characteristic relationships (K(ψ), Θ(ψ), and C(ψ)) used are listed in table 2. Note that custom soil characteristic functions are used for this example 

t = 0, z > 0, θ = θ0

t > 0, z = 0, θ = θ1.

z( ) θ, t φ1t

1/2 φ2t φ3t3/2 φ4t2 = ++ + + ...
Simulation Examples 27

to obtain the best results; see Appendix 3 for a discussion on using custom functions. Model results for both simulations are plotted 

in relation to the Phillip solution at selected times in figure 11. In both the sand and the clay cases, the model matches the analytical 

solution values very closely despite some minor underprediction of the infiltration front for sand (fig. 11B) during the later time 

(t = 0.8 hour).

Case 2 —One-Dimensional Field Study

Experimental results from the infiltration study of Warrick and others (1971) are used here to demonstrate the utility of the 

VSF Process in representing a real physical system. This field experiment was performed in a 6.1-m2

 soil plot containing tensiometers at 30-, 60-, 90-, 120-, 150-, 180-cm depths. A conservative tracer solution (0.2 N CaCl2) was used to track solute transport 

through the system. Soil tensiometer measurements were collected for 17.5 hours during which the soil surface was kept at constant 

saturation. 

The initial soil-moisture profile is approximately linear from Θ = 0.394 cm3/cm3 at the surface to Θ = 0.525 cm3/cm3 at 

z < 60 cm. The soil characteristic functions were fit3 to the soil hydraulic conductivity and capillary saturation data. A summary 

of the soil parameters and model conditions used for this VSF simulation are listed in table 3. Results for the flow simulation after 

2 hours and 9 hours of infiltration are shown in figure 12 plotted in relation to the experimental data as well as simulation results 

from the 1-D unsaturated flow finite-element model, SUMATRA (van Genuchten, 1978). The VSF results represent the 

experimental results well; there is some overprediction of the infiltration front at t = 2 hours that can be attributed to the overprediction of the experimental values for kr at low Θ in the van Genuchten relationship.

Case 3 —Two-Dimensional Variably Saturated Transient Infiltration Study

The VSF Process was tested against laboratory measurements from the 2-D transient water-table experiment of Vauclin and 

others (1979), as presented by Clement and others (1994) and VanderKwaak (1999). The experiment consisted of a 6.0-m by 

2.0-m box containing a sandy soil with the initial water table located at 0.65 m from the bottom. A constant flux of 

q = 3.55 m/d was applied across the center 1.0 m of the soil surface for 8 hours, and the remainder of the surface was covered to 

prevent evaporative losses. 

Due to the experiment’s symmetry, only one-half of the box was modeled. The model domain was 3.0 m by 2.0 m, with the 

constant-flux boundary condition (q = 3.55 m/d) applied across the left 0.5 m of the defined ground surface. The initial total head 

of all cells was set to 0.65 m, and the right boundary cells were constrained to this initial water-table position throughout the 

8-hour simulation. The model parameters for the VSF simulation are listed in table 4. The modeled water-table positions for 

selected times were compared with the experimental data of Vauclin and others (1979) in figure 13 and showed that the model 

satisfactorily captured the transient water-table dynamics.

Case 4—Two-Dimensional Variably Saturated Transient Drainage Study

The experiment of Vauclin and others (1975) examined transient drainage in a fine sandy soil through a seepage face within 

a rectangular 3.0-m by 2.0-m by 2-cm thick soil slab. The transient positions of the water table measured by Vauclin and others

(1975) (as presented by Clement and others, 1994) were used to test the seepage-face boundary condition described in the VSF 

Process. The soil slab initially was held in hydrostatic equilibrium at h = 1.45 m. The experiment involved measuring the change 

in moisture content and pressure head as the constant-head boundary imposed on one side of the soil slab instantaneously dropped 

to h = 0.75 m. The measurements continued until the water-table position equilibrated to the boundary condition of 0.75 m.

The VSF model domain was constrained by a constant head of 0.75 m on the right boundary cell at and below the water table. 

Above the water table, the right boundary cells were set as a seepage face. The first 5 hours of the experiment were simulated. 

The model parameters are listed in table 5. Note that custom soil characteristic functions were used to obtain the best results; see 

Appendix 3 for details on using custom functions.

3

 The interpolation utility RETC (van Genuchten and others, 1997) was used to determine the best-fit van Genuchten-Mualem soil parameters for this soil.
28 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW Table 2. Model parameters for constant-head infiltration in a soil column.

[Kr(ψ), relative hydraulic conductivity, dimensionless; 

ψ, pore-water pressure, in centimeters; 

Ksat, saturated hydraulic conductivity, in centimeters per hour; 

Θ(ψ), soil saturation as a function of pore-water 

pressure, dimensionless; 

Θr, residual soil saturation, dimensionless; ne, porosity, dimensionless; 

C(ψ), specific moisture capacity, in centimeters-1; t, simulation time, in seconds; z, depth from surface, in centimeters; cm3/cm3, cubic centimeter per cubic centimeter;

Δx, cell dimension, in centimeters; cm, centimeter]

Yolo light clay (Haverkamp and others, 1977) Sand (Haverkamp and others, 1977)

Kr(ψ)

;

Kr(ψ)

;

Θ(ψ)

;

Θ(ψ)

;

C(ψ) C(ψ)

Initial conditions t = 0, z > 0, 

Θ0 = 0.48 cm3/cm3 Initial conditions t = 0, z > 0, 

Θ0 = 0.348 cm3/cm3

Boundary conditions t > 0, z = 0, 

Θ = 1.0 cm3/cm3 Boundary conditions t > 0, z = 0, 

Θ = 0.930 cm3/cm3

Discretization

Δx = 1 cm;

Δt = 40 seconds (t < 130 hours);

Δt = 500 seconds (t > 130 hours)

Discretization

Δx = 1 cm;

Δt = 5 seconds

kr 124.6

124.6

ψ

1.77

+

----------------------------------

=

Ksat 4.428

E 2 cm/h

–

=

kr 1.175

E 6 + ( )

1.175

E 6 + ( )

ψ

4.74

+

----------------------------------------------------

=

Ksat 34 cm/h

=

Θ 739 1

Θr – ( )

739

ψ4 ln

+

------------------------------ Θr +

=

Θr 0.251

=

Θ 1.611

E 6 + ( ) 1 Θr – ( ) ⋅

1.611

E 6 + ( )

ψ

3.96

+

-------------------------------------------------------- Θr +

=

Θr 0.261

=

C 739 4⋅ ( )

ψ

-------------------- ne ne Θr ⋅ – ( )

ψ3 ln ⋅

739

ψ4 ln

+ [ ]2

--------------------------------------------------; =

ne 0.495

=

C 3.96 1.611

E 6 + ⋅ ( )

ψ

---------------------------------------------- ne ne Θr ⋅ – ( )

ψ

2.96 ⋅

1.611

E 6 ψ3.96 + + [ ]2

------------------------------------------------------

=

ne 0.287

=

;
Simulation Examples 29

(A) Yolo light clay

(B) Sand

0

10

20

30

40

50

60

70

80

90

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

t = 0.1 hour

t = 0.2 hour

t = 0.8 hour

Model 

results

Analytical 

solution

Saturation, Θ, dimensionless

0

20

40

60

80

100

120

140

160

180

200

0.35 0.45 0.55 0.65 0.75 0.85 0.95 1.05

Depth, in centimeters

t = 28.8 hours

t = 278.8 hours

t = 833.3 hours

Model 

results

Analytical 

solution

Figure 11. Comparison between simulated water saturation (Θ) profiles and Phillip’s analytical solution for

one-dimensional, constant-head infiltration in (A) Yolo light clay and (B) sand.
30 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Table 3. Model parameters for the Warrick and others (1971) field study.

[t, simulation time, in seconds; z, depth from surface, in centimeters; cm, centimeters; Θ, soil saturation, dimensionless; Θ0 , initial soil saturation, dimensionless; Θ1, constant soil saturation boundary condition, dimensionless; Ksat , saturated hydraulic conductivity, in centimeters per day; cm/d, centimeters per day; ne porosity, dimensionless; Θr , residual soil saturation, dimensionless; α, soil parameter for van Genuchten-Mualem relationships, in 

centimeters-1; cm-1, centimeters-1; nv, soil parameter for van Genuchten-Mualem relationships, dimensionless]

Parameter Parameter value

Model conditions

Initial conditions t = 0, 0 < z < 60 cm, 

Θ0= 0.394 + 0.00219*z

t < 0, z > 60 cm, Θ0 = 0.525 

Boundary conditions t > 0, z = 0, Θ = 0.997 

Soil properties 

Medium coarse soil Ksat = 38.16 cm/d

ne = 0.381, Θr = 0.394

van Genuchten-Mualem parameters:

α = 0.016 cm-1, nv = 2.7

Saturation, Θ, dimensionless

Depth, in centimeters

Initial profile

t = 2 hours

t = 9 hours

0

20

40

60

80

100

120

140

160

180

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

VSF results

Experimental data

SUMATRA results

Figure 12. Model results from the one-dimensional simulation plotted in relation to experimental data 

and the results of a comparable model. 
Simulation Examples 31

Table 4. Model parameters for the Vauclin and others (1979) infiltration study.

[h, hydraulic head, in meters; z, elevation, in meters; m, meters; x, horizontal distance, in meters; q, constant surface flux, in meters per day; m/d, 

meters per day; Δx, cell width, in meters; Δz, cell thickness, in meters; Δt, simulation time-step size, in minutes; Ksat, saturated hydraulic conductivity, 

in meters per day; ne, porosity, dimensionless; Θr , residual soil saturation, dimensionless; α, soil parameter for van Genuchten-Mualem relationships, 

in meters-1; m-1, meters-1; nv, soil parameter for van Genuchten-Mualem relationships, dimensionless]

Parameter Parameter value

Model conditions

Initial conditions h = 0.65 m

Boundary conditions Surface cells: z = 2.0 m, 0 < x < 0.5 m,

q = 3.55 m/d

Right boundary cells: 

x = 3.0 m, 0 < z < 0.65 m,

h = 0.65 m

Discretization Δx = 0.1 m;

Δz = 0.05 m;

Δt = 1 minute

Soil properties

Sandy soil Ksat = 8.4 m/d

ne = 0.3, Θr = 0.033

van Genuchten-Mualem parameters:

α = 3.3 m-1, nv = 4.1

Initial water table

Model results

Experimental data

Constrained at 

0.65 meter

Distance, in meters

Elevation, in meters

t = 8 hours

t = 4 hours

t = 3 hours

t = 2 hours

q = 3.55 meters per day

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

0123

Figure 13. Modeled water-table positions for two-dimensional infiltration study plotted in relation to the 

experimental data of Vauclin and others (1979).
32 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

The modeled water-table positions for selected times were compared with the experimental data of Vauclin and others (1975) 

in figure 14. The model satisfactorily captured the seepage-face position along the boundary as the soil slab drained, indicating 

proper implementation of the boundary condition.

Case 5 —Two-Dimensional Variably Saturated Model Comparison

A 2-D comparison case simulating the seasonal dynamics along a hypothetical catchment flow path is presented here to evaluate the performance of the new VSF Process by comparison to a similar model, VS2DT (Healy, 1990). Two simulation runs are 

presented; run 1 uses the RZE1 Package to simulate root zone evapotranspiration during variable specified surface flux (that is, 

precipitation, via the RCH6 Package) with the PND1 Package simulating surface ponding, and run 2 uses the SEV1 Package to 

simulate surface evaporation under identical initial and boundary conditions. 

Table 5. Model parameters for the Vauclin and others (1975) drainage study.

[h, hydraulic head, in meters; x, horizontal distance, in meters; z, elevation, in meters; m, meters; Δx, cell width, in meters; Δz, cell thickness, 

in meters; Δt, simulation time-step size, in minutes; Kr (ψ), relative hydraulic conductivity as a function of pore-water pressure, dimensionless; ψ, 

pore-water pressure, in centimeters; Ksat, saturated hydraulic conductivity, in meters per hour; m/h, meters per hour; Θ(ψ), soil saturation as a function 

of pore-water pressure, dimensionless; ne, porosity, dimensionless; Θr , residual soil saturation, dimensionless]

Parameter Parameter value

Model conditions

Initial conditions h = 1.45 m

Boundary conditions Right boundary cells: 

x = 3.0 m, 0 < z < 0.75 m,

h = 0.75 m

x = 3.0 m, 0.75 m < z < 2.0 m,

Seepage face

Discretization Δx = 0.1 m;

Δz = 0.05 m;

Δt = 0.6 minute

Soil properties (Vauclin and others, 1975)

kr(ψ)

Θ(ψ)

kr

3.6E + 5

( ) 3.6E + 5 ψ

4.5

+ = --------------------------------------------;

Ksat = 0.4 m/h

Θ ( ) 4.0E + 4 1 Θr ⋅ ( ) –

( ) 4.0E + 4 ψ

2.9

+

-------------------------------------------------- Θr = + ;

Θr = 0.0
Simulation Examples 33

The model domain, shown in figure 15, consists of a 530-m long by 12.2-m deep slice with a 3.0-percent slope. At the downslope end of the domain, a 1-m deep constant-stage irrigation canal is represented by constant-head cells in the top seven layers of 

the model grid, constraining the head to the ground-surface elevation at that end. The soil media consisted of a homogeneous 

medium-coarse soil (Ksat = 1.0 m/d). The soil parameters, discretization, and boundary conditions are listed in table 6. A year of 

daily precipitation and potential evapotranspiration totals from Fort Collins, Colorado (National Climate Data Center, 2002), are 

used in the simulation (fig. 16). Both simulation runs use the same surface source totals (that is, precipitation) for the RCH6 input 

file and potential surface sink totals (that is, potential evapotranspiration or evaporation) for the RZE1 (run 1) and SEV1 (run 2) 

input files.

To minimize any flow bias caused by the initial head distribution, several consecutive simulation runs were made (using the 

heads from the end of the previous run as the initial heads) until the water table reached a general equilibrium. The overall hydraulic behavior is similar for run 1 and run 2; for brevity, only water-table, soil-saturation, and ponding results for run 1 are displayed, 

and only cumulative soil evaporation volume results are displayed for run 2. 

For run 1, the time series of the depth to the water table for two locations (355 and 180 m upslope from the irrigation canal, 

respectively) is plotted with precipitation in figure 17. The simulated soil saturation for run 1 at the bottom of the root zone 

(depth = 0.6 m) for the same two locations (180 m and 355 m upslope) is plotted in figure 18. The saturation results display similar 

behavior as the depth to water results; the soil saturation in the upslope location is slightly less responsive to big precipitation, 

whereas the downslope location is more affected by evapotranspiration losses during the summer months due to the shallower 

water table. The model results of run 1 match up well for both saturation and depth to water in the upslope location (figs. 17 and 

18); however for the downslope location, VSF slightly overestimates the VS2DT results for soil saturations and underestimates 

the depth to water results through most of the year. For run 2, the simulated cumulative soil evaporation volumes for both models 

are plotted in figure 19. The excellent match between the two models during a complex long-term simulation verifies the performance of the SEV1 and RZE1 Packages.

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

0123

Constrained at 

0.75 meter

Final water table

Initial water table

t = 0.1 hour

t = 0.5 hour

t = 1.0 hour

t = 2.0 hour

t = 5.0 hour

Distance, in meters

Elevation, in meters

Model results

Experimental data

Figure 14. Modeled water-table positions for two-dimensional drainage study plotted in relation to 

experimental data of Vauclin and others (1975).
34 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Table 6. Model parameters for the two-dimensional hypothetical comparison case.

[x, horizontal distance, in meters; z, elevation, in meters; m, meters; h, hydraulic head, in meters; Δx, cell width, in meters; Δy, cell length, in meters; 

Δz, cell thickness, in meters; Δt, simulation time-step size, in minutes; TMAX, maximum time-step size, in days; TMIN, minimum time-step size, in 

days; Ksat, saturated hydraulic conductivity, in meters per day; ne, porosity, dimensionless; Θr , residual soil saturation, dimensionless; α, soil parameter for van Genuchten-Mualem relationships, in meters-1; m-1, meters-1; m/d, meters per day; nv, soil parameter for van Genuchten-Mualem relationships, dimensionless]

Parameter Parameter value

Model conditions

Initial conditions Final heads of third simulation run

Boundary conditions Constant-head cells:

x = 530.0 m, 9.03 m < z < 8.03 m,

h = 9.03 m

BOTH RUNS:

 Temporally variable recharge with surface ponding

RUN 1: Root zone evapotranspiration

Extinction depth = 0.6 m

Root activity at surface = 9E-5 m-2

Root activity at 0.6 m depth = 2E-5 m-2

Pressure head in roots = -80 m

Temporally variable evapotranspiration rates 

RUN 2: Surface evaporation

Atmospheric pressure: -0.001 m

Surface resistance: 2.0E-5 m-1

Temporally variable surface evaporation rates

Discretization Δx = Δy = 5 m;

Δz = 0.1 to 6 m;

Adaptive time-stepping: 

Initial Δt = 1 hour

TMAX = 0.5 day, TMIN = 0.01 day

Soil properties

Medium-coarse soil Ksat = 1.0 m/d

ne = 0.381, Θr = 0.3937

van Genuchten-Mualem parameters:

α = 1.6 m-1, nv = 2.7

530 meters

12.2 meters

x

z

Constanthead cells

Slopex = 3.0 percent

Figure 15. Conceptual model for a catchment flow path that feeds a steady irrigation canal.
Simulation Examples 35

-0.010

-0.005

0

0.005

0.010

0.015

0.020

0.025

4565

105145

165

185

225

245

305325

345

365

25

125

205

265

285

85

5

Precipitation

Potential 

evapotranspiration

Daily total, in meters

Simulation day

Figure 16. Daily precipitation and potential evapotranspiration totals used for catchment flow-path 

simulations. 

Simulation day

Depth to water, in meters

Precipitation, in meters

4565

105145

165

185

225

245305

325

345

365

25

125

205

265285

85

5

180 meters upslope

355 meters upslope

VSF results VS2DT results Precipitation

0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0 0

0.005

0.010

0.015

0.020

0.025

EXPLANATION

Figure 17. Comparison of VSF and VS2DT model results for depth to water at locations 180 and 355 meters 

upslope from the irrigation canal, plotted with daily precipitation. The REF1, RZE1 and PND1 Packages were used 

in the VSF simulation for run 1 (shown here).
36 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Simulation day

Saturation, dimensionless

Precipitation, in meters

4565

105145

165

185225

245305

325

345

365

25

125

205

265285

85

5

180 meters upslope

355 meters upslope

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0

0.005

0.010

0.015

0.020

0.025

VSF results VS2DT results Precipitation

EXPLANATION

Figure 18. Comparison of VSF and VS2DT model results for soil saturation at locations 180 and 355 meters 

upslope from the irrigation canal, plotted with daily precipitation. The REF1, RZE1 and PND1 Packages were 

used in the VSF simulation for run 1 (shown here). 4565 105 145

165

185225245

305325

345365

25

125

205

265

285

85

5

Cumulative soil evaporation, in cubic meters

Simulation day

-0.90

-0.80

-0.70

-0.60

-0.50

-0.40

-0.30

-0.20

-0.10

0

VS2DT results

VSF results

Figure 19. Comparison of VSF and VS2DT model results for cumulative surface soil evaporation. The 

REF1, SEV1, and PND1 Packages were used in the VSF simulation for run 2 (shown here). 
Simulation Examples 37

Case 6 —Three-Dimensional Variably Saturated Hypothetical Case

The 2-D domain of simulation example case 5 is expanded to 3-D here to examine the effects of a transverse surface slope on 

the water table, soil saturation, and flow behavior. The conceptual model is presented in figure 20 and consists of a 530-m by

125-m by 12-m deep rectangular domain representing a small agricultural plot that outlets to an irrigation canal. The subsurface 

domain is 12 m thick with the surface relief consisting of a 3.0-percent slope down toward the canal (x-direction), increasing to 

5.4 percent over the last 50 m to level out at the canal’s edge. The surface plot also has a 1.0-percent slope parallel to the canal 

(y-direction). A 1-m deep steady-flowing irrigation canal is represented by a row of constant-head cells (layers 1 through 7) constraining the water table to the ground surface at the downslope end. The soil media consists of a homogeneous medium-coarse 

soil. The soil parameters, discretization, and boundary conditions are listed in table 7, and the daily totals for precipitation (using 

the REF1 Package) and potential evapotranspiration (using the RZE1 Package) are identical to the values used in simulation example case 5, figure 16 (climatic data for Fort Collins, Colorado, from National Climate Data Center, 2002).

To minimize any flow bias caused by the initial head distribution, several consecutive simulation runs were made (using the 

heads from the previous run as the initial heads) until the water table reached a general equilibrium. The resultant time series of 

the depth to the water table for two locations (355 and 180 m upslope from the irrigation canal, respectively) in grid row 25 is 

plotted with precipitation in figure 21. The results reveal that the water table near the canal (180 m) responds more rapidly and 

with greater magnitude to precipitation than farther upslope. The total number of ponded cells for the simulation is plotted as a 

time series with the actual evapotranspiration computed by the RZE1 Package in figure 22. Large precipitation amounts during 

the late summer and early fall (days 290 to 330) created large ponded areas that disappeared relatively quickly due to the increase 

in evapotranspiration during the summer and fall months. To examine the flow dynamics and soil moisture distribution, a snapshot 

of day 299 (the day after the largest precipitation occurred) is shown in figures 23 and 24. Figure 23 is a cross-section view from 

the front of the model grid of the soil saturation distribution and flow vectors. The flow vectors show steep gradients near the 

ground surface from the large rainfall flux and reveal how the soil-moisture distribution was controlled by the prevailing water 

flow. Figure 24 is the plan view of the soil saturation, reflecting the contours of the water table and the general shape of a variable 

source area that contributed to the irrigation canal.

530 meters

125 meters

Constant-head cells 

Slopey = 

1.0 percent Slopex =

3.0 percent

12 meters

x

y

z

Figure 20. Conceptual model for hypothetical agricultural drainage that feeds a steady irrigation canal. The vertical 

scale is exaggerated 10 times to show surface topography.
38 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Table 7. Model parameters for the three-dimensional hypothetical case.

[x, horizontal distance, in meters; z, elevation, in meters; m, meters; h, hydraulic head, in meters; Δx, cell width, in meters; Δy, cell length, in meters; 

Δz, cell thickness, in meters; Δt, simulation time-step size, in minutes; TMAX, maximum time-step size, in days; TMIN, minimum time-step size, in 

days; Ksat, saturated hydraulic conductivity, in meters per day; ne, porosity, dimensionless; Θr , residual soil saturation, dimensionless; α, soil 

parameter for van Genuchten-Mualem relationships, in meters-1; m/d, meters per day; m-1, meters-1; nv, soil parameter for van Genuchten-Mualem 

relationships, dimensionless]

Parameter Parameter value

Model conditions

Initial conditions Final heads of third simulation run

Boundary conditions Constant-head cells:

x = 530.0 m, 9.03 m < z < 8.03 m,

h = 9.03 m

 Temporally variable recharge and evapotranspiration (RZE1 

Package) with surface ponding (PND1 Package)

Discretization Δx = Δy = 5 m;

Δz = 0.1 to 6 m;

Adaptive time-stepping: 

Initial Δt = 2 minutes

TMAX = 0.1 day, TMIN = 1.0 E-5 day

Soil properties

Medium-coarse soil Ksat =1.0 m/d

ne = 0.381, Θr = 0.3937

van Genuchten-Mualem parameters:

α = 1.6 m-1, nv = 2.7

Simulation day

Depth to water, in meters

Precipitation, in meters

4565

105145

165

185

225

245

305

325

345

365

25

125

205

265

285

85

5

180 meters upslope

Daily precipitation

355 meters upslope

0

0.5

1.0

1.5

2.0

2.5

3.0

3.5 0

0.005

0.010

0.015

0.020

0.025

Figure 21. Simulated depth to water under drainage swale for locations 180 and 355 meters upslope from the 

irrigation canal plotted with daily precipitation. 
Simulation Examples 39

0

100

200

300

-0.0010

-0.0008

-0.0006

-0.0004

-0.0002

0

4565

105145

165

185

225

245

305

325

345

365

25

125

205

265285

85

5

Simulation day

Actual evapotranspiration,

in meters

Total number of ponded cells 

Figure 22. Simulated ponded-surface cell totals plotted with daily actual evapotranspiration computed by 

RZE1 Package. Soil saturation, dimensionless

Constant-headcell

x

z

1.00

0.90

0.80

0.70

0.60

0.50

0.40

Figure 23. Longitudinal cross-section view showing simulated soil saturations and flow vectors along grid row 1 at the start 

of day 299. Note that the flow vectors are scaled relative to magnitude, and the z scale is exaggerated 10 times.
40 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Summary

A new process for simulating three-dimensional (3-D) variably saturated flow (VSF) using Richards’ equation has been added 

to the 3-D modular finite-difference ground-water model MODFLOW. Five new packages are presented here as part of the VSF 

Process, and the Adaptive Time-Stepping (ATS1) Package also is presented for use with both the VSF Process and Ground-Water 

Flow (GWF) Process. Model performance is evaluated by comparison to an analytical solution for one-dimensional (1-D) constant-head infiltration (Dirichlet boundary condition), field experimental data for a 1-D constant-head infiltration, laboratory 

experimental data for two-dimensional (2-D) constant-flux infiltration (Neumann boundary condition), laboratory experimental 

data for 2-D transient drainage through a seepage face, and numerical model results (VS2DT) of a 2-D flow-path simulation using

realistic surface boundary conditions. A hypothetical 3-D example case also is presented to demonstrate the new capability using 

periodic boundary conditions and varied surface topography over a larger spatial scale (0.133 km2

). 

This new MODFLOW process was developed with a primary focus of retaining the flexibility that is at the core of the 

MODFLOW concept. The resulting code retains the modular structure characteristic of MODFLOW and retains compatibility 

with most existing features and packages. The code is easily adapted to function with a commercial pre- and post-processing software package, thus greatly expanding the model’s capability to simulate the flow above and below the water table. Soil saturation, dimensionless

x

y

Constant-head cell

Ponded cell

1.00

0.90

0.80

0.70

0.60

0.50

0.40

Figure 24. Plan view of simulated soil saturations and ponded cells for day 299.

References Cited

Anderman, E.R., and Hill, M.C., 2000, MODFLOW–2000, the 

U.S. Geological Survey modular ground-water model—documentation of the Hydrogeologic-Unit Flow (HUF) Package: 

U.S. Geological Survey Open-File Report 00–342, 89 p.

Baker, D.L., 1995, Darcian weighted interblock conductivity 

means for vertical unsaturated flow: Ground Water, v. 33, 

no. 3, p. 385–390.

Broadbridge, P., and White, I., 1988, Constant rate rainfall infiltration—a versatile nonlinear model, part 1. Analytic solution: Water Resources Research, v. 24, p. 145–154.

Brooks, R.H., and Corey, A.T., 1966, Properties of porous 

media affecting fluid flow: Journal of Irrigation Drainage 

Division, v. 92 (IR2), p. 61–88.

Brutsaert, W.F., 1971, A function iteration technique for solving the Richards equation applied to two dimensional 

infiltration problems: Water Resources Research, v. 7, no. 6, 

p. 1516–1583.

Campbell, G.S., 1977, An introduction to environmental biophysics: New York, Springer Verlag, 159 p.

Celia, M.A., Bouloutas, E.T., and Zarba, R.L., 1990, A general 

mass-conservative numerical solution for the unsaturated 

flow equation: Water Resources Research, v. 26, p. 1483–

1496.

Clement, T.P., 1997, RT3D—a modular computer code for simulating reactive multi-species transport in 3-dimensional 

groundwater systems: Pacific Northwest National Lab–

11720, 59 p.

Clement, T.P., Wise, W.R., and Molz, F.J., 1994, A physically 

based, two-dimensional, finite-difference algorithm for modeling variably-saturated flow: Journal of Hydrology, v. 161, 

p. 71–90.
References Cited 41

Desbarats, A.J., 1995, An interblock conductivity scheme for 

finite difference models of steady unsaturated flow in 

heterogeneous media: Water Resources Research, v. 31, 

no. 11, p. 2883–2889.

Dingman, S.L., 2002, Physical hydrology (2nd ed.): Upper 

Saddle River, New York, Prentice Hall, 646 p.

Downer, C.W., and Ogden, F.L., 2004, Appropriate vertical discretization of Richards’ equation for two-dimensional watershed-scale modeling: Hydrological Processes, v. 18, p. 1–22.

Freeze, R.A.,1978, Mathematical models of hillslope hydrology, in Kirkby, M.J. ed., Hillslope hydrology: New York, 

Wiley, p. 227–293.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., 

2000, MODFLOW–2000, the U.S. Geological Survey modular ground-water model—user guide to modularization concepts and the Ground-Water Flow Process: U.S. Geological 

Survey Open-File Report 00–92, 121 p.

Haverkamp, R., Vauclin, R.M., Touma, J., Wierenga, P.J., and 

Vachaud, G., 1977, A comparison of numerical simulation 

models for one-dimensional infiltration: Journal of Soil Science Society of America, v. 41, p. 285–294.

Healy, R.W., 1990, Simulation of solute transport in variably 

saturated porous media with supplemental information on 

modifications to the U.S. Geological Survey’s computer program VS2D: U.S. Geological Survey Water-Resources 

Investigations Report 90–4025, 125 p.

Hill, M.C., Banta, E.R., Harbaugh, A.W., and Anderman, E.R., 

2000, MODFLOW–2000, the U.S. Geological Survey modular ground-water model—user guide to the Observation, 

Sensitivity, and Parameter-Estimation Processes and three 

post-processing programs: U.S. Geological Survey OpenFile Report 00–184, 210 p.

Huyakorn, P.S., and Pinder, G.F., 1983, Computational methods in subsurface flow: Orlando, Florida, Academic Press 

Inc., 473 p.

Hydrogeologic, 1997, MODFLOW–SURFACT, version 2.1, 

Code documentation and user’s guide: Herndon, Virginia, 

Hydrogeologic Inc., 231 p. 

Jensen, M.E., Burman, R.D., and Allen, R.G., 1990, Evapotranspiration and irrigation water requirements: ASCE Manuals and Reports on Engineering Practice, No. 70, 332 p.

Lappala, E.G., Healy, R.W., and Weeks, E.P., 1987, Documentation of computer program VS2D to solve the equations of 

fluid flow in variably saturated porous media: U.S. Geological Survey Water-Resources Investigations Report 83–4099, 

184 p.

Mehl, S.W., and Hill, M.C., 2002, Development and evaluation 

of a local grid refinement method for block-centered finitedifference groundwater models using shared nodes: 

Advances in Water Resources, v. 25, p. 497–511.

Millington, R.J., and Quirk, J.P., 1961, Permeability of porous 

solids: Transactions of Faraday Society, v. 57, p. 1200–1206.

Mualem, Y., 1976, A new model for predicting the hydraulic 

conductivity of unsaturated porous media: Water Resources 

Research, v. 12, p. 513–522.

National Climate Data Center, 2002, Local climatological data 

(edited): Information available on the Web, accessed September 3, 2002, at http://www.ncdc.noaa.gov/oa/ncdc.html

Niswonger, R.G., Prudic, D.E., and Regan, R.S., 2006, Documentation of the Unsaturated Zone Flow (UZF1) Package for 

modeling percolation between the land surface and the water 

table with MODFLOW–2000: U.S. Geological Survey Techniques and Methods 6A–19, 74 p.

Paniconi, C., and Putti, M., 1994, A comparison of Picard and 

Newton iteration in the numerical solution of multidimensional variably saturated flow problems: Water Resources 

Research, v. 30, no. 12, p. 3357–3374.

Phillip, J.R., 1969, Theory of infiltration: Advances in Hydroscience, v. 5, p. 215–305.

Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Simunek, 

J., and Flerchinger, G.N., 2002, Intercode comparisons for 

simulating water balance of surficial sediments in semiarid 

regions: Water Resources Research, v. 38, no. 12, p. 59-1–

59-16.

Simunek, J.K., Huang, K., and van Genuchten, M., 1995, The 

SWMS_3D code for simulating water flow and solute transport in three-dimensional variably-saturated porous media, 

version 1.0: Riverside, California, U.S. Department of Agriculture, Salinity Laboratory, Research Report No. 139, 155 p.

Thoms, R.B., 2003, Simulating fully coupled overland and variably saturated subsurface flow using MODFLOW: Portland, 

Oregon, Oregon Health and Science University, Department 

of Environmental Science and Engineering, master’s thesis, 

141 p.

VanderKwaak, J.E., 1999, Numerical simulation of flow and 

chemical transport in integrated surface-subsurface hydrologic systems: Waterloo, Ontario, Canada, University of 

Waterloo, Department of Earth Sciences, Ph.D. dissertation, 

217 p.

van Genuchten, M. Th., 1978, Mass transport in saturatedunsaturated media—one-dimensional solutions: Princeton, 

New Jersey, Princeton University, Research Report 78–WR–

11, 102 p.

van Genuchten, M. Th., 1980, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils: Journal of Soil Science Society of America, v. 44, p. 892–898.

van Genuchten, M.Th., Simunek, J., Leij, F.J., and Sejna, 

M.,1997, RETC (RETention Curve)—code for quantifying 

the hydraulic functions of unsaturated soils, version 6.0: Riverside, California, U.S. Salinity Laboratory, 85 p.

Vauclin, M., Khanji, J., and Vachaud, G., 1979, Experimental 

and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem: Water 

Resources Research, v. 15, no. 5, p. 1089–1101.

Vauclin, M., Vachaud, G., and Khanji, J., 1975, Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils, in Vansteenkiste, G.C., ed., Modeling and simulation of water resources systems: Amsterdam, 

North-Holland, p. 299–323.

Warrick, A.W., 1991, Numerical approximations of Darcian 

flow through unsaturated soil: Water Resources Research, 

v. 27, p. 1215–1222.

Warrick, A.W., Biggar, J.W., and Nielsen, D.R., 1971, Simultaneous solute and water transfer for an unsaturated soil: 

Water Resources Research, v. 7, p. 1216–1225.
42 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Zaidel, J., and Russo, D., 1992, Estimation of finite difference 

interblock conductivities for simulation of infiltration into 

initially dry soils: Water Resources Research, v. 28, no. 9, 

p. 2285–2295.

Zheng, Chunmiao, and Wang, P.P., 1999, MT3DMS, A modular three-dimensional multi-species transport model for 

simulation of advection, dispersion and chemical reactions of 

contaminants in groundwater systems—documentation and 

user’s guide: Vicksburg, Mississippi, U.S. Army Engineer 

Research and Development Center, Contract Report 

SERDP–99–1, 202 p.
Appendix 1 43

Appendix 1

Modified Picard Iteration Technique

The governing equation of MODFLOW was changed from a “head-based” form to a “mixed” form to incorporate the 

modified Picard iteration technique (Celia and others, 1990). This technique is used to minimize errors resulting from temporal 

discretization and to maintain mass balance for a range of time-step sizes. This technique was chosen because it is easier to implement and yet performs just as well as the Newton-Raphson iteration under most conditions (Paniconi and Putti, 1994). The modified Picard linearization technique is used in the backward Euler approximation through the Taylor series expansion of θ n+1,m+1

with respect to total head (h):

, (14)

, (15)

where

n = time-step level;

m = iteration level;

Δt = t

n

-t

n-1; 

Hn,m = approximate value for h at nth discrete time level (t = tn) and mth iteration; and

O( ) = infinite sum of partial derivatives of second order and higher.

Neglecting all terms higher than linear in equation 15 and combining equations 14 and 15 with the storage terms of equation 3 

results in the approximation:

, (16)

In the MODFLOW model formulation, the finite-difference form of equation 5 for a grid cell is (Harbaugh and others, 2000, p. 11):

(17)

,

where

CR, CC, CV = hydraulic conductances between node (i,j,k) and an adjacent node [L2/T];

= total head in cell (i,j,k) at time step n [L];

Pi,j,k = sum of head coefficients from source and sink terms [L2/T];

Qi,j,k = sum of constants from source and sink terms [L3/T];

SSi,j,k = specific storage [L-1];

∂θ

∂t

------ C∂h

∂t

----- θn m, θn – 1 –

Δt = ≅ ------------------------------

θn m, θn m, – 1 dθ

dh

------

n m, – 1

Hn m, Hn m, – 1 ( ) – O ∂

2 = + + ( )

∂

2

C∂h

∂t

----- ΘSs

∂h

∂t

----- θn m, – 1 Cn m, – 1 Hn m, Hn m, – 1 ( ) – θn – 1

+ –

Δt ------------------------------------------------------------------------------------------------------ Θn – 1/2Ss

Hn m, Hn – 1 –

Δt

+ ≅ + --------------------------------

CRi j , – 1/2, k Hi j , – 1, k

n Hijk , ,n ( ) – CRi j , + 1/2, k Hi j , + 1, kn Hijk , ,n + ( ) –

CCi–1/2, ,j k Hi – 1, ,j k

n Hijk , ,n ( ) – CCi + 1/2, ,j k Hi + 1, ,j kn Hijk , ,n + + ( ) –

CVijk , , – 1/2 Hijk , , – 1

n Hijk , ,n ( ) – CVijk , , + 1/2 Hijk , , + 1n Hijk , ,n + + ( ) +

Pijk , , Hijk , ,

n + + Qijk , , SSijk , , DELRj × DELCi × THICKijk , , ( )

Hijk , ,

n Hijk , ,n – 1 –

Δt = -----------------------------

, ,

n Hi jk
44 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

DELRj = cell width of column j in all rows [L];

DELCi = cell width of row i in all columns [L];

THICKi,j,k = vertical thickness of cell (i,j,k) [L]; and

Δt = change in time between time step n and time step n-1 [T].

This is simplified further for solution by computer into:,

, (18)

where

; and

.

To implement the modified Picard iteration technique into MODFLOW, the HCOF and RHS terms in equation 18 are changed to 

incorporate the unsaturated storage terms:

; and (19)

. (20)

CVijk , , – 1/2Hijk , , – 1

n CCi – 1/2,j k, Hi – 1, ,j kn CRi j , – 1/2,kHi j , – 1, kn + +

–CVijk , , – 1/2 CCi – 1/2, ,j k – CRi j , – 1/2, k – CVijk , , + 1/2, j k, – CCi + 1/2, j k, – CRi j , + 1/2, k – + HCOFijk , , ( )Hijk , ,

n +

CVijk , , + 1/2Hijk , , + 1

n CCi + 1/2, j k, Hi + 1, ,j kn CRi j , + 1/2, kHi j , + 1, kn =RHS + + ijk , ,

HCOFijk , , Pijk , ,

SSijk , , DELCi × DELRj × THICKk ( )

Δt = – ----------------------------------------------------------------------------------------

RHSijk , , Qijk , , – SSijk , , DELCi × DELRj × THICKk ( )

Δt ----------------------------------------------------------------------------------------Hijk , ,

n – 1 = –

HCOF Pijk , ,

Cn m, – 1 Θn – 1/2SS + ijk , ,

Δt ------------------------------------------------------- DELCi × DELRj × THICKk = – ( )

RHS Qijk , , – θn m, – 1 θn – 1 –

t

n tn – 1 –

------------------------------------ Cn m, – 1 Θn – 1/2SS + ijk , ,

Δt -------------------------------------------------------Hijk , ,

n – 1 – DELCi × DELRj × THICKk = + ( )
Appendix 2 45

Appendix 2 

REF1 Package Example Input and Output

Example input and resultant output files for the REF1 Package are shown for simulation example case 3 (2-D infiltration) in 

figures 25 and 26, respectively. The model budget written by the Basic (BAS6) Package is shown for the final time step of the

simulation in figure 27. For the input file (fig. 25), the output file number is set as 39 (IREFCB=39); one simulation time is 

requested for output of soil saturations (IREFOC=1, ISATOC=1); there is only one soil type (SLNUM=1), and 1,000 specified 

pressure values are requested for the internal table (NTAB=1,000); the intercell averaging method is upstream weighted, and the

initial saturations are determined from initial head distribution in BASIC input file (IAVG=3, ISATFL, IHYDST); output is 

requested in the first stress period for the last time step (REFOC(1)=1, REFOC(2)=-500); and the pore-water pressure head interval 

is from PTABA = -0.001 to PTABB = -200 m. The van Genuchten-Mualem soil parameters are singular (ISC=1) and are defined 

as α = 3.3 m-1 and n = 4.1. Listed next are the residual soil saturation value (Θr=0.033) and the porosity (n=0.3). The soil type is 

the same for all 40 layers (STYPE=1). The output file (fig. 26) contains the soil saturation distribution for each node location at 

the final time step in the first stress period, total elapsed time of 0.333 day.
46 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

# IREFCB

39

# IREFOC ISATOC

1 1

# SLNUM NTAB

1 1000

# IAVG ISATFL IHYDST

3 0 0

# REFOC

1 -500

# PTABA PTABB

-0.001 -200.

# ISC

1

# ALPHA VGN

3.3 4.1

# RSAT EFFP

0.033 0.3

CONSTANT 1 STYPE Layer 1

CONSTANT 1 STYPE Layer 2

CONSTANT 1 STYPE Layer 3

CONSTANT 1 STYPE Layer 4

CONSTANT 1 STYPE Layer 5

CONSTANT 1 STYPE Layer 6

CONSTANT 1 STYPE Layer 7

CONSTANT 1 STYPE Layer 8

CONSTANT 1 STYPE Layer 9

CONSTANT 1 STYPE Layer 10

CONSTANT 1 STYPE Layer 11

CONSTANT 1 STYPE Layer 12

CONSTANT 1 STYPE Layer 13

CONSTANT 1 STYPE Layer 14

CONSTANT 1 STYPE Layer 15

CONSTANT 1 STYPE Layer 16

CONSTANT 1 STYPE Layer 17

CONSTANT 1 STYPE Layer 18

CONSTANT 1 STYPE Layer 19

CONSTANT 1 STYPE Layer 20

CONSTANT 1 STYPE Layer 21

CONSTANT 1 STYPE Layer 22

CONSTANT 1 STYPE Layer 23

CONSTANT 1 STYPE Layer 24

CONSTANT 1 STYPE Layer 25

CONSTANT 1 STYPE Layer 26

CONSTANT 1 STYPE Layer 27

CONSTANT 1 STYPE Layer 28

CONSTANT 1 STYPE Layer 29

CONSTANT 1 STYPE Layer 30

CONSTANT 1 STYPE Layer 31

CONSTANT 1 STYPE Layer 32

CONSTANT 1 STYPE Layer 33

CONSTANT 1 STYPE Layer 34

CONSTANT 1 STYPE Layer 35

CONSTANT 1 STYPE Layer 36

CONSTANT 1 STYPE Layer 37

CONSTANT 1 STYPE Layer 38

CONSTANT 1 STYPE Layer 39

CONSTANT 1 STYPE Layer 40

Figure 25. Example input for

Richards’ Equation Flow (REF1) Package. 
Appe

ndix 2 47

500 1 3.333337E-01 3.333337E-01 SATURATION 30 1 1 ( 30F8.4)0.8322 0.8304 0.8254 0.8123 0.7712 0.5721 0.4236 0.3187 0.2437 0.1855 0.1313 0.0763 0.0446 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.0430 0.04300.8315 0.8289 0.8219 0.8044 0.7559 0.6115 0.4757 0.3667 0.2831 0.2163 0.1559 0.0953 0.0502 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.0442 0.04420.8311 0.8277 0.8187 0.7978 0.7469 0.6318 0.5100 0.4032 0.3158 0.2440 0.1795 0.1154 0.0605 0.0459 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.0457 0.04570.8313 0.8270 0.8163 0.7928 0.7420 0.6449 0.5351 0.4324 0.3441 0.2691 0.2018 0.1358 0.0740 0.0482 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.0474 0.04740.8326 0.8276 0.8153 0.7899 0.7401 0.6550 0.5551 0.4570 0.3691 0.2924 0.2229 0.1549 0.0890 0.0518 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.0495 0.04950.8360 0.8302 0.8165 0.7898 0.7414 0.6647 0.5728 0.4793 0.3925 0.3148 0.2435 0.1737 0.1046 0.0573 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.0520 0.05200.8427 0.8361 0.8212 0.7935 0.7465 0.6759 0.5904 0.5011 0.4157 0.3374 0.2648 0.1934 0.1216 0.0652 0.0551 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.0550 0.05500.8542 0.8471 0.8311 0.8028 0.7570 0.6911 0.6107 0.5249 0.4410 0.3624 0.2886 0.2160 0.1424 0.0791 0.0594 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.0587 0.05870.8724 0.8648 0.8483 0.8198 0.7755 0.7132 0.6371 0.5542 0.4714 0.3924 0.3178 0.2450 0.1725 0.1094 0.0744 0.0660 0.0640 0.0634 0.0633 0.0632 0.0632 0.0632 0.0632 0.0632 0.0632 0.0632 0.0632 0.0632 0.0632 0.06320.8977 0.8902 0.8740 0.8465 0.8044 0.7459 0.6736 0.5934 0.5114 0.4320 0.3566 0.2842 0.2148 0.1562 0.1167 0.0949 0.0821 0.0748 0.0713 0.0698 0.0692 0.0690 0.0689 0.0688 0.0688 0.0688 0.0688 0.0688 0.0688 0.06880.9280 0.9214 0.9071 0.8826 0.8448 0.7916 0.7244 0.6475 0.5668 0.4870 0.4107 0.3387 0.2720 0.2155 0.1745 0.1464 0.1257 0.1094 0.0968 0.0877 0.0821 0.0789 0.0773 0.0766 0.0762 0.0760 0.0760 0.0759 0.0759 0.07590.9580 0.9531 0.9423 0.9233 0.8930 0.8487 0.7899 0.7193 0.6420 0.5630 0.4861 0.4138 0.3475 0.2895 0.2451 0.2114 0.1845 0.1621 0.1432 0.1273 0.1142 0.1040 0.0965 0.0916 0.0887 0.0870 0.0861 0.0856 0.0853 0.08520.9811 0.9783 0.9718 0.9599 0.9397 0.9080 0.8629 0.8044 0.7358 0.6614 0.5861 0.5132 0.4454 0.3834 0.3326 0.2916 0.2577 0.2287 0.2035 0.1817 0.1627 0.1465 0.1328 0.1219 0.1134 0.1073 0.1032 0.1007 0.0993 0.09860.9942 0.9930 0.9902 0.9846 0.9742 0.9564 0.9281 0.8876 0.8351 0.7729 0.7049 0.6351 0.5667 0.5019 0.4429 0.3928 0.3497 0.3124 0.2796 0.2507 0.2251 0.2027 0.1832 0.1664 0.1523 0.1408 0.1320 0.1256 0.1215 0.11950.9991 0.9988 0.9981 0.9963 0.9927 0.9854 0.9721 0.9505 0.9186 0.8760 0.8241 0.7652 0.7027 0.6392 0.5769 0.5186 0.4660 0.4188 0.3766 0.3387 0.3049 0.2748 0.2482 0.2248 0.2045 0.1872 0.1731 0.1621 0.1545 0.15051.0000 1.0000 0.9999 0.9997 0.9990 0.9972 0.9931 0.9850 0.9709 0.9489 0.9182 0.8786 0.8314 0.7783 0.7215 0.6632 0.6055 0.5504 0.4989 0.4512 0.4076 0.3680 0.3324 0.3007 0.2727 0.2484 0.2280 0.2116 0.1998 0.19331.0000 1.0000 1.0000 1.0000 1.0000 0.9999 0.9993 0.9976 0.9937 0.9860 0.9731 0.9536 0.9267 0.8922 0.8508 0.8035 0.7519 0.6976 0.6425 0.5881 0.5360 0.4869 0.4412 0.3996 0.3622 0.3291 0.3007 0.2775 0.2605 0.25101.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999 0.9995 0.9982 0.9949 0.9887 0.9781 0.9622 0.9401 0.9111 0.8754 0.8335 0.7864 0.7355 0.6823 0.6284 0.5754 0.5247 0.4776 0.4346 0.3968 0.3653 0.3417 0.32851.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9997 0.9988 0.9966 0.9920 0.9842 0.9720 0.9544 0.9306 0.9003 0.8635 0.8208 0.7731 0.7220 0.6693 0.6166 0.5660 0.5195 0.4793 0.4484 0.43081.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999 0.9994 0.9981 0.9951 0.9896 0.9807 0.9671 0.9481 0.9228 0.8909 0.8526 0.8088 0.7607 0.7106 0.6612 0.6160 0.5797 0.55831.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9998 0.9991 0.9974 0.9938 0.9875 0.9774 0.9625 0.9418 0.9148 0.8814 0.8425 0.8001 0.7578 0.7214 0.69891.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9997 0.9988 0.9966 0.9924 0.9851 0.9737 0.9573 0.9353 0.9081 0.8777 0.8487 0.82941.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999 0.9995 0.9983 0.9956 0.9906 0.9823 0.9702 0.9544 0.9372 0.92451.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9998 0.9992 0.9977 0.9945 0.9892 0.9821 0.97591.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9997 0.9989 0.9973 0.99521.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999 0.99961.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

Figure 26. Example output for Richards’ Equation Flow (REF1) Package.
48 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

The model budget printed to the ‘.OUT’ file by the BAS6 Package is shown in figure 27. The REF1 Package adds to the model 

budget entries for cumulative soil-moisture change and unsaturated storage, which are shown to be substantial components of the

total flow at the end of the simulation for simulation example case 3.

SPF1 Package Example Input and Output

Example input and resultant output files for the SPF1 Package are shown for simulation example case 4 (2-D drainage) in 

figures 28 and 29, respectively. Within the model budget, seepage-face flow is included in MODFLOW’s computation of constant-head flow. For the input file (fig. 28), the output file unit number is set to 83 (ISPFCB=83); output is requested for five simulation times (ISPFOC=5); output is written to the ‘.SPO’ file for time steps 1, 20, 100, 400, and 1,000 in the first stress period 

(SPFOC); and the seepage-face location is defined as the rightmost column of cells between layers 1 and 25 (SEEP Flag array).

The output file (.SPO) generated by the above input file is shown in figure 29. It shows that the active seepage face initially 

conducts large fluxes out of the domain and occurs between layer 16 and layer 25. As the simulation proceeds, the fluxes gradually 

decrease, and the seepage face shrinks and eventually disappears by the final time step.

VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP500 IN STRESS PERIOD 1 

 ---------------------------------------------------------------------------------------------------------------------------------------

 CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T 

 ---------------------------------- ----------------------------------------

IN: IN:

--- ---

 STORAGE = 0.0000 STORAGE = 0.0000

CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000

 RECHARGE = 5.9167 RECHARGE = 17.7500

 SOIL MOISTURE = 2.3804E-02 SOIL MOISTURE = 0.0000

 UNSAT STORAGE = 5.9567E-08 UNSAT STORAGE = 6.4528E-12

 TOTAL IN = 5.9405 TOTAL IN = 17.7500 

 OUT: OUT:

 ---- ----

 STORAGE = 0.0000 STORAGE = 0.0000 

CONSTANT HEAD = 2.6329 CONSTANT HEAD = 15.8951 

 RECHARGE = 0.0000 RECHARGE = 0.0000 

 SOIL MOISTURE = 3.3069 SOIL MOISTURE = 1.8537 

 UNSAT STORAGE = 1.0784E-05 UNSAT STORAGE = 7.3730E-06

 TOTAL OUT = 5.9399 TOTAL OUT = 17.7488 

 IN - OUT = 6.0892E-04 IN - OUT = 1.2169E-03

PERCENT DISCREPANCY = 0.01 PERCENT DISCREPANCY = 0.01 

Figure 27. Example model budget for Richards’ Equation Flow (REF1) Package.
Appendix 2 49

# ISPF1CB

83

# ISPF1OC

5

# SPF1OC

1 -1

1 -20

1 -100

1 -400

1 -1000

INTERNAL 1 (free) 0 SEEP Flag Array Layer 1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 2

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 3

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 4

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 5

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 6

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 7

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 8

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 9

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 10

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 11

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 12

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 13

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 14

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 15

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 16

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 17

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 18

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 19

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 20

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 21

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 22

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 23

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 24

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

INTERNAL 1 (free) 0 SEEP Flag Array Layer 25

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

CONSTANT 0 SEEP Flag Array Layer 26

CONSTANT 0 SEEP Flag Array Layer 27

CONSTANT 0 SEEP Flag Array Layer 28

CONSTANT 0 SEEP Flag Array Layer 29

CONSTANT 0 SEEP Flag Array Layer 30

CONSTANT 0 SEEP Flag Array Layer 31

CONSTANT 0 SEEP Flag Array Layer 32

CONSTANT 0 SEEP Flag Array Layer 33

CONSTANT 0 SEEP Flag Array Layer 34

CONSTANT 0 SEEP Flag Array Layer 35

CONSTANT 0 SEEP Flag Array Layer 36

CONSTANT 0 SEEP Flag Array Layer 37

CONSTANT 0 SEEP Flag Array Layer 38

CONSTANT 0 SEEP Flag Array Layer 39

CONSTANT 0 SEEP Flag Array Layer 40

Figure 28. Example input for Seepage Face (SPF1) Package.
50 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

PND1 Package Example Input and Output

Example input and resultant output files for the Surface Ponding (PND1) Package are shown for simulation example case 5 

(2-D flow-path comparison) in figures 30 and 31, respectively. Within the model budget, surface ponding flow is included in 

MODFLOW’s computation of constant-head flow. For the input file (fig. 30), the output file unit number is set to 38 (IPNDCB=

38); output is requested for six simulation times (IPNDOC=6); output is written to the ‘.PDO’ file for the first time step of stress 

period 297, the first and last time steps of stress periods 298 and 299, and the first time step of stress period 300 (PNDOC); and 

the surface ponding depth is set to 0.01 m (POND). The output file (.PDO) generated by the above input file is shown in figure 31. 

It shows that surface ponding develops during stress period 298 and recedes during stress period 299.

SEEPAGE FACE PERIOD 1 STEP 1 SIMULATION TIME 0.500E-02 # OF SEEPAGE CELLS: 10

LAYER 16 ROW 1 COL 30 FLOW -527.332

LAYER 17 ROW 1 COL 30 FLOW -40.0232

LAYER 18 ROW 1 COL 30 FLOW -79.2796

LAYER 19 ROW 1 COL 30 FLOW -115.791

LAYER 20 ROW 1 COL 30 FLOW -151.955

LAYER 21 ROW 1 COL 30 FLOW -189.760

LAYER 22 ROW 1 COL 30 FLOW -231.333

LAYER 23 ROW 1 COL 30 FLOW -279.488

LAYER 24 ROW 1 COL 30 FLOW -338.557

LAYER 25 ROW 1 COL 30 FLOW -415.878

SEEPAGE FACE PERIOD 1 STEP 20 SIMULATION TIME 0.100 # OF SEEPAGE CELLS: 6

LAYER 20 ROW 1 COL 30 FLOW -409.895

LAYER 21 ROW 1 COL 30 FLOW -75.2486

LAYER 22 ROW 1 COL 30 FLOW -138.303

LAYER 23 ROW 1 COL 30 FLOW -201.392

LAYER 24 ROW 1 COL 30 FLOW -271.296

LAYER 25 ROW 1 COL 30 FLOW -356.758

SEEPAGE FACE PERIOD 1 STEP 100 SIMULATION TIME 0.500 # OF SEEPAGE CELLS: 4

LAYER 22 ROW 1 COL 30 FLOW -22.7561

LAYER 23 ROW 1 COL 30 FLOW -42.1152

92

LAYER 24 ROW 1 COL 30 FLOW -141.372

LAYER 25 ROW 1 COL 30 FLOW -246.855

SEEPAGE FACE PERIOD 1 STEP 400 SIMULATION TIME 2.00 # OF SEEPAGE CELLS: 1

LAYER 25 ROW 1 COL 30 FLOW -195.397

SEEPAGE FACE PERIOD 1 STEP1000 SIMULATION TIME 5.00 # OF SEEPAGE CELLS: 0

Figure 29. Example output for Seepage Face (SPF1) Package.

# IPNDCB

38

# IPNDOC

6

# PNDOC

297 -1

298 -1

298 -24

299 -1

299 -24

300 -1

CONSTANT 0.01 Surface Ponding Depth

Figure 30. Example input for Surface Ponding (PND1) Package.
Appendix 2 51

SEV1 Package Example Input and Output

Example input and resultant output files for the Surface Evaporation (SEV1) Package are shown for run 2 of simulation example case 5 (2-D flow-path comparison) in figures 32 and 33, respectively. For brevity, the input for only the first 3 of 365 stress 

periods is shown here. The model budget written by the Basic (BAS6) Package is shown for the first time step of the simulation 

in figure 34. 

In the input file (fig. 32), the output file unit number is set to 37 (ISEVCB=37); output is requested for three simulation times 

at two cell locations (ISEVOC=3, ISEVCL=2); output is written to the ‘.SEV’ file for last time steps in the first and second stress 

periods and first time step of the fourth stress period (SEVOC); SEV1 output is recorded for all layers within the root zone for the 

horizontal cell locations of column 80, row 1 and column 105, row 1 (SEVCL); all parameters are read in for the first stress period 

(INPEV=INHA=INSRES=INSEL=1); the potential soil evaporation rate for the first stress period is 0.003494 m/d 

(PEV=0.003494); the atmospheric pressure is -0.001 m (HA=-0.001); the surface resistance is 2.0E-4 m-1 (RES=2.0E-4); and 

evaporation is extracted from the second model layer downward (SEL=2).

The resultant output file is shown in figure 33. The evaporation flow rates for cell locations (80,1) and (105,1) are listed for 

the last time steps of the first and second stress periods and for the first time step of the fourth stress period.

The model budget listing for the first time step of the fourth stress period is listed in figure 34. The total soil evaporation rate 

from the second model layer is 0.003178 m3/d, and the cumulative volume up to that time step is 0.010354 m3.

PONDED FLOW PERIOD 297 STEP 1 SIMULATION TIME 296.042 # OF PONDED CELLS: 2

LAYER 1 ROW 1 COL 104 FLOW -0.309360E-02

LAYER 1 ROW 1 COL 105 FLOW -0.396652E-01

PONDED FLOW PERIOD 298 STEP 1 SIMULATION TIME 297.042 # OF PONDED CELLS: 3

LAYER 1 ROW 1 COL 103 FLOW 0.188034E-01

LAYER 1 ROW 1 COL 104 FLOW -0.129519E-01

LAYER 1 ROW 1 COL 105 FLOW -0.425892E-01

PONDED FLOW PERIOD 298 STEP 24 SIMULATION TIME 298.000 # OF PONDED CELLS: 10

LAYER 1 ROW 1 COL 96 FLOW 0.158992E-01

LAYER 1 ROW 1 COL 97 FLOW 0.123461E-01

LAYER 1 ROW 1 COL 98 FLOW 0.961116E-02

LAYER 1 ROW 1 COL 99 FLOW 0.776160E-02

LAYER 1 ROW 1 COL 100 FLOW 0.588828E-02

LAYER 1 ROW 1 COL 101 FLOW 0.319300E-02

LAYER 1 ROW 1 COL 102 FLOW -0.144910E-02

LAYER 1 ROW 1 COL 103 FLOW -0.993651E-02

LAYER 1 ROW 1 COL 104 FLOW -0.255280E-01

LAYER 1 ROW 1 COL 105 FLOW -0.491849E-01

PONDED FLOW PERIOD 299 STEP 1 SIMULATION TIME 298.042 # OF PONDED CELLS: 2

LAYER 1 ROW 1 COL 104 FLOW -0.171942E-01

LAYER 1 ROW 1 COL 105 FLOW -0.444869E-01

PONDED FLOW PERIOD 299 STEP 24 SIMULATION TIME 299.000 # OF PONDED CELLS: 2

LAYER 1 ROW 1 COL 104 FLOW 0.130985E-02

LAYER 1 ROW 1 COL 105 FLOW -0.378838E-01

PONDED FLOW PERIOD 300 STEP 1 SIMULATION TIME 299.042 # OF PONDED CELLS: 1

LAYER 1 ROW 1 COL 105 FLOW -0.347702E-01

Figure 31. Example output for Surface Ponding (PND1) Package.
52 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

# ISEVCB

37

# ISEVOC ISEVCL

3 2

# SEVOC (STRESS PERIOD, TIMESTEP)

1 -24

2 -24

4 -1

# SEVCL (COLUMN, ROW)

80 1

105 1

# INPEV INHA INSRES INSEL

1 1 1 1

CONSTANT 0.003494 PEV

-0.001 HA

CONSTANT 2.0E-4 SRES

CONSTANT 2 SEL

# INPEV INHA INSRES INSEL

1 1 1 1

CONSTANT 0.002487

-0.001

CONSTANT 2.0E-4

CONSTANT 2

# INPEV INHA INSRES INSEL

1 1 1 1

CONSTANT 0.002381

-0.001

CONSTANT 2.0E-4

CONSTANT 2

Figure 32. Example input for Surface Evaporation (SEV1) Package.

SOIL EVAPORATION PERIOD 1 STEP 24 ELAPSED TIME 1.000000 TOTAL SEV1 FLOW FOR THIS TIME STEP = 0.311571E-02

 LAYER 2 ROW 1 COL 80 FLOW 0.00000

 LAYER 2 ROW 1 COL 105 FLOW -0.259478E-03

SOIL EVAPORATION PERIOD 2 STEP 24 ELAPSED TIME 2.000000 TOTAL SEV1 FLOW FOR THIS TIME STEP = 0.338443E-02

 LAYER 2 ROW 1 COL 80 FLOW -0.110939E-03

 LAYER 2 ROW 1 COL 105 FLOW -0.112732E-03

SOIL EVAPORATION PERIOD 4 STEP 1 ELAPSED TIME 3.041667 TOTAL SEV1 FLOW FOR THIS TIME STEP = 0.317823E-02

 LAYER 2 ROW 1 COL 80 FLOW 0.00000

 LAYER 2 ROW 1 COL 105 FLOW -0.259179E-03

Figure 33. Example output for Surface Evaporation (SEV1) Package.
Appendix 2 53

RZE1 Package Example Input and Output

Example input and resultant output files for the RZE1 Package are shown for simulation example case 5 (2-D flow-path comparison) in figures 35 and 36, respectively. For brevity, the input for only the first 3 of 365 stress periods is shown here. The model 

budget written by the Basic (BAS6) Package is shown for the first time step of the simulation in figure 37. 

 In the input file (fig. 35), the output file unit number is set to 37 (IRZECB=37); output is requested for two simulation times 

at two cell locations (IRZEOC=2, IRZECL=2); output is written to the ‘.RZO’ file for last time step in the first stress period and 

the first time step of the third stress period (RZEOC); RZE1 output is recorded for all layers within the root zone for the horizontal 

cell locations of column 1, row 1 and column 105, row 1 (RZECL); all parameters are read in for the first stress period (INPET=

INRTDPTH=INRTBOT=INRTTOP=INHROOT=INRZEL=1); the potential evapotranspiration rate for the first stress period is 

0.003494 m/d (PET=0.003494); the root depth is 0.6 m (RTDPTH=0.6); the root activity at the bottom depth is 

2.0E-5 m-2 (RTBOT=2.0E-5); the root activity at the top depth is 9.0E-5 m-2 (RTTOP=9.0E-5); the pressure in the roots is -80.0 m 

(HROOT=-80.0); and evapotranspiration is extracted from the second model layer downward (RZL=2).

The resultant output file is shown in figure 36. The evapotranspiration flow rates for cell locations (1,1) and (105,1) are listed 

for the last time step of the first stress period and for the first time step of the third stress period.

The model budget listing for the first time step of the third stress period is listed in figure 37. The total evapotranspiration 

flow rate from the root zone is 0.3817 m3/d, and the cumulative volume up to that time step is 0.7606 m3.

 VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 1 IN STRESS PERIOD 4 

 --------------------------------------------------------------------------------------------------------------------------------

 CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T

 ---------------------------------- --------------------------------------

 IN: IN: 

 --- ---

STORAGE = 0.0000 STORAGE = 0.0000

CONSTANT HEAD = 3.6747E-02 CONSTANT HEAD = 0.0000

SOIL EVAPORATION = 0.0000 SOIL EVAPORATION = 0.0000

RECHARGE = 1.8926 RECHARGE = 0.0000

SOIL MOISTURE = 2.2538 SOIL MOISTURE = 0.7906

UNSAT STORAGE = 5.3151E-03 UNSAT STORAGE = 9.5146E-04

TOTAL IN = 4.1885 TOTAL IN = 0.7916

OUT: OUT:

 ---- ----

STORAGE = 0.0000 STORAGE = 0.0000

CONSTANT HEAD = 1.0416 CONSTANT HEAD = 0.3294

SOIL EVAPORATION = 1.0354E-02 SOIL EVAPORATION = 3.1782E-03

RECHARGE = 0.0000 RECHARGE = 0.0000

SOIL MOISTURE = 3.1330 SOIL MOISTURE = 0.4587

UNSAT STORAGE = 4.3073E-03 UNSAT STORAGE = 1.3255E-04

TOTAL OUT = 4.1892 TOTAL OUT = 0.7914

IN - OUT = -7.6151E-04 IN - OUT = 2.2542E-04

PERCENT DISCREPANCY = -0.02 PERCENT DISCREPANCY = 0.03

Figure 34. Example model budget for Surface Evaporation (SEV1) Package.
54 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

# IRZECB

37

# IRZEOC IRZE1CL

2 2

# RZEOC (STRESS PERIOD, TIMESTEP)

1 -24

3 -1

# RZECL (COLUMN, ROW)

1 1

105 1

# INPET INRTDPTH INRTBOT INRTTOP INHROOT INRZL

1 1 1 1 1 1

CONSTANT 0.003494 PET

CONSTANT 0.6 RTDPTH

CONSTANT 2.0E-5 RTBOT

CONSTANT 9.0E-5 RTTOP

CONSTANT -80.0 HROOT

CONSTANT 2 RZL

# INPET INRTDPTH INRTBOT INRTTOP INHROOT INRZL

1 1 1 1 1 1

CONSTANT 0.002487

CONSTANT 0.6

CONSTANT 2.0E-5

CONSTANT 9.0E-5

CONSTANT -80.0

CONSTANT 2

# INPET INRTDPTH INRTBOT INRTTOP INHROOT INRZL

1 1 1 1 1 1

CONSTANT 0.002381

CONSTANT 0.6

CONSTANT 2.0E-5

CONSTANT 9.0E-5

CONSTANT -80.0

CONSTANT 2

Figure 35. Example input for Root Zone Evapotranspiration (RZE1) Package.

STRESS PERIOD 1 STEP 24 ELAPSED TIME 1.000000 TOTAL RZE FLOW FOR THIS TIME STEP = 0.283808 

LAYER 2 ROW 1 COL 1 FLOW 0.389412E-05

LAYER 3 ROW 1 COL 1 FLOW 0.316735E-05

LAYER 4 ROW 1 COL 1 FLOW 0.260399E-05

LAYER 5 ROW 1 COL 1 FLOW 0.588196E-06

LAYER 2 ROW 1 COL 105 FLOW 0.267203E-02

LAYER 3 ROW 1 COL 105 FLOW 0.220719E-02

LAYER 4 ROW 1 COL 105 FLOW 0.226088E-02

LAYER 5 ROW 1 COL 105 FLOW 0.736164E-03

STRESS PERIOD 3 STEP 1 ELAPSED TIME 2.041667 TOTAL RZE FLOW FOR THIS TIME STEP = 0.303337

LAYER 2 ROW 1 COL 1 FLOW 0.665321E-05

LAYER 3 ROW 1 COL 1 FLOW 0.446539E-05

LAYER 4 ROW 1 COL 1 FLOW 0.340623E-05

LAYER 5 ROW 1 COL 1 FLOW 0.804811E-06

LAYER 2 ROW 1 COL 105 FLOW 0.267204E-02

LAYER 3 ROW 1 COL 105 FLOW 0.220720E-02

LAYER 4 ROW 1 COL 105 FLOW 0.226090E-02

LAYER 5 ROW 1 COL 105 FLOW 0.736175E-03

Figure 36. Example output for Root Zone Evapotranspiration (RZE1) Package.
Appendix 2 55

ATS1 Package Example Input and Output

Example input for the ATS1 Package is shown in figure 38. The parameters for example simulation case 5 are shown; an 

output file is not requested (ITIM=0), the maximum number of steps in each stress period is 8,640 (MXSTP=8640), the maximum 

time-step size is 0.5 day (TMX=0.5), the minimum time-step size is 0.01 day (TMN=0.01), the time-step multiplier is 1.5 

(TSM=1.5), and the time-step reducer is 2.0 (TSD=2.0).

If an output file is request by the user, the ATS1 Package creates a binary file with the ‘.TIM’ extension. At the end of each

stress period, two entries are written into the output file; a header and an array of REAL values. The header contains the stress 

period number (KPER, INTEGER), the stress period length (PERLEN, REAL), and the total number of time steps (NSTP, 

INTEGER) used by ATS1 to reach the end of the stress period. The 1-D array of REAL values lists each time-step size in order 

from the first to the last time step (1 to NSTP).

 VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 1 IN STRESS PERIOD 3 

-------------------------------------------------------------------------------------------------------------------------------- 

CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T 

----------------------------------- --------------------------------------- 

IN: IN:

--- --- 

STORAGE = 0.0000 STORAGE = 0.0000 

CONSTANT HEAD = 0.1162 CONSTANT HEAD = 0.0000 

ROOT ZONE ET = 0.0000 ROOT ZONE ET = 0.0000 

RECHARGE = 2.1263 RECHARGE = 0.0000 

SOIL MOISTURE = 1.5536 SOIL MOISTURE = 1.5399 

UNSAT STORAGE = 7.0756E-03 UNSAT STORAGE = 4.3591E-02 

TOTAL IN = 3.8032 TOTAL IN = 1.5835 

OUT: OUT:

----- ----- 

STORAGE = 0.0000 STORAGE = 0.0000 

CONSTANT HEAD = 0.5961 CONSTANT HEAD = 0.2857 

ROOT ZONE ET = 0.6024 ROOT ZONE ET = 0.3033 

RECHARGE = 0.0000 RECHARGE = 0.0000 

SOIL MOISTURE = 2.5996 SOIL MOISTURE = 0.9941 

UNSAT STORAGE = 5.0985E-03 UNSAT STORAGE = 3.0185E-04 

TOTAL OUT = 3.8032 TOTAL OUT = 1.5835 

IN - OUT = 7.8678E-06 IN - OUT = -1.9789E-05 

PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00

Figure 37. Example model budget for the Root Zone Evapotranspiration (RZE1) Package. 

# Adaptive time stepping parameters for Case 5 2D Flowpath comparison

#ITIM, MXSTP

0 8640

#TMX, TMN, TSM, TSD 

0.5 0.01 1.5 2.0 

Figure 38. Example input for Adaptive Time-Stepping (ATS1) Package. 
56 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

Appendix 3

Custom Soil Characteristic Functions

This section is intended for advanced users with knowledge of FORTRAN who would like to define soil characteristic functions of a different form than the default van Genuchten-Mualem relationships (equations 6–8). The FORTRAN FUNCTIONS 

defining the soil pore pressure, relative permeability, soil saturation, and specific moisture capacity relationships, are listed in figure 39. This code can be found in the file VSF1_REF1_SOILCHAR.F. Each FORTRAN FUNCTION represents a unique relation; SVSF1REF1PS calculates soil pore-water pressure for a given effective saturation, SVSF1REF1KP calculates relative permeability for a given soil pore-water pressure, SVSF1REF1SP calculates soil saturation for a given pore-water pressure, and 

SVSF1REF1CP calculates specific moisture capacity for a given pore-water pressure. 

The integer flag, ISC (read from the REF1 input file), is used to provide flexibility in defining the soil characteristic functions. 

The value of the integer flag, ISC, determines the form of the relationship used within each FORTRAN FUNCTION. For example, 

simulation example case 1 uses a logarithmic relationship between soil pore pressure and relative permeability for the clay soil 

(Haverkamp and others, 1977) with a dual set of parameters (that is, a unique set for kr and Θ, respectively); setting the ISC value 

in the REF1 input file to -3 will access this custom soil characteristic function. Simulation example case 4 also uses a custom 

relationship with a dual set of soil parameters (Vauclin and others, 1975), which is accessed by setting ISC = -2. ISC values greater 

than zero use the default van Genuchten-Mualem functions, and ISC values less than zero access custom functions. Dual parameter 

sets can be used by setting the absolute value of ISC to greater than 1. For simplicity, an error is returned if the absolute value of 

ISC is less than 1 or greater than 3.

To insert a custom set of soil characteristic functions, the user must (1) determine the appropriate ISC value (for example, 

singular parameter set, ISC = -1, dual parameter set, ISC = -2), (2) insert the relationships into the corresponding FORTRAN 

FUNCTION, and (3) re-compile MODFLOW with the newly created file containing the custom set of relationships.
Appendix 3 57

C======================================================================= 

 FUNCTION SVSF1REF1PS(EFSAT,ALPHA,VGN,ISC) 

C-----------VERSION 20AUGUST2005 RBT 

C ******************************************************************

C CALCULATE PORE-WATER PRESSURE AS A FUNCTION OF EFFECTIVE 

C SATURATION USING VAN GENUCHTEN-MUALEM FUNCTION (ISC>1), 

C VAUCLIN FUNCTION (ISC=-2) OR HAVERKAMP LOG FUNCTION (ISC=-3) 

C ******************************************************************

C SPECIFICATIONS:

 REAL EFSAT,ALPHA,VGN,M,PRESS 

 INTEGER ISC 

C ******************************************************************

C 

 IF (ISC.GT.0) THEN 

 M=1-(1/VGN) 

 PRESS=(1./ALPHA)*(1./(EFSAT**(1/M))-1)**(1./VGN) 

 SVSF1REF1PS=-PRESS 

 ELSEIF(ISC.EQ.-2) THEN 

 PRESS=((ALPHA/EFSAT)-ALPHA)**(1/VGN) 

 SVSF1REF1PS=-PRESS 

 ELSE 

 P=((ALPHA/EFSAT)-ALPHA)**(1/VGN) 

 SVSF1REF1PS=-EXP(P) 

 ENDIF 

C 

 END 

C 

C======================================================================= 

 FUNCTION SVSF1REF1KP(PRESS,ALPHA,VGN,ISC) 

C-----------VERSION 20AUGUST2005 RBT 

C ******************************************************************

C CALCULATE RELATIVE PERMEABILITY AS A FUNCTION OF 

C PORE-WATER PRESSURE USING VAN GENUCHTEN FUNCTION (ISC>1) 

C OR VAUCLIN FUNCTION (ISC<0), 

C ******************************************************************

C SPECIFICATIONS:

 REAL PRESS,ALPHA,VGN,M 

 INTEGER ISC 

C ******************************************************************

C 

 IF (ISC.GT.0) THEN 

 M=1.-(1./VGN) 

 PRESS=-PRESS 

 TOP=(1-(ALPHA*PRESS)**(VGN-1)*(1+(ALPHA*PRESS)**VGN)**(-M))**2 

 SVSF1REF1KP=TOP/((1+(ALPHA*PRESS)**VGN)**(M/2)) 

 PRESS=-PRESS 

 ELSE 

 PRESS=-PRESS 

 SVSF1REF1KP=ALPHA/(ALPHA+(PRESS)**VGN) 

 PRESS=-PRESS 

 ENDIF 

C 

 END 

C 

C======================================================================= 

 FUNCTION SVSF1REF1SP(PRESS,RSAT,ALPHA,VGN,ISC) 

C-----------VERSION 20AUGUST2005 RBT 

C ******************************************************************

C CALCULATE SOIL WATER SATURATION AS A FUNCTION OF PORE-WATER 

C PRESSURE USING VAN GENUCHTEN-MUALEM FUNCTION (ISC>1), 

C VAUCLIN FUNCTION (ISC=-2), OR HAVERKAMP LOG FUNCTION (ISC=-3) 

C ******************************************************************

C SPECIFICATIONS:

 REAL PRESS,RSAT,ALPHA,VGN,M 

 INTEGER ISC 

C ******************************************************************

C 

 

Figure 39. Code listing of FORTRAN FUNCTIONS defining soil characteristic relations.
58 User’s Guide to the Variably Saturated Flow (VSF) Process for MODFLOW

 IF (ISC.GT.0) THEN 

 M=1.-(1./VGN) 

 PRESS=-PRESS 

 SVSF1REF1SP=RSAT+(1.-RSAT)/((1.+(ALPHA*PRESS)**VGN)**M) 

 PRESS=-PRESS 

 ELSEIF (ISC.EQ.-2) THEN 

 PRESS=-PRESS 

 SVSF1REF1SP=RSAT+(ALPHA*(1-RSAT))/(ALPHA+PRESS**VGN) 

 PRESS=-PRESS 

 ELSEIF (ISC.EQ.-3) THEN 

 PRESS=-PRESS 

 SVSF1REF1SP=RSAT+(ALPHA*(1-RSAT))/(ALPHA+LOG(PRESS)**VGN) 

 PRESS=-PRESS 

 ENDIF 

 END 

C 

C====================================================================== 

 FUNCTION SVSF1REF1CP(PRESS,RSAT,POR,ALPHA,VGN,ISC) 

C-----------VERSION 20AUGUST2005 RBT 

C *****************************************************************

C CALCULATE SPECIFIC MOISTURE CAPACITY AS A FUNCTION OF PORE-WATER 

C PRESSURE USING VAN GENUCHTEN FUNCTION (ISC>1), 

C VAUCLIN FUNCTION (ISC=-2), OR HAVERKAMP LOG FUNCTION (ISC=-3) 

C *****************************************************************

C SPECIFICATIONS 

 REAL SAT,RSAT,POR,ALPHA,VGN,M,EFSAT,BUFF,PRESS,RSWC 

 INTEGER ISC 

C ******************************************************************

C 

 IF (ISC.GT.0) THEN 

 M=1.-(1./VGN) 

 PRESS= -PRESS 

 EFSAT=(1./(1.+(PRESS*ALPHA)**VGN))**M 

 BUFF=(1./(1.-M))*ALPHA*M*(POR-RSAT)*(EFSAT**(1./M)) 

 SVSF1REF1CP=BUFF*((1.-EFSAT**(1./M))**M) 

 PRESS=-PRESS 

 ELSEIF (ISC.EQ.-2) THEN 

 PRESS=-PRESS 

 RSWC=RSAT*POR 

 BUFF=ALPHA*VGN*(POR-RSWC)*PRESS**(VGN-1) 

 SVSF1REF1CP=BUFF/((ALPHA+PRESS**VGN)**2.) 

 PRESS=-PRESS 

 ELSEIF (ISC.EQ.-3) THEN 

 PRESS=-PRESS 

 RSWC=RSAT*POR 

 BUFF=ALPHA*VGN*(POR-RSWC)*(LOG(PRESS))**(VGN-1) 

 SVSF1REF1CP=BUFF/(PRESS*((ALPHA+(LOG(PRESS))**VGN)**2.)) 

 PRESS=-PRESS 

 ENDIF 

C 

 END 

Figure 39. Code listing of FORTRAN FUNCTIONS defining soil characteristic relations.—Continued