Water Resources of the United States

Home News Getting Started Overview Gallery FAQ Glossary Special Support Software

Time-dependent boundary conditions



Two ways of specifying time-dependent boundary conditions

SUTRA Version 2.1 and earlier required that time-dependent boundary conditions be implemented by programming subroutine BCTIME. Beginning with Version 2.2, time-dependent boundary conditions can alternatively be specified using one or more (".bcs") input files; no programming is required.

When time-dependent boundary conditions are specified via subroutine BCTIME, the SUTRA code must be recompiled using a Fortran 90 compiler. Time-dependent boundary conditions specified via the optional ".bcs" files do not require any code modifications, and they take precedence over corresponding specifications in the ".inp" file or subroutine BCTIME.

Preparing the main input file

In SUTRA, boundary conditions are applied at nodes. All boundary condition nodes must be listed in datasets 17 - 20 of the main input (".inp") file, whether or not their specifications change with time. See Appendix B of the SUTRA documentation for details.

For a boundary condition that does not vary with time (e.g., withdrawal of ground water at a constant rate throughout the entire simulation), or whose time-dependence is defined using only ".bcs" files (not subroutine BCTIME), the specification in the ".inp" file consists of the node number followed by the constant, default values of the parameters associated with that boundary condition. For a boundary condition whose time-dependence includes specifications programmed in subroutine BCTIME, the specification in the ".inp" file consists only of the node number preceded by a minus sign. The following example shows one constant-valued and one time-dependent boundary condition (defined via BCTIME) in dataset 17 (where "#" indicates a comment):

# Dataset 17: Sources/sinks of fluid mass
#
 417  4.75E-3   # Constant recharge to node 417 at 4.75E-3 kg/s
-543            # Time-dependent source/sink at node 543

In this example, the default specifications are that water is injected at a constant rate of 4.75x10-3 kg/s at node 417, and injected or withdrawn at node 543 at a time-dependent rate set in subroutine BCTIME (as indicated by the negative node number). If a source/sink of fluid mass is specified at node 417 in a ".bcs" file on a certain time step, the ".bcs" specification will supersede the default ".inp" specifications beginning on that time step. Similarly, if a source/sink of fluid mass is specified at node 543 in a ".bcs" file on a certain time step, the ".bcs" specification will supersede any fluid source/sink specification made in subroutine BCTIME up to that point.

Preparing the ".bcs" input files

The optional ".bcs" files need contain only those specifications that the user wishes change on any given time step. Specifications made in the ".bcs" files are persistent: they remain in force unless/until changed by a subsequent ".bcs" specification. The input formats for ".bcs" datasets parallel those for ".inp" datasets 17 - 20. Any boundary condition node listed in a ".bcs" dataset must also be listed (with either a positive or negative node number) in the corresponding ".inp" dataset. See Appendix B of the SUTRA documentation for details.

Programming subroutine BCTIME

The subroutine is divided into four main sections that correspond to the four types of boundary conditions supported by SUTRA:

  1. Specified pressures
  2. Specified concentrations or temperatures
  3. Sources/sinks of fluid mass
  4. Sources/sinks of solute mass or energy

Each of the four sections of code performs a similar sequence of steps:

  • Loop through the list of nodes to which a particular type of boundary condition has been assigned (i.e., the nodes listed in dataset 17, 18, 19, or 20).
  • For each node designated with a negative node number, set the user-programmed boundary conditions.

Guidelines for customizing each section of BCTIME are given below. Note that the following information is automatically passed to BCTIME and can be used in formulating the time-dependent boundary conditions:

  • IT = number of the current time step
  • TSEC = time at the end of the current time step, in seconds
  • TMIN = time at the end of the current time step, in minutes
  • THOUR = time at the end of the current time step, in hours
  • TDAY = time at the end of the current time step, in days
  • TWEEK = time at the end of the current time step, in weeks
  • TMONTH = time at the end of the current time step, in months
  • TYEAR = time at the end of the current time step, in years
  • X = array of x-coordinates of nodes
  • Y = array of y-coordinates of nodes
  • Z = array of z-coordinates of nodes
  • GRAVX = x-component of gravity
  • GRAVY = y-component of gravity
  • GRAVZ = z-component of gravity

Also note that the simple examples given below are structured to maximize clarity of exposition rather than computational efficiency. Depending on the nature and complexity of the boundary conditions, it may or may not be possible to improve computational efficiency by streamlining or rearranging the sequence of operations. 

Section (1): Specified pressures

The section of subroutine BCTIME in which specified-pressure boundary conditions are set reads as follows (with additional comments in red):

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

