DESCRIPTION OF THE NUMERICAL MODEL

A three-dimensional numerical model of the hypothetical basin was constructed by selecting boundary conditions and estimating initial values of hydraulic characteristics and recharge. Hydraulic characteristics were adjusted until the model simulated predevelopment conditions within acceptable tolerances. Simulated ground-water levels and discharge to streams and springs agreed with conditions expected in a typical basin.

Approach

Once the conceptual model of the hypothetical basin was defined, the next step was to create a mathematical representation of the basin using a simulation model. To develop the simulation model, first a three-dimensional grid was designed, then it was populated with data on hydraulic characteristics, then boundary conditions were specified, and finally, parameters were hydrogeo.htmladjusted. Model parameters were adjusted until simulated hydrologic conditions were comparable with those that would be expected in the hypothetical basin. The hydrologic conditions considered included the direction and magnitude of hydraulic head gradients (both horizontal and vertical), the rate of seepage to streams, and the rate of discharge to springs. Previous investigators have measured or estimated ranges for these conditions in many parts of the Puget Sound Lowland (as described in previous sections) and these results were relied upon to evaluate how well the numerical model represented the conceptual model. Five parameters were adjusted: horizontal hydraulic conductivity, ratio of horizontal to vertical hydraulic conductivity, streambed hydraulic conductance, hydraulic conductance of springs, and recharge rates. Parameter adjustment continued until the simulated hydrologic conditions were determined to represent typical conditions as defined by previous investigations. Throughout the remainder of the report, the model defined by these parameters and boundary conditions is referred to as the baseline model.

The baseline model was then used to test the effects of one or more discharging wells on ground-water discharge to streams and springs. Several series of simulations were made; each series was designed to show the effect of changes in one variable on ground-water discharge. Variables considered included well depth, distance between well and stream, well discharge rate, and others.

General Features of the Numerical Model

The U.S. Geological Survey's numerical model for simulating ground-water flow, MODFLOW, was used to represent the conceptual model of the hypothetical basin. MODFLOW simulates ground-water flow in three dimensions using finite-difference techniques to solve the partial differential equation describing ground-water movement (McDonald and Harbaugh, 1988). An ancillary program, MODFLOWARC (Orzol and McGrath, 1992), was used to move data directly into MODFLOW from the GIS database.

The MODFLOW program requires that the ground-water system be subdivided, vertically and horizontally, into a finite difference grid of rectangular blocks or cells. The hydraulic properties of the flow system are assumed to be homogeneous within each cell. The saturated flow system of the hypothetical basin was subdivided vertically into 13 layers and horizontally into a 50 column by 70 row grid of 3,500 square cells, each having dimensions of 1,500 ft per side (fig. 6).

Each of the unconsolidated hydrogeologic layers of the conceptual model of the hypothetical basin is represented by an individual layer in the model; bedrock is not explicitly included as a layer in the model, but is represented numerically as a boundary to the system as discussed below. The three glacial sequences, each consisting of four hydrogeologic layers, are represented by the upper 12 model layers, where these layers are present beneath the drift plain. The undifferentiated Qu, are included in the model as the bottom layer (13) beneath the drift plain. The alluvial sediments, Qa, which occur only beneath the major stream valley, were represented in layers 9 through 13 in that part of the model.

Transmissivity, which is the product of the horizontal hydraulic conductivity of the aquifer and its saturated thickness, was specified for each cell. Transmissivity was not changed in the simulations because simulated pumping did not significantly reduce the saturated thickness of aquifers.

Boundaries

Boundary conditions were specified in the numerical model according to the concepts of the flow system used to define the hypothetical basin. An important part of the conceptual model of the basin is that the Tertiary bedrock underlying the basin is a low-permeability unit that does not store or transmit significant quantities of ground water. This was the basis for treating the contact between the Tertiary bedrock and the Quaternary glacial and alluvial sediments as a no-flow boundary in the numerical model. This no-flow boundary extends from where the bedrock is exposed along ridges on the eastern and southern boundaries of the basin, beneath the glacial and alluvial deposits, to the northern and western boundaries of the model (figs. 3 and 4).

