Time-Variant Conductance Package:
Documentation of a Computer Program to Include Horizontal Wells, Multiaquifer Wells, and Faults in MODFLOW
Richard B. Winston
2163 Tulip St.
Baton Rouge, La, 70806
Table of Contents
Table of Contents*
Conceptualization and Implementation*
Explanation of Parameters Read by the TVC1 package.*
Sample Input and Output*
Modified Main Program*
The Time-Variant Conductance Package for MODFLOW enables the user to vary the conductance between specific cells in each stress period. This capability is useful for more easily modeling active faults, horizontal wells and multiaquifer wells. The latter two applications can be implemented by specifying high conductances between the cells penetrated by the horizontal or multiaquifer well for those stress periods for which the well is present. Pumpage for the well would then be specified using the Well package. Pumpage could be withdrawn from any of the cells making up the horizontal or multiaquifer well and the high conductance between cells would ensure that the difference in head between cells in the well would be low.
Conceptualization and Implementation
The Time-Variant Conductance package (TVC1) is intended to address groundwater flow problems in which the hydraulic properties of a few grid cells change as a function of time. With the most recently released public version of MODFLOW (Harbaugh and McDonald, 1996a, 1996b) this could only by accomplished by creating a series of models that have different hydraulic properties and then using the final heads computed by one model as the starting heads of another model. Although feasible, this is hardly convenient if only a few changes need to be made.
There are at least three situations for which it would be convenient to be able to vary the conductance between cells as a function of time: (1) active faults which change their conductance as a function of time (2) horizontal wells, and (3) multiaquifer wells. In both of the latter two applications, the conductance between the cells that make up the well can be changed to a value sufficiently high to ensure that all those cells maintain a nearly identical head. Pumpage from any one of the cells would then, in effect, be applied to all of the cells taken together. The model would then take care of computing the appropriate flow into or out of each cell making up the well. Even in a homogeneous aquifer, this might result in differing flows to different cells because of the differences in those cells positions. For example, in a long horizontal well, in a single uniform aquifer, flow toward the center of the well would be nearly perpendicular to the well whereas flow toward the ends of the wells, water would flow radially. The end result would be that the flow into the well per unit length would be greater near the ends of the well than near the center. To model a horizontal well with the current public version of MODFLOW would require the modeler to compute the appropriate pumpage for each cell and place a well with that pumpage in each cell making up the horizontal well. The Time-Variant Conductance package allows the modeler to avoid this onerous task.
In the Block Centered Flow package (BCF5, McDonald and Harbaugh,1988; McDonald and others, 1992; Goode and Appel, 1992; Harbaugh and McDonald 1996a, 1996b) horizontal conductance along rows is calculated as an average transmissivity times cell height divided by the distance between the cells. In cases where hydraulic conductivity rather than transmissivity is specified, the trasmissivity is first calculated from the hydraulic conductivity and saturated thickness (McDonald and Harbaugh,1988). Several different methods are available for calculating the average the hydraulic conductivity (Goode and Appel, 1992). The horizontal conductance along columns is calculated by multiplying the horizontal conductance along rows by the anisotropy. The vertical conductance is calculated by multiplying the input parameter "Vcont" by the horizontal area of the cell.
The Time-Variant Conductance package uses the same formulas as the Block Centered Flow package (BCF5) to compute conductance between adjacent cells either horizontally or vertically. However, unlike the Block Centered Flow package, the conductances are only calculated between specified cells. The modeler must specify a pair of adjacent cells and give a hydraulic characteristic for each. The hydraulic characteristic will be either hydraulic conductivity, transmissivity or Vcont depending on the type of layer and the relative positions of the cells. TVC1 will then apply those hydraulic characteristics to calculate the correct conductance between those cells. For transmissivities (layer types 0 and 2 in the Block Centered Flow Package) and vertical conductances (VCONT), those changes will remain in effect for the remainder of the model or until explicitly changed. Thus, if the modeler wishes to restore the original model conductances in a later stress period for layer types 0 and 2, TVC1 must be used to modify the conductances between those cells whose original conductance are to be restored. In addition, if any horizontal flow barriers (Hsieh and Freckleton, 1993) were implemented, TVC1 does not restore them directly. Therefore, if you wish to restore a horizontal flow barrier, you must do so by using TVC1 to specify hydraulic properties of the cells in such a way that the correct conductance between them, including the effects of the horizontal flow barrier, is restored.
TVC1 can be used to implement variable anisotropy within a layer. However, most users would probably not find it convenient for that purpose because the size of the input file would be large compared with a package specifically designed for implementing variable anisotropy such as the one by Kladias and Ruskauff (1997). The present version of TVC1 can not by used with Kladias and Ruskauff’s variable anisotropy package because of differences in the way anisotropy is stored. However, the input file reserves spaces for uses an option to indicate that Kladias and Ruskauff’s variable anisotropy package is being used. A future version of the Transient Leakage Package will work correctly with either the USGS version of the Block-Centered Flow Package or Kladias and Ruskauff’s variable anisotropy package. At present, if you wish to use TVC1 with Kladias and Ruskauff’s variable anisotropy package, the source code of TVC1 must be modified to access the anisotropy array correctly or to eliminate all references to the anisotropy array. The latter approach would have the advantage of making TVC1 compatible with both BCF5 and Kladias and Ruskauff’s variable anisotropy package but would also require the user to incorporate the effects of anisotropy into the hydraulic conductivities or transmissivities specified in TVC1.
Models with large differences in conductance between cells may be insoluble because rounding errors in flow between cells with large conductances may be larger than the flows between cells with small conductances. With single precision arithmetic, as in MODFLOW, this occurs when flows between cells vary by a factor of approximately 108. If transmissivities or hydraulic conductivities vary by more than 108 differences in flows may make the model insoluble. For example, if a transmissivity, specified by TVC1 is 10 orders of magnitude larger than those of other cells in the same layer, the model will likely prove insoluble. Thus, when using TVC1 to model horizontal wells and multiaquifer wells and the hydraulic characteristic must be high enough to keep the head in the well the same in all cells comprising the well but still low enough not to cause numerical problems.
To test that TVC1 was implemented correctly, the results of several models using TVC1 were compared to models that did not use TVC1 but which were designed to have the same net effect. These tests are used as sample input for TVC1 in a following section. It was used to test that TVC1 did change the conductivity of selected cells in the intended manner. The model had twenty rows, twenty columns, and two layers. All layers were confined. Transmissivity of each layer was uniform, and set to 10-4 m2/s. The primary storage parameter of both layers was 10-3. VCONT (vertical hydraulic conductivity divided by thickness) was 10-1 s-1. The layers were isotropic with initial heads of 1000 m. Grid cells were 100 m x 100 m in size. There were two stress periods each 1 year in length and divided into 40 time steps of equal length. All the cells in the leftmost column in both layers were constant head cells. There was a single well in layer 2 column 10, row 20 with a pumping rate of –1.0 x 10-1 m3/s in both stress periods. SIP was used as a solver with a closure criterion of 1.0 x 10-4 m.
TVC1 was used in both stress periods to set the transmissivity along the 10’th column. However, in the first stress period, all the specified transmissivities were identical to those specified in the Block-Centered Flow package. In the second stress period, the transmissivities set with TVC1 were 10-1 m2/s: 1000 times larger than those set by with the Block-Centered Flow package. This was intended to simulate a horizontal well along the tenth column during the second stress period.
For comparison, the TVC1 package from the model. For the first stress period, this model was expected to produce results identical to the model using TVC1 for the first stress period but to diverge from it in the second stress period.
The second test was intended to show that TVC1 could correctly simulate a horizontal well. The second test of TVC1, was identical to the first test except that the transmissivities set by TVC1 were 1 m2/s and there was only one stress period in the model. For comparison, TVC1 was eliminated and wells were placed in all twenty cells in the 10’th column of layer 2. Each such well had a discharge of –5.0 x 10-3 m3/s: 1/20’th the discharge of the single well in the model with TVC1. (In this particular model, the withdrawal from each cell making up the horizontal well in the model with TVC1 was expected to be nearly the same because of symmetry.)
The third test was intended to show that TVC1 could correctly simulate a multiaquifer well. The third test of TVC1, was identical to the second test except that the instead of horizontal transmissivity, Vertical conductance was set by TVC1 between the cells in column 10, row 20, layer 1 and in column 10, row 20, layer 2. The vertical conductance used was 1.0 x 10-5 s-1 For comparison, TVC1 was eliminated and wells were placed in both cells that were connected with TVC1 in the test model. Each such well had a discharge of –5.0 x 10-2 m3/s: one half the discharge of the single well in the model with TVC1.
All three of the tests had the expected result: the models that did not use TVC1 had the same or nearly the same result as those which did for the first stress period. In the second stress period in model 1, the model with and without TVC1 diverged as was expected. This was intended to show that TVC1 was actually changing the conductances.
In the first test, there was no difference in the recorded heads at the end of stress period 1. Heads ranged from –82.39 to +1000 m. In the second test, the maximum difference in the recorded heads in time step 40 was 0.6 m. Heads ranged from 847.3 to 1000 m in layer 1 and from 713.0 to 1000 m in layer 2. In the third test, the maximum difference in the recorded heads in time step 40 was 0.3 m. Heads ranged from 375.7 to 1000 m with heads in both layers being nearly equal.
Input to the Time-Variant Conductance (TVC) Package is read from the file that has file type "TVC".
For each simulation
For each stress period
I10, I10, I10, F10.0 I10, I10, I10, F10.0
Input Item 3 consists of one record for each pair of adjacent cells between which the conductance will be altered in the current stress period. If ITEMP is negative or 0 item 3 is skipped. However any hydraulic characteristics which are transmissivities or VCONT remain in effect unless changed by TVC1 (or some other package). Hydraulic conductivities revert to the values assigned with the Block-Centered flow package.
Explanation of Parameters Read by the TVC1 package.
MXTVC—is the maximum number of hydraulic characteristics that will be changed in any stress period.
IVARAN—is not functional in TVC1. It is reserved as an option to be used in a future version of TVC1 to indicate that the Block-Centered Flow package incorporates variable anisotropy within layers.
ITMP—is a flag and counter
If ITMP > 0, ITMP is the number of hydraulic characteristics that will be read in the current stress period. (Note: horizontal transmisivities and vertical conductances (VCONT) specified in TVC1 will remain in effect for the remainder of the model or until changed).
If ITMP < 0, hydraulic characteristics from previous stress period will be reused.
IRow1—is the row containing the first of the two cells between which the conductance will be altered.
ICol1—is the column containing the first of the two cells between which the conductance will be altered.
ILay1—is the layer containing the first of the two cells between which the conductance will be altered.
Hydch1—is the hydraulic characteristic of the first of the two cells between which the conductance will be altered. If the two cells are on different layers (ILay1¹ ILay2), Hydch1 is the appropriate value of VCONT for that layer. (Only the value for the higher of the two cells is used, a dummy value must be specified for the other cell.) If the two cells are on the same layer, Hydch1 is either the hydraulic conductivity of transmissivity of that layer depending on whether hydraulic conductivity or transmissivity was used in the Block-Centered Flow package. Transmissivity is specified for layers in which LAYCON in the Block-Centered Flow package is 0 or 2. Hydraulic conductivity is specified for layers in which LAYCON in the Block-Centered Flow package is 1 or 3.
IRow2—is the row containing the second of the two cells between which the conductance will be altered.
ICol2—is the column containing the second of the two cells between which the conductance will be altered.
ILay2—is the layer containing the second of the two cells between which the conductance will be altered.
Hydch2—is the hydraulic characteristic of the second of the two cells between which the conductance will be altered. See Hydch1 for more details.
Cheng, X. and Anderson, M.P. 1994, Simulating the influence of lake position on groundwater, Water Resources Research 30(7), 2041-49.
Goode, D.J. and Appel, C.A., 1992, Finite-difference interblock transmissivity for unconfined aquifers and for aquifers having smoothly varying transmissivity: U.S. Geological Survey Water Resources Investigations Report 92-4124, 79 p.
Harbaugh, A.W., and McDonald, M.G. 1996a. User's Documentation for MODFLOW-96, An Update to the U.S. Geological Survey Modular Finite-Difference Ground-Water Flow Model.: U.S. Geological Survey Open-File Report 96-485, 56 p.
Harbaugh, A.W., and McDonald, M.G. 1996b. Programmer's Documentation for MODFLOW-96, An Update to the U.S. Geological Survey Modular Finite-Difference Ground-Water Flow Model.: U.S. Geological Survey Open-File Report 96-486, 220 p.
Hsieh, P.A. and Freckleton, J.R., 1993, Documentation of a computer program to simulate horizontal-flow barriers using the U.S. Geological Survey modular three-dimensional finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 92-477, 32 p.
Kladias, M.P., and Ruskauff, G.J., 1997. Implementing Spatially Variable Anisotropy in MODFLOW. Ground Water, 35(2) 368-370.
McDonald, M.G. and Harbaugh, A.W., 1988, A modular three-dimensional finite-difference ground-water flow model. U.S. Geological Survey Techniques of Water-Resources Investigations Book 6, Chapter A1, 586 p.
McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1992, A method of converting no-flow cells to variable-head cells for the U.S. Geological Survey modular finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 91-536, 99 p.