Water Resources of the United States

Home News Getting Started Overview Gallery FAQ Glossary Special Support Software

Unsaturated flow functions



Preparing the main input file

When simulating saturated-unsaturated flow, SUTRA repeatedly evaluates the unsaturated flow properties at various points in the model by calling subroutine UNSAT, which must be programmed by the user. To allow UNSAT to set different properties in different regions of the model domain, each node and element is assigned an "unsaturated flow property region number." These region numbers are specified in datasets 14B (for nodes) and 15B (for elements). See Appendix B of the SUTRA documentation for details.

The following example shows a portion of dataset 14B in which nodes 1 - 5 are assigned to unsaturated flow property region 1, NREG = 1, and nodes 6 - 10 are assigned to region 2, NREG = 2 (where "#" indicates a comment):

# Dataset 14B
# [II]  [NREG]     [X] [Y]    [Z]    [POR]
     1       1      0.  3.  21.00     0.35 
     2       1      0.  3.  19.95     0.35 
     3       1      0.  3.  18.90     0.35 
     4       1      0.  3.  17.85     0.35 
     5       1      0.  3.  16.80     0.35 
     6       2      0.  3.  15.75     0.35 
     7       2      0.  3.  14.70     0.35 
     8       2      0.  3.  13.65     0.35 
     9       2      0.  3.  12.60     0.35 
    10       2      0.  3.  11.55     0.35 

The region numbers for elements, LREG, are assigned similarly in dataset 15B.

Whenever SUTRA calls UNSAT to evaluate the unsaturated flow properties, the region number of the node or element (NREG or LREG) is passed to UNSAT, where it is named KREG. If programmed to do so by the user, UNSAT can use the value of KREG to determine the unsaturated flow properties.

Programming subroutine UNSAT

The unsaturated flow functions must be programmed by the user in subroutine UNSAT. The subroutine is divided into three main sections that correspond to the three functions required by SUTRA:

  1. Saturation, SW, as a function of pressure
  2. Derivative of saturation with respect to pressure, DSWDP
  3. Relative permeability, RELK, as a function of pressure or saturation

Guidelines for customizing each section of UNSAT are given below. Note that the following information is automatically passed to UNSAT and can be used in formulating the unsaturated flow functions:

  • PRES = pressure
  • KREG = unsaturated region number for the node or element (NREG or LREG)

Note also that, by default, subroutine UNSAT is programmed for the Van Genuchten model [1], with separate (but identical) values of the model parameters in each of two unsaturated flow property regions. The user must either modify the model parameters or replace the unsaturated flow functions entirely, to suit the problem being solved.

Setting model parameters

In subroutine UNSAT, the Van Genuchten model [1] is programmed by default. The model has three parameters: SWRES, AA, and VN. The values of these parameters depend on the region number, KREG, of the node or element for which the unsaturated functions are being computed. In the default implementation, there are two regions, KREG = 1 and KREG = 2.

The values of SWRES, AA, and VN in region 1 are named SWRES1, AA1, and VN1, respectively, and are set in the following statements:

C     DATA FOR REGION 1:                                                 UNSAT.........3600
      DATA   SWRES1/0.30E0/,   AA1/5.0E-5/,   VN1/2.0E0/                 UNSAT.........3700
      SAVE SWRES1, AA1, VN1                                              UNSAT.........3800

The SAVE statement ensures that the parameter values are remembered between calls to UNSAT.  Similarly, the values of SWRES, AA, and VN in region 2 are named SWRES2, AA2, and VN2, respectively, and are set in the following statements:

C     DATA FOR REGION 2:                                                 UNSAT.........3900
      DATA   SWRES2/0.30E0/,   AA2/5.0E-5/,   VN2/2.0E0/                 UNSAT.........4000
      SAVE SWRES2, AA2, VN2                                              UNSAT.........4100

Before the unsaturated flow functions are computed, the values of the model parameters SWRES, AA, and VN are set depending on the region number, KREG, as follows (with additional comments in red):

C     SET PARAMETERS FOR CURRENT REGION, KREG                            UNSAT.........5500
      GOTO(10,20),KREG                                                   UNSAT.........5600
   10 SWRES=SWRES1       ! If KREG=1, set parameters to region 1 values  UNSAT.........5700
      AA=AA1                                                             UNSAT.........5800
      VN=VN1                                                             UNSAT.........5900
      GOTO 100                                                           UNSAT.........6000
   20 SWRES=SWRES2       ! If KREG=2, set parameters to region 2 values  UNSAT.........6100
      AA=AA2                                                             UNSAT.........6200
      VN=VN2                                                             UNSAT.........6300
  100 CONTINUE                                                           UNSAT.........6400

The sections of code described above illustrate one method for setting the parameters for a particular set of unsaturated flow functions (i.e., the Van Genuchten model).  They can be replaced with whatever code is appropriate for the unsaturated flow functions you want to use.

Note that if you have modified UNSAT, you must recompile the SUTRA code using a standard FORTRAN-90 compiler before running it.

Section (1): Saturation as a function of pressure

The section of subroutine UNSAT that computes the value of the saturation as a function of pressure reads as follows (with additional comments in red):

