SutraSuite

https://water.usgs.gov/nrp/gwsoftware/sour/
Admin
SutraSuite is a groundwater flow and transport modeling software package featuring the USGS SUTRA...
Tokens:12,403
Snippets:144
Trust Score:10
Update:2 months ago
Show doc for...
Tokens:
Raw
### Implementing Time-Dependent Boundary Conditions

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

Example implementation for setting time-dependent injection and withdrawal rates based on the simulation year (TYEAR).

```Fortran
C     NOTE: A FLOW AND TRANSPORT SOLUTION MUST OCCUR FOR ANY            BCTIME.......15700
C     TIME STEP IN WHICH QIN( ) CHANGES.                                BCTIME.......15800
      IF (TYEAR.LT.1980.) THEN   ! Before 1980.
         IF (-I.EQ.652) THEN
            QIN(-I) = 10.           ! Injection rate = 10 kg/s.
         ELSE IF (-I.EQ.798) THEN
            QIN(-I) = -10.          ! Withdrawal rate = 10 kg/s.
         END IF
      ELSE                       ! During and after 1980.
         IF (-I.EQ.652) THEN
            QIN(-I) = 30.           ! Injection rate = 30 kg/s.
         ELSE IF (-I.EQ.798) THEN
            QIN(-I) = -30.          ! Withdrawal rate = 30 kg/s.
         END IF
      END IF
C     NOTE: A TRANSPORT SOLUTION MUST OCCUR FOR ANY                     BCTIME.......16000
C     TIME STEP IN WHICH UIN( ) CHANGES.                                BCTIME.......16100
      UIN(-I) = 0.                  ! Solute conc of injectate = 0.
```

--------------------------------

### Implement Time-Dependent Temperature Boundary Conditions

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

Example implementation for a temperature rising at a constant rate per year.

```Fortran
      TRISE = 2.*TYEAR         ! Temp rising at 2 deg-C/yr
      UBC(IUP) = 18. + TRISE   ! Set time-dependent temperature
```

--------------------------------

### Implementing a Time-Dependent Pulse Injection

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

Example implementation for injecting a solute pulse at node 58 for one hour starting at t=2 hours.

```Fortran
      IF (-I.EQ.58) THEN
         INJECT = ((THOUR.GE.2.).AND.(THOUR.LT.3.))
         IF (INJECT) THEN  ! Within injection period.
            QUIN(-I) = 1E-5    ! Inject solute at 0.00001 kg/s.
         ELSE              ! Not within injection period.
            QUIN(-I) = 0.      ! Do not inject solute.
         END IF
      END IF
```

--------------------------------

### Implement Time-Dependent Pressure Boundary Conditions

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

Example implementation for hydrostatic pressure at the bottom of a rising lake. Requires recompilation of the SUTRA code after modification.

```Fortran
      STAGE = 50. + 0.01*TDAY     ! Lake stage rising at 0.01 m/day
      DEPTH = STAGE - Z(-I)       ! Depth = stage - (node elevation)
      PBC(IP) = 1000.*GRAVZ*DEPTH ! Hydrostatic P = (density)*g*(depth)
      UBC(IP) = 0.                ! Lake water solute conc = 0.
```

--------------------------------

### Define Boundary Conditions in .inp File

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

Example of specifying a constant-valued boundary condition and a time-dependent boundary condition (via BCTIME) in dataset 17 of the .inp file.

```text
101  0.001  0.0  # Constant-valued boundary condition at node 101
-102  0.000  0.0  # Time-dependent boundary condition at node 102 (defined in BCTIME)
```

--------------------------------

### Assigning Unsaturated Flow Property Regions in Dataset 14B

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Example of assigning unsaturated flow property region numbers to nodes in the SUTRA input file.

```text
# Dataset 14B
```

--------------------------------

### BCTIME Subroutine Specified Concentration/Temperature Logic

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

The default structure for setting specified-concentration or temperature boundary conditions in BCTIME.

```Fortran
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......12200
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......12300
C.....SECTION (2): SET TIME-DEPENDENT SPECIFIED                         BCTIME.......12400
C     CONCENTRATIONS (TEMPERATURES)                                     BCTIME.......12500
C                                                                       BCTIME.......12600
  250 CONTINUE                                                          BCTIME.......12700
      DO 400 IU=1,NUBC   ! Loop over all specified-conc/temp nodes.     BCTIME.......12800
      IUP=IU+NPBC                                                       BCTIME.......12900
      I=IUBC(IUP)        ! Set I = number of current node.              BCTIME.......13000
      IF(I) 300,400,400  ! If I is negative, set boundary conditions;   BCTIME.......13100
  300 CONTINUE           !    otherwise, skip to the next node.         BCTIME.......13200
C     NOTE: A TRANSPORT SOLUTION MUST OCCUR FOR ANY TIME STEP IN WHICH  BCTIME.......13300
C     UBC( ) CHANGES. IN ADDITION, IF FLUID PROPERTIES ARE              BCTIME.......13400
C     SENSITIVE TO 'U', THEN A FLOW SOLUTION MUST OCCUR AS WELL.        BCTIME.......13500
C     UBC(IUP) = (( ))   ! Enter expression for conc/temp here.         BCTIME.......13600
400   CONTINUE                                                          BCTIME.......13700
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......13800
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......13900
```

--------------------------------

### Calculate Saturation in UNSAT

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Computes the saturation value (SW) based on pressure (PRES) and Van Genuchten model parameters.

```Fortran
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
```

--------------------------------