The northern and western boundaries of the numerical model coincide with the major stream valleys bounding the drift plain. The streams within these valleys are conceptualized as the discharge areas for local- and intermediate-scale ground-water flow systems in the lowlands; all ground water that enters the system as recharge within the hypothetical basin is assumed to leave the system as discharge within the boundaries of the model. The northern and western boundaries of the model were also specified as no-flow boundaries to reflect the assumption that there is no subsurface ground-water flow out of the basin.

Head-dependent flux boundaries were used to represent ground-water discharge from within the basin by seepage to streams and springs (fig. 6). The underlying equations of these boundary conditions are described by McDonald and Harbaugh (1988).

Hydraulic Characteristics

In order to simulate a ground-water flow system with a numerical model, the hydraulic characteristics of the aquifers and confining beds must be specified for each model cell. The hydraulic characteristics normally required to simulate a ground-water system are thickness, hydraulic conductivity, and specific storage. Specific storage was not required for this analysis because all simulations were for steady-state conditions. The thickness of each hydrogeologic layer, specified on the basis of previous investigations, was variable. Ranges in thickness were 10-50 ft for recessional outwash and till, 10-115 ft for advance outwash, and 10-150 ft for interglacial sediments (table 2).

Horizontal Hydraulic Conductivity
The hydraulic conductivity of deposits in the Puget Sound Lowland spans several orders of magnitude. Even within an individual hydrogeologic layer, hydraulic conductivity values show high variability due to small scale erosional and depositional features. This type of spatial variability in hydraulic conductivity was not represented in the model primarily because the scale of the features that cause heterogeneity in hydraulic conductivity are probably not large enough to significantly affect the simulated response of the system on a basin-wide scale. Also, if this heterogeneity were incorporated in the model, it would have made it difficult or impossible to separate effects of variables such as distance between the well and the stream from effects of spatial variability on hydraulic conductivity on the capture of discharge to streams and springs. Therefore, horizontal hydraulic conductivity was assumed to be uniform throughout the basin for each hydrogeologic layer. The recessional, Qr, and advance, Qa, outwash aquifers were assumed to have the same horizontal hydraulic conductivity, 100 ft/d (table 2). The confining layers (till, Qt, and interglacial sediments, Qf) were assigned horizontal hydraulic conductivities of 0.25 and 1.0 ft/d respectively. The undifferentiated deposits, Qu, are composed mainly of fine-grained sediments, but have coarse-grained interbeds; the relatively high hydraulic conductivity of 25 ft/d assigned to the layer reflected the influence of these coarse-grained beds. Transmissivity varied within model layers due to variations in the thickness of layers. The mean transmissivity of the aquifer layers (Qal, Qr, and Qa) ranged from 2,900 ft2/d for recessional outwash aquifers to 20,000 ft2/d for the alluvial aquifer.

Vertical Hydraulic Conductivity
The heterogeneity of the sedimentary layering in the basin, and particularly of the glacial deposits, imparts considerable anisotropy to their hydraulic conductivity. Anisotropy is the condition of having different properties in different directions. Expressed as the ratio of horizontal to vertical hydraulic conductivities, the anisotropy ratio of unconsolidated sediments is usually greater than one because of the preferred orientation of grains and clasts during deposition. Anisotropy ratios were assigned to each hydrogeologic layer ranging from 10 for the coarse-grained aquifer layers, (Qal, Qr, Qa) to 200 for the mostly fine-grained interglacial layers (Qf). Till layers, Qt, were assigned an intermediate ratio of 100 and the undifferentiated sediments, Qu, were assigned a ratio of 150.

Based on the horizontal hydraulic conductivity value of 0.25 and the anisotropy ratio of 100, the effective vertical hydraulic conductivity of the till was 0.0025 ft/d. This value is within the range of 0.001 to 0.01 ft/d reported by Vaccaro (J.J. Vaccaro, U.S. Geological Survey, written commun., 1993). The vertical hydraulic conductivity of the interglacial deposits was about 0.005 ft/d, or twice that of the till layers.

Recharge

Recharge to the ground-water system was specified at rates computed with the regression equations 1 and 2 described in the previous section on ground-water recharge. Based on average annual precipitation of 44 in/yr, rates of 18 in/yr and 27 in/yr were computed for areas where till and outwash are exposed. The distribution of recharge is shown on figure 6. Mean annual recharge for the basin is about 20 in/yr.