C*********************************************************************** UNSAT.........6700
C*********************************************************************** UNSAT.........6800
C.....SECTION (1):                                                       UNSAT.........6900
C     SW VS. PRES   (VALUE CALCULATED ON EACH CALL TO UNSAT)             UNSAT.........7000
C     CODING MUST GIVE A VALUE TO SATURATION, SW.                        UNSAT.........7100
C                                                                        UNSAT.........7200
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT.........7300
C     THREE PARAMETER MODEL OF VAN GENUCHTEN(1980)                       UNSAT.........7400
      SWRM1=1.E0-SWRES                                                   UNSAT.........7500
      AAPVN=1.E0+(AA*(-PRES))**VN                                        UNSAT.........7600
      VNF=(VN-1.E0)/VN                                                   UNSAT.........7700
      AAPVNN=AAPVN**VNF                                                  UNSAT.........7800
      S W   =   DBLE (SWRES+SWRM1/AAPVNN)    ! Set value of SW           UNSAT.........7900
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT.........8000
C*********************************************************************** UNSAT.........8100
C*********************************************************************** UNSAT.........8200

The Van Genuchten model [1] is programmed by default. There are several lines of code, but the ultimate result is that the saturation, SW, is computed based on the pressure, PRES, and the model parameters SWRES, AA, and VN. The values of the model parameters depend on the region number, KREG, and are set in an earlier section of the code. To program a different function for SW, replace the code that sets the model parameters and the code in Section (1) with appropriate expressions.

Note that if you have modified UNSAT, you must recompile the SUTRA code using a standard FORTRAN-90 compiler before running it.

Section (2): Derivative of saturation with respect to pressure

The section of subroutine UNSAT that computes the value of the derivative of saturation with respect to pressure reads as follows (with additional comments in red):

C*********************************************************************** UNSAT.........9000
C*********************************************************************** UNSAT.........9100
C.....SECTION (2):                                                       UNSAT.........9200
C     DSWDP VS. PRES, OR DSWDP VS. SW   (CALCULATED ONLY WHEN IUNSAT=1)  UNSAT.........9300
C     CODING MUST GIVE A VALUE TO DERIVATIVE OF SATURATION WITH          UNSAT.........9400
C     RESPECT TO PRESSURE, DSWDP.                                        UNSAT.........9500
C                                                                        UNSAT.........9600
  600 CONTINUE                                                           UNSAT.........9700
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT.........9800
      DNUM=AA*(VN-1.E0)*SWRM1*(AA*(-PRES))**(VN-1.E0)                    UNSAT.........9900
      DNOM=AAPVN*AAPVNN                                                  UNSAT........10000
      D S W D P   =   DBLE (DNUM/DNOM)     ! Set value of DSWDP          UNSAT........10100
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT........10200
      GOTO 1800                                                          UNSAT........10300
C*********************************************************************** UNSAT........10400
C*********************************************************************** UNSAT........10500

The Van Genuchten model [1] is programmed by default. There are several lines of code, but the ultimate result is that the derivative of saturation with respect to pressure, DSWDP, is computed based on the pressure, PRES, and the model parameters SWRES, AA, and VN. [Note that SWRM1 = 1. - SWRES, AAPUN = 1. + (-AA*PRES)**VN, and AAPUNN = AAPUN**VN.] The values of the model parameters depend on the region number, KREG, and are set in an earlier section of the code. To program a different function for DSWDP, replace the code that sets the model parameters and the code in Section (2) with appropriate expressions. If desired, the function for DSWDP can be formulated in terms of saturation, SW, instead of pressure, PRES.

Note that if you have modified UNSAT, you must recompile the SUTRA code using a standard FORTRAN-90 compiler before running it.

Section (3): Relative permeability as a function of pressure or saturation

The section of subroutine UNSAT that computes the value of the relative permeability reads as follows (with additional comments in red):

C*********************************************************************** UNSAT........11200
C*********************************************************************** UNSAT........11300
C.....SECTION (3):                                                       UNSAT........11400
C     RELK VS. P, OR RELK VS. SW   (CALCULATED ONLY WHEN IUNSAT=2)       UNSAT........11500
C     CODING MUST GIVE A VALUE TO RELATIVE PERMEABILITY, RELK.           UNSAT........11600
C                                                                        UNSAT........11700
 1200 CONTINUE                                                           UNSAT........11800
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT........11900
C     GENERAL RELATIVE PERMEABILITY MODEL FROM VAN GENUCHTEN(1980)       UNSAT........12000
      SWSTAR=(SW-SWRES)/SWRM1                                            UNSAT........12100
      R E L K   =   DBLE (SQRT(SWSTAR)*        ! Set value of RELK       UNSAT........12200
     1                   (1.E0-(1.E0-SWSTAR**(1.E0/VNF))**(VNF))**2)     UNSAT........12300
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  UNSAT........12400
C                                                                        UNSAT........12500
C*********************************************************************** UNSAT........12600
C*********************************************************************** UNSAT........12700

The Van Genuchten model [1] is programmed by default. There are several lines of code, but the ultimate result is that the relative permeability, RELK, is computed based on the saturation, SW, and the model parameters SWRES and VN. [Note that SWRM1 = 1. - SWRES, and VNF = (1. - VN)/VN]. The values of the model parameters depend on the region number, KREG, and are set in an earlier section of the code. To program a different function for RELK, replace the code that sets the model parameters and the code in Section (3) with appropriate expressions. If desired, the function for RELK can be formulated in terms of pressure, PRES, instead of saturation, SW.

Note that if you have modified UNSAT, you must recompile the SUTRA code using a standard FORTRAN-90 compiler before running it.

References

[1]  Van Genuchten, M.Th., 1980, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils: Soil Science Society of America Journal, v. 44, no. 5, p. 892-898.

 


USGS Home :: Biology :: Geology :: Geography :: Water Intranet :: USGS Intranet :: Site Map

Accessibility FOIA Privacy Policies and Notices

Take Pride in America logo USA.gov logo U.S. Department of the Interior | U.S. Geological Survey
URL: http://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm
Page Contact Information: Alden Provost <aprovost@usgs.gov>
Page Last Modified: Monday, 07-Jan-2013 11:58:07 EST