### Calculate Derivative of Saturation in UNSAT

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Computes the derivative of saturation with respect to pressure (DSWDP) using the Van Genuchten model parameters.

```Fortran
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
```

--------------------------------

### Define Van Genuchten Parameters for Region 1

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Sets the default model parameters for the first unsaturated flow region.

```Fortran
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
```

--------------------------------

### BCTIME Solute Mass or Energy Source/Sink Loop

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

The standard loop structure in BCTIME for iterating over solute/energy source/sink nodes and identifying nodes requiring boundary condition updates.

```Fortran
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......17300
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......17400
C.....SECTION (4): SET TIME-DEPENDENT SOURCES/SINKS                     BCTIME.......17500
C     OF SOLUTE MASS OR ENERGY                                          BCTIME.......17600
C                                                                       BCTIME.......17700
  650 CONTINUE                                                          BCTIME.......17800
      DO 800 IQU=1,NSOUI  ! Loop over solute/energy source/sink nodes.  BCTIME.......17900
      I=IQSOU(IQU)        ! Set I = number of current node.             BCTIME.......18000
      IF(I) 700,800,800   ! If I is negative, set boundary conditions;  BCTIME.......18100
      700 CONTINUE        !    otherwise, skip to the next node.        BCTIME.......18200
C     NOTE: A TRANSPORT SOLUTION MUST OCCUR FOR ANY                     BCTIME.......18300
C     TIME STEP IN WHICH QUIN( ) CHANGES.                               BCTIME.......18400
C     QUIN(-I) = (( ))    ! Enter expression for source/sink rate here. BCTIME.......18500
  800 CONTINUE                                                          BCTIME.......18600
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......18700
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......18800
```

--------------------------------

### Define Van Genuchten Parameters for Region 2

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Sets the default model parameters for the second unsaturated flow region.

```Fortran
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
```

--------------------------------

### Relative Permeability Calculation in UNSAT

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Computes relative permeability using the Van Genuchten model based on saturation and model parameters.

```Fortran
 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
```

--------------------------------

### Set Time-Dependent Specified Pressures in BCTIME

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

This section of the BCTIME subroutine iterates through specified-pressure nodes and applies user-defined expressions for pressure and inflow concentration or temperature when the node index is negative.

```Fortran
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME........9700
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME........9800
C.....SECTION (1): SET TIME-DEPENDENT SPECIFIED PRESSURES OR            BCTIME........9900
C     CONCENTRATIONS (TEMPERATURES) OF INFLOWS AT SPECIFIED             BCTIME.......10000
C     PRESSURE NODES                                                    BCTIME.......10100
C                                                                       BCTIME.......10200
   50 CONTINUE                                                          BCTIME.......10300
      DO 200 IP=1,NPBC   ! Loop over all specified-pressure nodes.      BCTIME.......10400
      I=IPBC(IP)         ! Set I = number of current node.              BCTIME.......10500
      IF(I) 100,200,200  ! If I is negative, set boundary conditions;   BCTIME.......10600
  100 CONTINUE           !    otherwise, skip to the next node.         BCTIME.......10700
C     NOTE: A FLOW AND TRANSPORT SOLUTION MUST OCCUR FOR ANY            BCTIME.......10800
C     TIME STEP IN WHICH PBC( ) CHANGES.                                BCTIME.......10900
C     PBC(IP) = (( ))    ! Enter expression for pressure here.          BCTIME.......11000
C     UBC(IP) = (( ))    ! Enter expression for inflow conc/temp here.  BCTIME.......11100
  200 CONTINUE                                                          BCTIME.......11200
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......11300
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......11400
```

--------------------------------

### BCTIME Fluid Source/Sink Loop Structure

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm

The standard loop structure for iterating over fluid source/sink nodes in the BCTIME subroutine. Nodes with negative indices are processed for boundary condition updates.

```Fortran
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......14700
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......14800
C.....SECTION (3): SET TIME-DEPENDENT FLUID SOURCES/SINKS,              BCTIME.......14900
C     OR CONCENTRATIONS (TEMPERATURES) OF SOURCE FLUID                  BCTIME.......15000
C                                                                       BCTIME.......15100
  450 CONTINUE                                                          BCTIME.......15200
      DO 600 IQP=1,NSOPI  ! Loop over all fluid source/sink nodes.      BCTIME.......15300
      I=IQSOP(IQP)        ! Set I = number of current node.             BCTIME.......15400
      IF(I) 500,600,600   ! If I is negative, set boundary conditions;  BCTIME.......15500
  500 CONTINUE            !    otherwise, skip to the next node.        BCTIME.......15600
C     NOTE: A FLOW AND TRANSPORT SOLUTION MUST OCCUR FOR ANY            BCTIME.......15700
C     TIME STEP IN WHICH QIN( ) CHANGES.                                BCTIME.......15800
C     QIN(-I) = (( ))     ! Enter expression for source/sink rate here. BCTIME.......15900
C     NOTE: A TRANSPORT SOLUTION MUST OCCUR FOR ANY                     BCTIME.......16000
C     TIME STEP IN WHICH UIN( ) CHANGES.                                BCTIME.......16100
C     UIN(-I) = (( ))     ! Enter expression for inflow conc/temp here. BCTIME.......16200
  600 CONTINUE                                                          BCTIME.......16300
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......16400
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - BCTIME.......16500
```

--------------------------------

### Select Parameters Based on Region Number

Source: https://water.usgs.gov/nrp/gwsoftware/sour/special/unsat/unsat.htm

Uses the KREG value to assign the appropriate model parameters for the current calculation.

```Fortran
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
```
Tokens:
Raw