Discharge

Ground-water discharge to streams and springs was represented in the numerical model with the MODFLOW packages RIVER and DRAIN (McDonald and Harbaugh, 1988).

Streams
Streams in the hypothetical basin, both on the upper drift plain and in the major stream valley, were simulated as head-dependent flux boundaries with the RIVER package in MODFLOW. Ground-water flux across head-dependent boundaries is calculated as a piecewise-linear function of the difference between the water levels in the aquifer and at the boundary, and the conductance of that boundary (McDonald and Harbaugh, 1988). Sediments between the stream and the aquifer commonly form distinct bed material differing in hydraulic character from the aquifer sediments themselves. This frequently occurs in low-gradient streams where fine sediment may accumulate to form a bed with low hydraulic conductivity. However, in many streams where gradients are sufficient to maintain high flow velocities or where the streambed is periodically flushed of fine-grained sediments by high flows, there may be little or no contrast between the hydraulic properties of the stream bottom and the aquifer itself.

Stream reaches are defined as the lengths of stream contained within one model cell; stream segments are groups of contiguous reaches that have no tributaries. Streams in the basin are identified by segment numbers, with nine stream segments in the basin: segments 1 through 7 on the upper drift plain and segments 8 and 9 in the lower stream valleys (fig. 6).

Head-dependent boundary conductances are analogous to the vertical conductance between model layers. However, boundary conductances are adjusted external to the model to account for the area of a cell that is covered by the boundary, whereas vertical conductances are adjusted internally. The generalized boundary-conductance equation is

(1) ,

where is the boundary conductance, is the vertical hydraulic conductivity of the sediments between the stream and the center of the model cell containing the stream, is the area of the stream within the model cell, and is the thickness of sediments between the stream and the center of the model cell containing the stream.

Values of boundary conductance for the streams were estimated by assuming that the vertical hydraulic conductivity () of sediments between the stream and the cell is equal to 2 percent of the horizontal hydraulic conductivity of the cell. The average stream area in a cell was estimated by measuring the length of each reach and assuming a width of 10 ft. Although no discrete, low-permeability streambed may be present, a thickness () must be assumed at each model cell to describe the equivalent thickness of the materials that restrict the flow of water between the stream and the ground-water system. This thickness was specified as 2 ft throughout the model. The conductance values resulting from these assumptions are essentially lumped parameters that represent a proportionality constant that determines the flux between the ground-water system and the stream for a given difference between stream stage and ground-water level. The weighted-mean streambed conductance values used in the model ranged from 3,300 ft2/d for segment 5 to 14,600 ft2/d for segment 2; the overall weighted mean conductance was 7,800 ft2/d.

Springs
The springs that occur on the bluffs above the major stream valley and in the canyons of the streams on the upper drift plain were simulated with the DRAIN package of MODFLOW (McDonald and Harbaugh, 1988). The DRAIN package allows simulation of head-dependent boundary flux in a manner similar to the RIVER package. However, the flux is allowed only in one direction; once the head in the aquifer falls below a specified elevation, discharge stops. Springs were specified at 443 cells in layers 3, 5, 7, 9, and 11 where the outwash aquifers are exposed in canyons and on bluffs (fig. 6). The conductances for the springs in the model were computed with the same general equation as for streams, but the assumptions were somewhat different. The area, , of the seepage face at each spring was assumed to span the entire width of the cell (1,500 ft) and the height of the seepage face was assumed to be half of the thickness of the cell. The hydraulic conductivity, , of the spring was estimated to be the same as that of the hydrogeologic layer of the cell. The length of the flowpath, , was the distance from the cell center to the seepage face (750 ft). Using these assumptions, mean spring conductances for spring cells ranged from 1,300 ft2/d for layer 11 to 2,900 ft2/d for layer 3, with an overall mean of 2,400 ft2/d. The discharge altitude for each spring cell was specified as the altitude of the bottom of the aquifer at the point where it is exposed on the bluff or canyon. The 164 spring cells in layer 3 discharged from a mean altitude of 340 ft, and the 25 spring cells in layer 11 discharged from a mean altitude of 66 ft.

Baseline Model Results