The code loops over all NPBC specified-pressure nodes, whose node numbers are listed in array IPBC. For each node in the list, the sign of the node number is checked. If the node number is negative, the boundary pressure, PBC, and the associated inflow concentration or temperature, UBC, are set for that node.  If the node number is positive, the code simply moves on to the next node in the list. 

The two lines in which the values of PBC and UBC are set are initially commented out (using "C" as the first character of each line). To implement time-dependent boundary conditions, you must uncomment these two lines (replace each leading "C" with a space), replace the symbol "(( ))" with an appropriate expression on each line, and supply any additional calculations that are needed. For example, to set time-dependent boundary conditions that represent hydrostatic pressure at the bottom of a fresh-water lake whose stage is rising at 0.01 m/day from an initial level of 50 m above datum, one could replace the two lines just discussed with the following lines:

      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.

This example makes use of TDAY, GRAVZ, and the array Z, which are automatically passed to BCTIME by SUTRA. For convenience, it also defines two new variables, STAGE and DEPTH.

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

Section (2): Specified concentrations or temperatures

The section of subroutine BCTIME in which specified-concentration or temperature boundary conditions are set reads as follows (with additional comments in red):

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

The code loops over all NUBC specified-concentration or temperature nodes, whose node numbers are listed in array IUBC (in entries NPBC+1 through NPBC+1+NUBC). For each node in the list, the sign of the node number is checked. If the node number is negative, the boundary concentration or temperature, UBC, is set for that node.  If the node number is positive, the code simply moves on to the next node in the list.

The lines in which the value of UBC is set is initially commented out (using "C" as the first character of the line). To implement time-dependent boundary conditions, you must uncomment this line (replace the leading "C" with a space), replace the symbol "(( ))" with an appropriate expression, and supply any additional calculations that are needed. For example, to set time-dependent boundary conditions that represent a specified temperature that rises 2 °C/yr from an initial temperature of 18 °C at time t = 0 yrs, one could replace the line just discussed with the following lines:

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

This example makes use of TYEAR, which is automatically passed to BCTIME by SUTRA. For convenience, it also defines the new variable TRISE.

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

Section (3): Sources/sinks of fluid mass

The section of subroutine BCTIME in which fluid mass source/sink boundary conditions are set reads as follows (with additional comments in red):

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

The code loops over all NSOPI fluid mass source/sink nodes, whose node numbers are listed in array IQSOP. For each node in the list, the sign of the node number is checked.  If the node number is negative, the source/sink rate, QIN, and the associated inflow concentration or temperature, UIN, are set for that node. If the node number is positive, the code simply moves on to the next node in the list.

The lines in which the values of QIN and UIN are set are initially commented out (using "C" as the first character of each line). To implement time-dependent boundary conditions, you must uncomment these two lines (replace each leading "C" with a space), replace the symbol "(( ))" with an appropriate expression on each line, and supply any additional calculations that are needed. For example, to set time-dependent boundary conditions that represent injection of fresh water at node 652 and withdrawal from node 798, at different constant rates before and after the start of year 1980, one could replace the two lines just discussed, and their accompanying comment lines, with the following:

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.

This example makes use of TYEAR, which is automatically passed to BCTIME by SUTRA.

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

Section (4): Sources/sinks of solute mass or energy

The section of subroutine BCTIME in which solute mass or energy source/sink boundary conditions are set reads as follows (with additional comments in red):

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

The code loops over all NSOUI solute mass or energy source/sink nodes, whose node numbers are listed in array IQSOU. For each node in the list, the sign of the node number is checked.  If the node number is negative, the source/sink rate, QUIN, is set for that node. If the node number is positive, the code simply moves on to the next node in the list.

The line in which the value of QUIN is set is initially commented out (using "C" as the first character of the line). To implement time-dependent boundary conditions, you must uncomment this line (replace the leading "C" with a space), replace the symbol "(( ))" with an appropriate expression, and supply any additional calculations that are needed. For example, to set a time-dependent boundary condition that represents injection of a pulse of solute at node 58 for one hour, starting at time t = 2 hours, one could replace the line just discussed with the following:

      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

This example makes use of THOUR, which is automatically passed to BCTIME by SUTRA. For convenience, it also defines a new logical variable, INJECT.

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

 


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: https://water.usgs.gov/nrp/gwsoftware/sour/special/bctime/bctime.htm
Page Contact Information: Alden Provost <aprovost@usgs.gov>
Page Last Modified: Tuesday, 20-Dec-2016 15:22:41 EST