1、GEOLOGY 724Final ProjectModel Calibration and PredictionFall 2008NOTE: (1) You will need to run the final project using GWV, version 5.16 Build 6 since the GWV basecase file you will be given was created with that version of GWV. (2) Use MODFLOW88/96.-The final project is an exercise in model calibr
2、ation and prediction. The problem is very loosely based on a field example from the Basin and Range Province in the Western U.S. documented in USGS Professional Paper 1409-E. The project has four steps: model calibration, drawdown prediction, particle tracking, and presentation of results.Your first
3、 task is to calibrate a steady-state model. You will work in groups; each group is asked to work independently. When we compare results we will see that each group will have a different, though equally valid, calibrated model and thereby we will demonstrate that the calibration process yields non-un
4、ique solutions.After calibration, the next step is to introduce pumping wells into the calibrated model and predict heads and drawdowns. Then you will use the particle tracking code, ModPath to check for possible contamination of the wells. When we compare results, we will see that results of the pr
5、ediction phase of the modeling are different for each group. You will be given a gwv file for input to GwVistas. This file is set up for an initial calibration run. You will need to change parameter values in order to calibrate the model. During calibration you will attempt to match a set of “field”
6、 measured head values and fluxes, known as calibration targets. Parameters that are adjusted during calibration are called calibration parameters.I. CALIBRATIONConceptual Model. Rabbit Brush Valley is a closed basin located in an arid region in the western U.S. The basin fill consists of coarse- to
7、fine-grained unconsolidated sediments. The aquifer is recharged by runoff from the mountains surrounding the basin. Recharge occurs only at the foot of the mountains and is concentrated in canyons. Water is discharged by evapotranspiration (ET) that occurs within and around the playa in the interior
8、 of the basin (Fig. 1).Model Design. The model design of 3 layers, 11 columns and 16 rows is shown in Figure 1. Horizontal grid dimensions are 5280 ft (1 mile) by 10560 ft (2 miles). The boundaries of the basin are simulated using no-flow boundary conditions. The recharge option (under properties in
9、 GwVistas) is used to enter recharge rates at the edge of the basin. The evaporation option (under properties) is used to simulate ET using head-dependent conditions in the playa area.We will use a three layer model in order to simulate vertical changes in lithology of the sediments and vertical gra
10、dients within the flow system. Layer types are defined in the BCF package (under MODFLOW). The first layer is unconfined (laycon=1) and includes the playa deposits. The second and third layers are simulated as “unconfined- transmissivity varies” (laycon=3). When 2laycon=3, MODFLOW will calculate lay
11、er transmissivities given the top and bottom layer elevations and the hydraulic conductivities for the layer. The third layer has a variable bottom, which is entered into the bottom option (under properties). The elevation of the bottom of the first layer is equal to zero (Figure 1). (The elevation
12、is actually 6000 ft above sea level, but for convenience in viewing the head array we will define a new datum.) The top of the second layer is also zero. The bottom of the second layer is -250 ft and the top of the third layer is -300. Note that there is a 50 ft gap between the layers. In fact, ther
13、e is a 50 ft confining bed (aquitard) between layers two and three, which is not shown in the figure and is not simulated directly by the model. Because there is a gap between the top of layer 3 and the bottom of layer 2, GwVistas knows that the gap is occupied by a confining bed (aquitard). The con
14、fining bed will be represented by entering a value for the vertical hydraulic conductivity of the confining bed in the leakance option (under properties).Solution Package. We will use the SIP package. HCLOSE is set to 1E-5 and 200 iterations are allowed to reach closure. Model Input. The time unit f
15、or this model is years (under MODFLOWBasic package). Therefore, hydraulic conductivities, recharge rates, ET rates, pumping rates, and fluxes are all entered/computed with years as the time unit.Hydraulic Conductivity: The lithology of the aquifer is entered via hydraulic conductivities in the prope
16、rties menu. From field information we know that layer 1 is heterogeneous and consists of fine to coarse sand with some gravel. Pumping tests at wells SW1 and NP2a gave K= 20,000 ft/year. Sediments under the playa are fine grained and characterized as “silt and clay”. Layers 2 and 3 are thought to ha
17、ve lower conductivity than layer 1. The vertical hydraulic conductivity of the aquitard between layers 2 and 3 is entered via the leakance option (under properties). Calculation of leakance is discussed below.To achieve calibration you might choose to introduce anisotropy and heterogeneity within th
18、e lithologies defined in the initial model. Hence, hydraulic conductivity values for all layers as well as the Kz of the aquitard should be treated as calibration parameters. That is, you should change, within reason, the parameters values already entered into the K database in order to achieve cali
19、bration. You may introduce anisotropy into the model by setting values of both Ky and Kz to be different from Kx for all units. You may assume that the ratio of horizontal to vertical conductivity for all lithologies falls between 1 and 1000. Leakance: Leakance (under properties) describes the verti
20、cal transmission properties between layers (A in the playa fringe the rate is 1.5 ft/day. You may accept these values as knowns; they are not calibration parameters. The extinction depth, which is also set under the ET option is a calibration parameter. It is initially set to 6 ft in both the playa
21、and the playa fringe. You should assume that the extinction depth is constant within each zone but that the depth may be greater in the playa fringe than in the playa, owing to the presence of deep-rooted phreatophytes in the playa fringe. According to the MODFLOW manual the extinction depth is “fre
22、quently assumed to be on the order of six to eight feet below land surface (although considerable 4variation can be introduced by climatic factors, the presence of deep-rooted phreatophytes or so on)“. In the arid setting of Rabbit Brush Canyon, the extinction depth could be greater than eight feet
23、in both ET zones. Initial heads. In the gwv file you will download from the website, the starting heads are set to 60 m everywhere in the model. With the initial parameter set, the model will converge with this choice of initial values within the maximum 200 iterations specified in the SIP input. Co
24、nvergence is sensitive to the initial head values. So, you should run the model with the initial parameter set to generate initial heads to use during your calibration runs. For your first (and subsequent) calibration runs, go to MODFLOWinitial heads and select “Read heads directly from binary head-
25、save file”. Then use the browse button to select the head file generated using the initial parameter set.A problem you may encounter when running the model is that ET may cause some of the cells in the 1st layer to go dry. That is, withdrawal of water by ET might cause head in the top layer to drop
26、below the bottom of the layer during iteration. Ideally, we would like the iterations to continue and allow the model to bring the head back up to layer 1 in the next iteration. But MODFLOW has the feature that when the head drops below the bottom of the layer, the cell dries up and is taken out of
27、the simulation. GWV will color in these nodes (dry cells) with dark purple. You dont want this to happen during your simulation. If you see dark purple nodes, you should adjust your parameter values to eliminate the dry cells in the next trial simulation.Calibration ProcessCalibration Parameters. Th
28、e calibration parameters are:1. Hydraulic conductivity values in all layers, including Kz values. (The vertical anisotropy ratio (Kx/Kz) for all lithologies in layers 1, 2 3. Recharge rates (recharge must total 7.1369 E08 ft3/yr);4. Extinction depths for two ET zones (equal to or greater than 6 ft)5
29、Calibration Targets. Water-level measurements in four shallow wells (SW), one deep test well (DW), and two nested piezometers (NP1 and NP2) form the calibration values for heads. Discharge into the playa and playa fringe via evapotranspiration was estimated from lysimeter and vegetation studies and
30、these measurements provide calibration values for flux. Calibrate your model to the calibration targets listed below.Name Node Head ET ratelayer row column (ft) (ft3 /yr)SW1 1 4 6 103.61 SW2 1 8 3 54.66 SW3 1 9 11 37.00SW4 1 10 6 35.99 NP1a 1 8 6 33.15 NP1b 2 8 6 38.05 NP1c 3 8 6 40.51 NP2a 1 14 4 8
31、1.87NP2b 2 14 4 81.83NP2c 3 14 4 80.51DW 3 10 8 35.92ET1 1 7 7 - 1.72 E7ET2 1 8 10 - 8.60 E6ET3 1 9 8 - 2.67 E6ET4 1 10 6 - 1.65 E7ET5 1 12 8 - 3.61 E7Calibration with GwVistasCalibration consists of adjusting the calibration parameters listed above within reasonable limits until you achieve a satis
32、factory match to the calibration targets. You will find that changing some of these parameters has little effect on the solution. You should do the calibration by trial-and-error adjustment of the parameters in multiple runs of the model. It is suggested that you develop a procedure to keep track of
33、 which parameters you have changed and the effect of each change on the calibration. In this way, you are doing a kind of sensitivity analysis as you perform the calibration. You can try to use the automated calibration options in GwVistas (e.g., PEST) but these codes can be tricky to manipulate.The
34、 calibration targets for head have already been entered into the gwv file that you will use for your initial calibration run. When you import results, check to make sure the box labelled “interpolate targets-observation data” is not checked. You will find this box at the lower left hand side of the
35、import results window. When the box is unchecked, GWV will compare the head at the center of the cell with the target value.6When you have a trial solution, go to Plotcalibration and look at the calibration statistics for the heads. Note that for the initial calibration run the absolute residual mea
36、n (ARM) is 17.89 ft and the sum of squares is 6630 ft2. During subsequent calibration runs, you should reduce these numbers and try to get them as close as possible to zero since for this exercise we are assuming that there is no error in the target values. Remember that you are unlikely to achieve
37、a perfect match to all the calibration values. Dont worry too much about getting a perfect calibration but try to get an ARM value of around 1 or lower if at all possible. Calibration will require luck, given the time limits for this project. You can produce a scatter plot of observed vs. calculated
38、 heads, by selecting the appropriate plot option under Plotcalibration. Although GwVistas will consider certain types of flux values for calibration targets, it does not yet support ET calibration targets. (In GwVistas, flux targets must be considered boundary fluxes, under the BC menu. In GwVistas,
39、 however, ET is considered a property and is listed under the properties menu.) Therefore, you will have to check the ET calibration targets manually. To check the ET rates, use the window option under Plotmass balance. Move the window to cover just one cell in the playa/playa fringe area for which
40、you have a calibration target. When the mass balance screen appears, check the numbers at the top of the screen to make sure that you have defined the window for only one cell. The ET discharge rate in ft3/yr is found in the outflow column under ET. Calculate the ARM error for the initial model and
41、attempt to reduce the error to zero in subsequent calibration runs.II. DRAWDOWN PREDICTION Suppose we try to capture some of the water “lost“ by ET by installing some pumping wells.Set up a new directory for your pumping simulation and copy your calibrated gwv file into the new directory.Simulate th
42、e effects of steady-state pumping from clusters of pumping wells located in the following five cells:Name ROW (I) COLU MN (J) LAYER (K)PW1 5 5 1PW2 8 6 1PW3 10 4 1PW4 13 5 1PW5 13 8 1The total pumping rate for each cell is -0.99E08 ft3/year. With all 5 pumping centers, we capture around 40% of the f
43、low that otherwise would discharge as ET to the playa. (1) Each pumping center is input using the well option (under BC). (2) Before you run MODFLOW, check the MODFLOWpackages screen. Be sure the well option is checked and assigned an IUNIT number (e.g., 12); also assign an IUNIT number for the cell
44、-by-cell flow for the Well Package. (3) You should use your steady-state calibrated heads as the initial heads (under MODFLOWinitial 7headsRead heads directly from binary head-save file) so that GV will calculate meaningful drawdown values. Use the browse button to specify the location of your calib
45、rated head file.III. PARTICLE TRACKINGSuppose there is a possibility that water recharging the aquifer in the canyons is contaminated and could contaminate the well water. We will use the program MODPATH, which is one of the codes supported by GwVistas, to insert imaginary “particles“ in the cells b
46、elow the canyons and trace flow paths. These are particles 1 - 6 below. Also suppose that there is a proposal to bury high-level radioactive waste in layer 3; the waste is represented by particle 7 (see below). We want to determine which, if any, of these “particles” will reach any of the wells and
47、which will be discharged into the playa. The z offset option governs the vertical placement of the particle. A value of 1 means the particle is placed at the top of the cell. A value of 0 means the particle is placed at the bottom of the cell.Particle # ROW COLUMN LAYER Z-Offset _1 1 5 1 12 5 1 1 13 10 11 1 14 14 2 1 15 12 3 1 16 4 7 1 17 8 3 3 0MODPATH is a particle tracking code that use