About 30 simulations were made in order to obtain a set of hydraulic characteristics and boundary conditions that simulated hydrologic conditions typical of small basins in the Puget Sound Lowland. The process used to obtain these baseline model parameters was trial and error. Following each simulation, simulated hydrologic conditions were compared with the expected conditions and parameters were adjusted to correct errors in the simulated conditions. The final set of model parameters provided the "best fit" between the expected and simulated conditions. Although there are other sets of parameters that could produce equal, and possibly better fits, all of the hydraulic characteristics and boundary conditions in this set are well within ranges expected for the conceptual model of the Puget Sound Lowland.

The baseline model simulates ground-water flow in the hypothetical basin prior to any withdrawal of ground water by wells. The system was assumed to be at equilibrium, or steady state, with discharge balanced by recharge. The expected hydrologic conditions used to calibrate the baseline model are based on data from studies of several areas in the Puget Sound Lowland. Although these data may not reflect steady-state conditions in all cases, the generalized expected conditions are probably typical of the conditions that would be found in an undeveloped basin.

In the following sections, ground-water levels, hydraulic head gradients, and ground-water discharge rates simulated with the baseline model are compared with conditions that would be expected in the typical Puget Sound Lowland basin.

Ground-Water Levels and Hydraulic Gradients
In the baseline model of the hypothetical basin, ground-water flow is generally from the southern, southeastern, and eastern parts of the basin toward the north, northwest, and west. In the shallow part of the flow system, hydraulic gradients and flow directions are strongly controlled by the streams and topography of the upper drift plain (fig. 7). Contours of ground-water altitude show that the shallow aquifers discharge to streams where the contours intersect the streams in V-shapes. Horizontal hydraulic gradients range from less than 15 ft/mi on the northwestern part of the plain to more than 40 ft/mi where they steepen in the uplands on the southeast corner of the plain. Flow directions in the deeper aquifers are more strongly influenced by the discharge areas on the bluffs and in the major stream valleys to the north and west (fig. 8). In the southwest corner of the drift plain, ground water in the shallow aquifer flows from southeast to northwest, whereas deeper ground-water flowpaths are nearly due west toward the springs on the bluff and the major stream. Simulated horizontal hydraulic gradients ranging from about 15 ft/mi to 40 ft/mi, compare closely with typical gradients of 35 ft/mi reported by Vaccaro (J.J. Vaccaro, U.S. Geological Survey, written commun., 1993).

Throughout the drift plain, the simulated vertical component of ground-water flow is downward, except in shallow aquifers near streams where there is upward flow toward the stream. Heads decrease with depth beneath the drift plain, which is a recharge area. Head differences between adjacent aquifers are greatest in the shallow aquifers underlying the plain near the bluff and decrease with depth and to the east and south. Differences in simulated heads are as much as 120 ft between shallow aquifers near the bluff and as little as 20 ft between deep aquifers near the eastern and southern edges of the basin. These simulated differences compare closely with expected differences based on data for southwest King County from Woodward and others (1995), who found differences ranging from about 40 to 150 ft between water levels in the uppermost confined outwash aquifers. Simulated heads increase with depth beneath the major stream valley, which is a discharge area. The maximum simulated vertical head differences between the shallow and deep aquifers beneath the major stream valley were 10 to 20 ft.

Stream and Spring Discharge
Total ground-water recharge to and discharge from the basin, as simulated by the baseline model, was 389 ft3/s (table 3). Seventy-six percent (295 ft3/s) of baseline ground-water discharge was by seepage to streams and 24 percent (94 ft3/s) was by spring discharge. Seventy-three percent (285 ft3/s) discharged to the major stream valley either by seepage to streams (202 ft3/s) or by spring discharge on the bluffs (83 ft3/s). The remaining 27 percent discharged to streams (93 ft3/s) and springs (11 ft3/s) on the drift plain. In the following discussion, springs have been grouped with the stream segment to which they would contribute if the spring discharge flowed to the stream; spring discharge has also been summarized this way in table 3. Total ground-water discharge to streams and springs by stream segment is shown in figure 9a.

Annual baseflow averages about 36 percent of precipitation in Puget Sound Lowland basins (J.J. Vaccaro, U.S. Geological Survey, written commun., 1993). The assumed annual precipitation of 44 in/yr in the hypothetical basin would result in baseflow of about 16 in/yr. Over the 262 mi2 basin this would result in 308 ft3/s of baseflow. The simulated total seepage to streams in the hypothetical basin is only 295 ft3/s (table 3); however, this figure does not include any spring discharge that might flow into the stream and contribute to baseflow. Of the 94 ft3/s of spring discharge, it is likely that at least 13 ft3/s, which is the difference between the simulated and expected baseflows, contributes to streamflow in the basin.

Another assessment of the simulated discharge to streams was made by comparing simulated discharge to streams on the drift plain with estimated baseflow in a small basin in southwest King County. Big Soos Creek (fig. 1) drains an area of 66.7 mi2 and has a mean monthly flow of about 35 ft3/s for October (1967-92), a time when flow in streams in the Puget Sound Lowland are supplied almost entirely by baseflow. This discharge represents the best available estimate of minimum annual baseflow to Big Soos Creek, since baseflow is typically greater during winter and spring months when the water table is higher. Based on this estimate, the minimum baseflow per square mile is 0.52 ft3/s. The drift plain of the hypothetical basin covers about 188 mi2 and therefore the expected minimum baseflow, based on the Big Soos Creek data, is about 100 ft3/s. Simulated discharge to streams on the drift plain was 93 ft3/s by seepage and 11 ft3/s from spring discharge, for a total of 104 ft3/s.

According to Vaccaro (J.J. Vaccaro, U.S. Geological Survey, written commun., 1993), estimated annual baseflow in 12 Puget Sound Lowland basins averages about 0.93 ft3/s/mi2, or nearly twice the minimum baseflow of 0.52 ft3/s/mi2 estimated above. Using the higher value, the expected annual baseflow from the 188 mi2 drift plain in the hypothetical basin would be 175 ft3/s. This value may be a more realistic estimate for comparison with the simulated baseflow, since the model simulates average annual baseflow. On this basis, the simulated baseflow of 104 ft3/s may be somewhat lower than what would be expected from an area this size. There are several sources of uncertainty in the estimate of annual baseflow. Considering only the uncertainty in determining the area contributing to baseflow, the uncertainty in the estimate of annual baseflow is probably percent or more. Because simulated baseflow from the drift plain (104 ft3/s) falls within the range of estimated baseflow (100 to 175 ft3/s), the distribution of ground-water discharge in the model is a reasonable representation of a typical basin.

The rate of ground-water discharge per mile of stream, or specific discharge, was computed for each stream segment (fig. 9b). Including streams and springs, stream segments 1 through 6 had simulated specific discharges ranging from 1.68 to 3.92 ft3/s/mi, while segments 7 through 9 had discharges ranging from 5.25 to 7.79 ft3/s/mi (fig. 9b). The large difference is due to the positions of segments 7 through 9 in the lower part of the flow system where they receive discharge from a broader recharge area than do the streams on the drift plain. If spring discharge is not included, specific discharge to streams on the drift plain ranges from 0.8 to 3.9 ft3/s/mi (segments 1 through 6 only); this range is similar to the range of 0.3 to 3 ft3/s/mi determined for reaches of Big Soos Creek in southwest King County (Woodward and others, 1995).

The specific discharge of simulated springs was computed by dividing the total discharge of springs contributing to each stream segment by the length of spring discharge area adjacent to that stream segment. Specific discharges ranged from about 0.3 to 1 ft3/s/mi. The values at the lower end of this range are comparable to values at the upper end of the range of 0.01 to 0.27 ft3/s/mi reported by Woodward and others (1995) for southwest King County. As many as four outwash aquifers are exposed on the bluffs above the major stream valley, and of the 453 spring cells in the model, 78 percent are located on these bluffs. Consequently, some of the highest specific discharges for springs (0.7 and
1.0 ft3/s/mi) are on segments 8 and 9 in the major stream valley below the bluffs. Discharge per spring, however, is not significantly higher; the springs on the bluffs comprise 78 percent of all springs in the model and account for 88 percent of the total spring discharge, whereas springs on the drift plain comprise 22 percent of all springs and discharge 12 percent of total spring discharge (table 3).