The three-dimensional flow model was developed and calibrated to be consistent with hydraulic property distributions (i.e., transmissivities) and boundary conditions used in calibrating the two-dimensional inverse model. The overall transmissivity of the aquifer, determined by inverse modeling with the two-dimensional model, was used with available data on hydraulic properties of each of the sub-water-table layers of the model to redistribute the calibrated transmissivity over the model layers used in the three-dimensional model. This aspect of the development of the three-dimensional model is described in Wurstner et al. (1995).
In this section, the development of the three-dimensional flow model is described, including:
The two-dimensional groundwater flow model, as described in Wurstner et al. (1995), was calibrated using a statistical inverse method as implemented in the CFEST Inverse (CFEST-INV) code. This code was designed to be compatible with CFEST and includes algorithms to calibrate the flow model to prior (measured) hydrogeologic parameters such as head, boundary conditions, and hydraulic conductivity or transmissivity. The code can be applied to calibrate a single-layer areal (x-y) model or a multilayered vertical (x-z or x-y-z) model. Insufficient data were available to calibrate the three-dimensional model directly, so the two-dimensional model was calibrated first. The resulting distribution of transmissivity was transferred to the three-dimensional model based on knowledge of subsurface hydrogeology. The procedures for calibrating the two-dimensional model and transferring the results to the three-dimensional model are described in Wurstner et al. (1995).
In the previous application of CFEST-INV to the Hanford Site's unconfined aquifer data set, the code was applied so that only the heads were constrained. The transmissivities were not constrained, and as a result, reached unrealistic values (greater than 10,000,000 m2/d) near Gable Mountain Pond. Consequently, the calibration of the two-dimensional model was reevaluated in FY 1997. The same 1979 time period, hydraulic head data, finite-element grid, boundary conditions, and artificial and natural recharges reported in Wurstner et al. (1995) were used in the recalibration effort. As previously described (see Section 3.1.1), Cold Creek Valley was treated as a constant-head boundary (148 m) for model calibration. Dry Creek Valley was treated as a constant-flux boundary (1,333 m3/d). In refining the inverse calibration, the initial transmissivities were adjusted by hand near Gable Mountain and B Ponds to be closer to the final calibrated distribution. A few of the transmissivity zones were adjusted to better reflect the trends in transmissivity developed in previous calibration efforts by Jacobson and Freshley (1990) and Cearlock et al. (1975).
The inverse calibration was performed in several steps to constrain the solution. The starting conditions were adjusted to be much closer to the final solution, and some transmissivity zones were modified to better reflect the distribution of transmissivity in the unconfined aquifer. The resulting transmissivity distribution is illustrated in Figure 4.1, and the hydraulic head solution predicted with the calibrated model is compared with the interpreted water table for 1979 in Figure 4.2. The inverse transmissivities vary from less than 250 to 125,000 m2/d. The average value is 6,780 m2/d. This is consistent with the range of transmissivities from the inverse calibration by Jacobson and Freshley (1990); the maximum transmissivity for their calibration that included areal recharge was 112,000 m2/d.
The hydraulic head distribution predicted with the calibrated model demonstrates good agreement with observed (measured) water levels over much of the site. Areas where agreement is poor include the gap between Gable Butte and Gable Mountain and between the 200-West Area and the boundary at Cold Creek Valley. In the latter region, the interpreted gradient is steeper than in the calibrated model.
The three-dimensional flow model was developed and calibrated based on the two-dimensional inverse transmissivities. Data on the hydrogeology of the individual layers in the model were preserved in the calibration. As in Wurstner et al. (1995), the database of geologic and hydrologic information was developed to be independent of the model grid. Thus, modifications to the finite-element grid (see Figure 3.2) can be made and the data resampled easily. This feature facilitated construction of the three-dimensional flow model.
To develop the regional numerical model of the Hanford Site, the Geological Finite-Element Synthesis Tool (GEOFEST) described in Foley et al. (1995) was used to transfer the interpreted water table and the extent and thickness of major hydrogeologic layers as input to CFEST. Bottom elevations of layers, elevation of the water table, and grids representing hydraulic conductivity zones were put into GEOFEST. The hydrogeologic structure represented in the numerical model is described in Wurstner et al. (1995).
The three-dimensional numerical model described in Wurstner et al. (1995) includes the nine hydrogeologic units above the basalt. The basalt layer (unit 10) described in the conceptual model of Wurstner et al. (1995) was excluded from the numerical model over most of the model area to reduce the number of nodes and, consequently, the model run times. The exception was for an area in the vicinity of the basalt outcrops at Gable Mountain and Gable Butte. A limited amount of basalt was retained in these areas so that the model could resolve the location of the intersection of the water table with basalt subcrops as the water table drops in response to the eventual cessation of Hanford discharges to the aquifer.
To further minimize the number of nodes and elements and model run times, two additional refinements and modifications were made:
Hydraulic conductivities were assigned to model layers based on data from aquifer pumping tests, slug tests, and some laboratory tests as described in Wurstner et al. (1995). Because data for each of the units are sparse, the general distribution of the inverse calibrated transmissivities was preserved during the three-dimensional model calibration.
The transmissivity distribution derived from the two-dimensional inverse calibration was preserved using the following steps:
Step 1. The hydrogeologic structure of the nine major units of the unconfined aquifer was constructed using the calibrated water-table surface predicted with the inverse calibration using the GEOFEST software processing utility.
Step 2. The transmissivity distribution derived from the inverse calibration was distributed among the major conductive hydrogeologic units (i.e., units 1, 3, 5, 7, and 9) at every node location using the ratios of average values of the major units given in Table 4.1. The hydraulic conductivity of the less conductive mud units (i.e., units 2, 4, 6, and 8) were assigned the average property values shown in Table 4.1. The transmissivity of these low-permeability units was then calculated using the calculated saturated thicknesses and these average hydraulic conductivity values. The transmissivities of these units were subtracted from the inverse, calibrated transmissivity before it was redistributed among the other water-bearing units (units 1, 3, 5, 7, and 9).
Table 4.1. Assumed Hydraulic Conductivities for Major Hydrogeologic Units Used in the Redistribution of Inverse Calibrated Transmissivities to the Three-Dimensional Model
| Hydrogeologic Unit | Estimated Range of Saturated Hydraulic Conductivities (m/d) | Reference(s) |
| Hanford Formation (Unit 1) | 1 to 1,000,000 | Wurstner et al. (1995): Thorne and Newcomer (1992) |
| Palouse Soils (Unit 2) | 0.057 to 0.12* | Khaleel and Freeman (1995) |
| Plio-Pleistocene Unit (Unit 3) | 0.005 to 17.3* | Khaleel and Freeman(1995) |
| Upper Ringold Unit (Unit 4) | 0.0003 to 0.09 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
| Middle Ringold - Unit E (Unit 5) | 0.1 to 200 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
| Middle Ringold - Unit C (Unit 6) | 0.0003 to 0.09 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
| Middle Ringold - Unit B and D (Unit 7) | 0.1 to 200 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
| Lower Ringold Unit (Unit 8) | 0.0003 to 0.09 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
| Basal Ringold Unit (Unit 9) | 0.1 to 200 | Wurstner et al. (1995); Thorne and Newcomer (1992) |
Step 3. The redistributed transmissivities were then used to develop areal distributions of hydraulic conductivities for each major unit, taking into consideration the calculated saturated thickness generated by GEOFEST from the hydrogeologic structure information.
Step 4. The resulting distributions of hydraulic conductivity were combined with the generated hydrogeologic structure of all major units and used as input to a steady-state run with CFEST using the confined aquifer option. Constant-head boundaries used in the inverse calibration were used to hold boundary nodes at and below the water-table surface in the same locations used in the inverse calibration ( i.e., Cold Creek Valley, North and East Dry Creek Valleys, and along the Rattlesnake Hills).
Step 5. The head distribution derived from the previous step was then used to repeat Steps 1 through 4 to arrive at a new calibrated head surface that would reflect the structure and distribution of hydraulic conductivity derived from the original inverse calibrated transmissivity distribution.
Step 6. The model run made in Step 5 was then replicated with CFEST, using the unconfined aquifer option and with the calculated constant fluxes derived at constant-head boundary nodes substituted at constant head nodes. Constant fluxes were the only input to the most significant water-bearing units at the specified boundary-node locations.
Step 7. As a final check, the model run made in Step 6 was then replicated with CFEST, using the unconfined aquifer option and the calculated constant fluxes derived at constant-head boundary nodes substituted at constant head nodes.
Using the steps given above, the final water table and hydraulic conductivity layering combination for the three-dimensional model was determined. Figure 4.4 shows the hydraulic conductivities used for the uppermost hydrogeologic units represented by surface elements in the model. Plots of the steady-state solution from the three-dimensional model for 1979 with the observed 1979 water table, shown in Figure 4.5, provide a comparison of model results with conditions interpreted for 1979. This figure and Figure 4.6, which provides a shaded map of the differences between these two surfaces, shows that in general, predicted water levels throughout most of the modeled region are within a meter or two of observed conditions. A statistical comparison of the difference between the predicted water table and the interpreted water-table surface, summarized in Table 4.2, provides additional information on the goodness of fit at all 1457 surface node locations.
Another measure of goodness of fit is a comparison of predicted water-table elevations with those measured in individual wells. This information is summarized in Figure 4.7, which provides a comparison of predicted versus observed water levels measured at 100 wells. Of the 100 wells plotted in this figure, predicted water levels at comparable locations were within a meter of observed water levels at 85 well locations and well within 5 m of observed values for all well locations.
Once the three-dimensional model was calibrated, the constant-head boundary conditions used to simulate fluxes for Cold Creek and Dry Creek Valleys and off the Rattlesnake Hills were changed to constant fluxes. Fluxes were distributed over all vertical nodes representing the boundary. The constant flux boundary conditions calibrated from the 1979 steady-state prediction were 2881 m3/d for Cold Creek Valley, 1207 m3/d for Dry Creek Valley North, 328 m3/d for Dry Creek Valley East, and 3104 m3/d for the Rattlesnake Hills area. The previous estimate of the Cold Creek boundary condition from the two-dimensional inverse calibration by Jacobson and Freshley (1990) was 8805 m3/d, and the estimate made with VAM3D by Law et al. (1996) was 10368 m3/d. The estimate of this same flux for the Variable Thickness Transient (VTT) flow model by Kipp et al. (1972) was 9104 m3/d. The current estimate for Cold Creek (2881 m3/d) is lower than previous estimates.
The flux at the Rattlesnake Hills boundary condition previously estimated with the two-dimensional inverse calibration by Jacobson and Freshley (1990) was 1331 m3/d, and the estimate made with VTT by Kipp et al. (1972) was 3105 m3/d. The flux used in the model by Jacobson and Freshley (1990) accounted for the boundary influx at the east valley of Dry Creek. In the model by Law et al. (1996), the flux along Rattlesnake Hills was assumed to be negligible. The current estimate for Rattlesnake Hills (3104 m3/d) is lower than previous estimates.
A previous estimate for the Dry Creek Valley boundary from the model by Law et al. (1996) in which a boundary flux used was about 44,064 m3/d. The current estimate for both the north and east valleys of Dry Creek (1535 m3/d) is significantly lower than this previous estimate.
In the three-dimensional model, the prescribed head for the Columbia River boundary was held constant at nodes located at the left edge of the river shore (looking upstream) and the node located in the middle of the river channel. Nodes located vertically below the node positioned in the middle of the river were assigned as no-flow boundaries. This was done to approximate the broad boundary of lateral and vertical discharge from the unconfined aquifer through the bed of the river, which is assumed to be a point of discharge for water and contaminants flowing toward the river from both the Hanford Site and Franklin County shorelines.
Table 4.2. Statistical Comparison of Differences Between Observed
and Predicted 1979 Water-Table Elevations
| Statistical Categories | Values, in meters |
| All Nodes (1457 nodes) | |
| Maximum Error | 5.13 |
| Minimum Error | -6.73 |
| Average Error | -0.63 |
| Absolute Value of Error | 1.08 |
| RMS of Error | 1.64 |
| Nodes with Positive Errors (509 nodes) | |
| Average Error | 0.65 |
| RMS of Error | 0.93 |
| Nodes with Negative Errors (948 nodes) | |
| Average Error | -1.31 |
| RMS of Error | 1.91 |
The distribution of natural recharge estimated by Fayer and Walters (1995) for 1979 (see Figure 3.1) was used as input for both the steady-state and transient-flow models. This distribution of natural recharge was assumed to remain constant during the entire simulation period. In reality, fires, land use, and natural progression of the ecosystem will change the distribution of vegetation types and, thus, future recharge conditions.
This section of the report describes the overall approach and results for 1) the transient calibration of the three-dimensional model from 1979 to 1996 conditions, and 2) the application of the three-dimensional model to simulating transient-flow conditions in the unconfined aquifer from 1996 to the year 4000. The simulation of transient-flow conditions provided the hydraulic framework and conditions used to make long-term forecasts of contaminant plumes currently being transported in the unconfined aquifer system.
After generating a reasonable steady-state solution for the three-dimensional model, the response of the model to transient flux input (e.g., B Pond waste discharges) was calibrated. Transient simulations with the two- and three-dimensional flow models required data on artificial recharge. Artificial recharge was applied at the nearest node to the actual discharge site. The artificial recharge values used in the transient-flow models and the Rattlesnake Mountain spring discharges are summarized in the Appendix (Table A.1). The yearly totals of all artificial recharge to the 200 Areas as a function of time are shown in Figure 4.8. Beyond 1995, the only facilities that are projected to discharge waste water to the unconfined aquifer are W-049H (the 200 Areas Treated Effluent Disposal Facility [TEDF]), SALDS, and the Richland North Area well field. The City of Richland recharges the groundwater with Columbia River water at the Richland North Area well field so that groundwater is available for withdrawal when the river pumping station is shut down for maintenance or additional water is needed to meet irrigation demands. As seen in Figure 4.8, discharges to the TEDF cease after the year 2026.
A number of calibration simulations were performed to evaluate the appropriate specific yield value to achieve the best overall match. The storage properties of the aquifer (specific yield) control how the model responds to changes in flux. Wurstner et al. (1995) estimated the range for specific yield of Unit 1 to be from approximately 0.1 to 0.3.
Three-dimensional transient simulations were conducted for the time frame from 1980 to 1995 using specific yields ranging from 0.1 to 0.35. The simulations used semi-annual time steps, semi-annual water-table updates, and yearly averaged fluxes. The best fit to the observed data was achieved when a specific yield of 0.1 was used where the Ringold Formation is found, and a specific yield of 0.25 was used where the Hanford Formation is found. The comparison of predicted with observed conditions for 1996 using this combination of specific yields is shown in Figure 4.9. Figures 4.10 4.11, 4.12, 4.13, and 4.14 show hydrographs comparing the model results with observed heads from 1979 to 1996 for 20 groups of wells. Figure 4.15 shows the locations of the wells used in the hydrographs.
Past analyses of post-Hanford unconfined aquifer conditions have considered land uses that could significantly alter the long-term behavior of the unconfined aquifer beneath the Hanford Site. An example of such use is large-scale irrigation of the Site, postulated for certain agricultural purposes. The potential for large-scale agricultural irrigation occurring on the Hanford Site in the future was also examined. Based on consultations with staff from the Agricultural Research Service at the Agricultural Experiment Station in Prosser, Washington, large-scale irrigation on the Hanford Site is considered to be unlikely for a variety of factors, including
Because of these conclusions regarding potential irrigation on the Hanford Site, projections of post-Hanford water-table conditions in this study focused mainly on predicting the impact of Hanford operations ceasing and the resulting changes in artificial discharges that have been used extensively as a part of site waste-management practices. Simulations of transient-flow conditions from 1996 through the year 4000 were conducted with the three-dimensional model using the distribution of hydraulic conductivities from the steady-state calibration and the distribution of specific yields developed from the transient calibration (0.25 for Hanford Formation and 0.1 for the Ringold Formation). The same natural recharge distribution that was used in the steady-state calibration was used in these future conditions simulations. A comparison of the initial conditions for this simulation with observed values was provided earlier in Figure 4.5. The water tables estimated with the three-dimensional model for the years 2000, 2100, 2200, and 2350 (Figures 4.16, 4.17, 4.18, and 4.19) show an overall decline in the hydraulic head and hydraulic gradient across the entire water table within the modeled region. From the selected hydrographs shown in Figures 4.20, 4.21, 4.22, 4.23, and 4.24, it can be observed that the different areas approach steady state with different rates of time. The areas along the Columbia River north of the gap between Gable Butte and Gable Mountain (areas 17 and 19 in Figure 4.24) have the shortest response times, and water levels in this region reach steady state by the year 2100. The area between Gable Butte and Gable Mountain (areas 18 and 20 in Figure 4.24) reach steady-state conditions sometime between the years 2200 and 2300. Area 3 (Figure 4.20) has the slowest response time, and water levels in this area reach steady-state conditions only after 2400. The remaining areas plotted in Figures 4.20, 4.21, 4.22, 4.23, and 4.24, including the area south of Gable Mountain and east of the 200-West Area plateau, all reach steady-state conditions between the years 2300 and 2350.
Flow modeling results also suggest that as water levels drop in the vicinity of central areas in the model where the basalt crops out above the water table, the saturated thickness of the unconfined aquifer greatly decreases and the aquifer may actually dry out. This thinning/drying of the aquifer is predicted to occur in the area between Gable Butte and the outcrop south of Gable Mountain, resulting in the northern area of the unconfined aquifer becoming hydrologically separated from the area south of Gable Mountain and Gable Butte. Therefore, flow from the 200-West Area and the northern half of the 200-East Area, which currently migrates through the gap between Gable Butte and Gable Mountain, will be effectively cut off in the next 200 to 300 years. In time, the overall water table (including groundwater mounds near the 200-East and -West Areas) will decline, and groundwater movement from the 200-Area plateau will shift to a dominantly west-to-easterly pattern of flow toward points of discharge along the Columbia River between the Old Hanford townsite and the Washington Public Power Supply System facility.
The decline in the water table from current conditions (January, 1997) to those predicted in the year 2350 is shown in Figure 4.25. In the 200-East Area, the maximum change in head is realized in the area of a current groundwater mound where the water table drops between 4 and 6 m. In the 200-West Area, model results indicate a maximum drop in the water table of about 11 m during this same period.
As a qualitative measure of the appropriateness of the predicted water table estimated with the three-dimensional model, the predicted water table for post-Hanford conditions for the year 2350 was compared to estimated water-table conditions in 1944 provided in Kipp and Mudd (1974) and shown in Figure 4.26. In general, the water table predicted for post-Hanford conditions agreed favorably with that estimated for 1944 conditions for the central part of the Hanford Site and the area north of Gable Mountain and Gable Butte. However, model results indicated a much different water table than was estimated for pre-Hanford conditions. This is apparent in two areas: 1) the area west of the 200-Area plateau, where higher predicted hydraulic heads reflect boundary conditions that consider the effect of irrigation from areas outside of the modeled region, and 2) the area north of Richland, where the model simulates the hydraulic effect of the North Richland well field. Discrepancies in these two areas of the water-table surface may be attributable to one or more of the following factors:
