Water Resources of the United States
Time-dependent boundary conditions
Two ways of specifying time-dependent boundary conditionsSUTRA 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 fileIn 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 filesThe 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 BCTIMEThe subroutine is divided into four main sections that correspond to the four types of boundary conditions supported by SUTRA:
Each of the four sections of code performs a similar sequence of steps:
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:
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 pressuresThe 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 temperaturesThe 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 massThe 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 energyThe 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