Water Resources of the United States
Unsaturated flow functionsPreparing the main input fileWhen 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 UNSATThe 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:
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:
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 parametersIn 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 pressureThe 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 pressureThe 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 saturationThe 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
U.S. Department of the Interior | U.S. Geological Survey
URL: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm
Page Contact Information: Alden Provost
<aprovost@usgs.gov>
Page Last Modified: Tuesday, 20-Dec-2016 15:22:40 EST
