Hydrogeology

The glaciated northeastern United States in areas underlain by crystalline rock is characterized by generally shallow depths to water (less than 30 ft bls) because the bulk permeability of the bedrock is low and recharge is relatively high (Harte and Winter, 1995). Bedrock exposures are few. A mantle of glacial till or recent alluvium overlying the bedrock in the uplands and stratified-drift deposits commonly overly the bedrock in the valleys. Flow is primarily through fractures that were created by a combination of cooling, uplift, and structural deformity. Nonfractured rock blocks have low (less than 5 percent) primary porosity (Zhou and others, 2007). Flow is often proportionally greater at shallower depths (less than 300 ft bls) than flow at deeper depths in the bedrock; although highly transmissive fractures at deep depths have been intercepted by wells. Nevertheless, the transmissivity of the bedrock is typically several orders of magnitude less than the overlying unconsolidated sediments. Transmissivities of crystalline-rock aquifers are typically less than 1,000 feet squared per day (Tiedeman and others, 1997), whereas transmissivities between 2,500 and 50,000 feet squared per day have been measured for the stratified-drift aquifers of the northeast (Randall, 2001).

Surface water and groundwater interactions are an important part of the groundwater system in the glaciated northeastern United States. Surface water recharges the groundwater where upland streams cross permeable deposits in the valleys. Surface water also can recharge the groundwater where there are large changes in the bulk permeability of an aquifer in valleys (Harte and Kiah, 2009). Groundwater discharges to surface water in uplands and valleys, but particularly in valleys.

In rural areas of the northeast, residential water-supply wells and monitoring wells open to the crystalline-rock aquifer intersect multiple fractures over the open interval of the well (open borehole). Groundwater, therefore, mixes in the borehole from the contribution of groundwater flow from multiple fractures (Harte, 2013).

In the model area, most of the uplands above altitudes of more than 330 ft mean sea level (msl) surrounding OU1 are served by residential wells set in the crystalline-rock aquifer (fig. 2). Most residential wells outside the MSGD aquifer are in the uplands. Yields of residential water wells range from 0.13 to 150 gallons per minute (gal/min) with an average of approximately 13 gal/min (New Hampshire Department of Environmental Services, written commun., 2011). The median well depth for residential water wells is 400 ft bls (an average of 421 ft bls). Most monitoring wells in OU1 are drilled to similar depths.

The crystalline-rock aquifer that underlies the model area consists of tonalite, granite, and gneiss (figs. 1 and 2; Bennett and others, 2006). Some metapelite strings are present in the uplands. Burton and Harte (2013) mapped the metapelite strings in the uplands coincident with a schist unit. Tonalite is mapped to the west and lies within the higher terrain. Tonalite is separated from the granite and gneiss to the east in the Milford-Souhegan River Valley (Bennett and others, 2006; Burton and Harte, 2013) by the Campbell Hill Fault (figs. 1 and 2). The fault appears to be primarily a normal fault that also may have a strike slip component (Burton and Harte, 2013). The position of the Campbell Hill fault is shown in Burton and Harte (2013) approximately 900 ft to the east of that shown in Bennett and others (2006) and in figure 1.

Fractures are predominantly northeast striking and steeply inclined (more than 45 degrees [°]; Weston Solutions, Inc., 2012; Burton and Harte, 2013) and dip to the northwest. The placement of the area of interest on the western limb of the Massabesic anticlinorium is likely a major factor in the number of foliation features and fractures that dip to the northwest (Burton and Harte, 2013). Water-yielding fractures generally occur every 30 to 80 ft in a borehole; fractures, whether water bearing, low permeability, or unconnected, occur every 5 to 10 ft. Although some wells drilled 500 to 600 ft into bedrock have negligible water-bearing yield.

Lineaments are linear features seen on the surface of the Earth with remotely sensed imagery. Lineaments are potentially related to structural features and high-water yielding zones within the bedrock (Mabee and Hardcastle, 1997). Lineaments were mapped in the area coincident with the area-wide model by Clark and others (1997). For the nearby Pinardville quadrangle (fig. 2), Moore and others (2002) documented statistical differences in well yield (higher) with wells along lineaments that coincide with the regional fracture directions observed at outcrops (primarily trending north-northeast).

The hydraulic conductivity of the rock is inferred from residential water well yields reported from driller logs because of widespread availability of driller logs in the State. Moore and others (2002) documented relatively high well yields associated with tonalite rocks in the USGS 7.5-minute Pinardville quadrangle directly to the northeast of the model area. They also determined tonalite rocks statewide had relatively high yields. Tonalite rocks extend from the Pinardville to the Milford quadrangle (fig. 2). The model area comprises the northern one-half of the Milford quadrangle. Wells in the gneiss rocks tended to have low yields.

Well yields normalized by the open borehole length of the well are shown in figure 2. Normalizing yields per open length provide insight into the bulk, nondirectional, hydraulic conductivity of the rock. Although not shown, lineaments do not appear to play a factor in well yield and by association the hydraulic conductivity of the rock except for the lineaments mapped coincident with the Campbell Hill Fault. The residential water wells north of OU1 appear to have a higher bulk hydraulic conductivity than other wells. These wells are proximal to the Campbell Hill fault and the contact between upland tonalite and valley granite and gneiss.

Harte (2013) sampled several residential wells near OU1 for chemical and isotopic analysis (fig. 1). These residential well data were used to infer the degree of mixing in boreholes and the origin of groundwater recharging the wells. Residential water well water was determined to have characteristics of being derived from short and long flow path lengths (Harte, 2013). Long flow paths were hypothesized as potentially being able to transport PCE given the distances of residential water wells from contamination areas. However, several of the wells with water characterized as originating from long flow paths were adjacent to the hillside where it was more likely that flow paths came from the western hillsides rather than from PCE-contaminated areas to the south and east.

The surficial aquifer of the Milford quadrangle, which covers the entire model area, was mapped by Koteff (1970). The Milford-Souhegan River Valley is directly underlain by alluvium and glacially derived sand and gravel that in turn overlie thin, discontinuous glacial till mantling bedrock (Harte, 2010). The sand and gravel and till units are discontinuous, but one of the two is generally present because bedrock exposure is limited. The highly permeable sand and gravel ranges from 50 to 130 ft in thickness and constitutes most of the MSGD aquifer (Harte and Mack, 1992; Harte, 2010). The horizontal hydraulic conductivity of the most permeable units was simulated at 450 feet per day (ft/d) in the MSGD aquifer model by Harte and others (1999). North and south of the Milford-Souhegan River Valley, the landscape is dominated by abundant glacial deposits, including till and glacial outwash, with relatively sparse bedrock outcrops (Koteff, 1970). The till has a hydraulic conductivity typically less than 10 ft/d (Harte, 2010).

Remedial Systems

The unconsolidated sediments of the MSGD aquifer have been undergoing cleanup by a pump and treat method to help contain and capture the dissolved PCE plume that is associated with the Savage Municipal Water-Supply Well Superfund site. Pump and treat cleanup started in 1999 in OU1 and in 2004 in OU2 (Gradient Corp., 2011). In addition to the OU1 pump and treat system, a low permeability, slurry barrier wall encircles areas of potential PCE source areas in the form of dense nonaqueous phase liquids or sorbed contaminated sediments. The barrier wall was installed in late 1998 and sits atop the bedrock or is keyed several feet into the till. The location of the barrier wall, OU1 pumping and injection points where treated effluent is injected back into the aquifer, and selected OU2 pumping and injection wells are shown in figure 1.

Previous and Current Numerical Modeling Investigations

Several numerical models of flow and transport have been constructed for the unconsolidated sediments of the MSGD aquifer (table 1). However, previous models did not explicitly simulate the underlying crystalline-rock aquifer. Previous models constructed for the Milford-Souhegan River Valley included Harte and Mack (1992) and for the western one-half of the valley that focused on the Savage Municipal Water-Supply Well Superfund site previous models included Harte and others (1999). The areal extent of previous models and the area-wide model of flow are shown in figure 3. All models use some form of the USGS three-dimensional finite-difference groundwater-flow model (MODFLOW; U.S. Geological Survey, 2014).

The areal extent of models is informed by the modeling objectives of a study. For previous models, simulations focused on the MSGD aquifer and the PCE plume contained therein; the areal and vertical extents were limited to the unconsolidated sediments in the Milford Souhegan River Valley (fig. 3; table 1). With the discovery of PCE in deep fractures in the underlying crystalline-rock aquifer, the need to understand and quantify flow and transport in the crystalline-rock aquifer became crucial in the assessment of strategies for remediation. The area-wide model simulates flow in the crystalline-rock aquifer and covers a much larger area than previous models with the goal of improving the computation of groundwater flow from distal locations to the residential water wells and the area shown in figure 1. The model is also aligned with the primary fracture strike direction as mapped by Burton and Harte (2013).

Simulated Units and Hydraulic Properties

Simulated units include the glacial stratified-drift deposits (termed sand and gravel in this report), glacial till, and the underlying bedrock. The areal extent of the sand and gravel unit was restricted to land surface altitudes less than 330 ft to match general patterns observed by Koteff (1970) in the Milford-Souhegan River Valley. The areal extent of the till was primarily continuous. The areal extent of the bedrock was simulated as continuous throughout the active model area. Some sand and gravel units have been mapped by Toppin (1987) in adjacent river valleys above altitudes of 330 ft; however, these units are not simulated in the area-wide model because they are assumed to have negligible effect on flow in the area of interest (fig. 1).

Each hydrogeologic unit (sand and gravel, till, and bedrock) was assigned unique hydraulic properties. Simulated hydraulic properties included hydraulic conductivity (horizontal and vertical hydraulic conductivity), storage coefficient for transient simulations, and porosity for simulations of advective transport.

The hydrogeologic units (sand and gravel and till) initially were specified as horizontally homogeneous. However, during the calibration process and to better replicate previously published models of the MSGD, spatial heterogeneity was introduced into the sand and gravel unit. The horizontal hydraulic conductivity (k_h)of the sand and gravel unit typically was set at least one order of magnitude greater than that of the till, following regional and local hydraulic information on stratified drift deposits and tills (Toppin, 1987; Melvin and others, 1992; Harte, 2010). The k_h of the till was initially set at 3 ft/d, similar to k_h rates for the till in previously published MSGD aquifer models. The k_h of the bedrock was initially set at one order of magnitude less than the till (0.3 ft/d) and initially specified as horizontally homogeneous. Similar to the sand and gravel unit, spatial heterogeneity was incorporated into the bedrock during the calibration process for selected simulations. This spatial heterogeneity was introduced into model layers of the bedrock by incorporating k_h values calculated from hydraulic testing with straddle packers of bedrock fractures detected with borehole geophysical logging (Weston Solutions, Inc., 2012). The k_h of the bedrock was assigned to appropriate model layers based on vertical positioning. Between bedrock wells, the assigned k_h was interpolated by using a point inverse distance algorithm (Winston, 2009). Bedrock k_h also was adjusted based on rock type (Bennett and others, 2006) and the Campbell Hill Fault location (Burton and Harte, 2013). Bedrock k_h was not adjusted based on mapped lineaments (Clark and others, 1997) because of the lack of observable trends between lineament locations and well yield for the Milford area, excluding lineament features associated with the Campbell Hill Fault.

The vertical hydraulic conductivity (k_v) of the sand and gravel was set at 0.10 to 0.25 of the k_h, respectively, to account for layered units found in the MSGD aquifer (Harte, 2010). The k_v for the till was set the same as the k_h of the till. The k_v for the bedrock layers was initially set equal to k_h. Many of the water-bearing fractures are steeply inclined (more than 45° slope), which facilitates vertical flow (Weston Solutions, Inc. 2012). During calibration, the k_v of the bedrock was increased to a multiple of the k_h to better simulate inclined fractures.

Simulated Conditions

A long-term average annual preremediation (before the barrier wall and remedial system were installed in OU1) period (pre-1998) was simulated under steady-state conditions (long-term mean simulation). A quasiequilibrium simulated period, which represents the time before an aquifer test (before December 6, 2010) but postremediation (after the barrier wall and remedial system were installed in OU1) period, was simulated also under steady-state conditions (early December 2010 simulation). Finally, a short-term dynamic period representing an aquifer test that was completed during December 6 and 9, 2010, when pumping of well BR–6 (fig. 1) occurred was simulated under transient conditions (short-pump test simulation). Results from the December 2010 simulations were used as the initial conditions for the short-pump test simulations.

Discretization

The model boundaries and subsequently discretization is aligned with the principal strike and dip direction of the fractured rock. Model columns are parallel to strike and rows parallel to dip direction. The horizontal grid spacing in the model ranges from 110 to 220 ft in the row direction and 73 to 220 ft in the column direction (fig. 4). The smaller cell sizes cover the area of interest (fig. 1). The total number of cells per layer is 36,966 (202 rows by 183 columns), which totals 295,728 cells (8 layers) for the model.

The model was constructed with eight model layers (fig. 5). The sand and gravel unit is assigned to layer 1 at altitudes less than 330 ft. The thickness of the sand and gravel unit varies from discontinuous (zero) to more than 100 ft. The till unit is assigned to layers 1 and 2 at altitudes equal to and greater than 330 ft. The thickness of the till unit typically varies from 1 to 3 ft in the valley. In some river valleys, the till is actually absent; for those areas in the model the till is specified as a thin layer (less than 1 ft) to avoid creating a discontinuity in simulated layers. Specifying a thin till layer has little effect on the vertical flow between the sand and gravel unit and the bedrock because flow is primarily controlled by the hydraulic properties of the upper bedrock. No adjustment of hydraulic properties of the till was specified for this condition besides specifying a thin layer.

The underlying bedrock is assigned to the bottom six layers. The thickness of the uppermost bedrock layer (model layer 3) is 30 ft, which represents the maximum length of well casing penetrating the bedrock; typically, casing is 20 ft long. The well casing is impermeable, and no residential water-supply wells are simulated as open to layer 3. Bedrock layers 2 through 5 (model layers 4 through 7) are 100 ft thick. The lowermost bedrock layer (model layer 8) is 200 ft thick. The maximum simulated thickness is approximately 800 ft.

Hydrologic Boundaries

The area of the model is larger than the area of interest so that all lateral model boundaries should have a negligible effect on groundwater flow in the area of interest (fig. 1). The capture zones (zones contributing flow) of the residential water-supply wells to the north of OU1 likely extend to the west, into the uplands. The ability to delineate capture zones would be adversely affected if capture zones intersect inaccurately assigned lateral boundaries, such as a no-flow boundary where appreciable flow might occur. If affected by model lateral boundaries, capture zones might inaccurately extend farther to the south increasing the likelihood of capturing PCE-contaminated groundwater. Therefore, extending the model allows for an improved assessment of the potential capture of PCE-contaminated groundwater by the residential water-supply wells. Risk susceptibility based on potential capture of contaminated groundwater by residential water wells can further inform decision makers on the appropriateness of certain remedies in reducing risk.

The upper boundary of the model is a specified-flux condition to simulate areal recharge. The lowermost boundary is considered a no-flow boundary based on the hypothesis that open, hydraulically active fractures below this depth are less frequent than above this depth and that these deep fractures do not contribute appreciable flow to the simulated flow system. In addition, the model assumes that no appreciable deep underflow from outside the simulated area discharges or interacts with flow in the model area. The model edges are also no-flow boundaries except for the southeastern corner of the model, which is a head-dependent flux boundary for all model layers (fig. 4). Most model edges are aligned with local topographic highs where it is hypothesized that the underlying water table is high, forming a divide, and a subdued reflection of the topography because of the low bulk permeability of the bedrock (Harte and Winter, 1995). For the southeast corner edge boundary, specifying all layers as a head-dependent flux boundary promotes lateral groundwater flow and reduces upward discharge from the bedrock to the overlying unconsolidated sediments in that area. The actively simulated area along the southeastern part of the model boundary coincides with drainages that flow into the Milford-Souhegan River Valley and adjacent river valleys to the south, therefore some lateral flow is accounted for in the model.

Streams are a hydraulic boundary. Perennial streams were simulated with the River package of MODFLOW–2005 (Harbaugh, 2005). Digital stream coverages were obtained from NH Granit (2012) databases. Stream stage, streambed altitude (determined from water depth), and conductance (channel width, length, bed thickness, and hydraulic conductivity of bed material) were assigned based on previously collected data (Harte and others, 1997). In general, streams above an altitude of 330 ft were assigned narrow channels (5 ft across), shallow depth of water (less than 2 ft), and thin streambeds (1 ft thick). Conversely, streams below an altitude of 330 ft were assigned wide channels (as much as 60 ft across), greater depths of water (as much as 4 ft), and thick streambeds (3 ft thick). The hydraulic conductivity of the streambed was set at 3 ft/d in most areas based on work by Harte and others (1999).

Recharge and Groundwater Withdrawals

Areal recharge from precipitation and snowmelt was specified based on altitude. Where the land surface altitude was more than 330 ft, recharge was initially set at 0.5 foot per year (ft/yr). Where the land altitude was less than or equal to 330 ft, recharge was initially set at 2 ft/yr. Rates of recharge match published rates from Harte and Winter (1995) for upland areas and Harte and others (1999) for valleys. Recharge rates were allowed to vary during calibration for long-term mean annual and December 2010 simulations but were constrained in that the valley recharge was at least three times that of the uplands. For short-pump test simulations, recharge was applied during the first stress period (steady-state) but not applied during the subsequent transient stress periods to match the dry weather conditions immediately before and during the aquifer test.

Groundwater withdrawals varied based on the period of the simulation (table 2). For long-term mean annual simulation, groundwater withdrawals were simulated from residential water-supply wells drilled into the bedrock, withdrawals from wells at the New Hampshire Fish and Game Department’s Milford Fish Hatchery at Milford, N.H., and withdrawals from several other commercial wells. Rates of withdrawals were specified from New Hampshire water use reports (New Hampshire Department of Environmental Services, undated; Genevieve Al-Egaily, New Hampshire Department of Environmental Services, written commun., August 2, 2010), for commercial wells or estimated from typical per capita water use amounts for residential water-supply wells (Carolyn Hayek, U.S. Geological Survey, written commun., 2009). For the December 2010 simulations, additional groundwater withdrawals were specified to simulate the operation of remedial wells. For the short-pump test simulations, an additional groundwater withdrawal was specified for several days to simulate an aquifer test at well BR–6 inside the OU1 barrier. The short-pump test simulation uses the starting heads from the December 2010 simulation.

The average withdrawal for a residential water-supply well is low (0.2 gal/min; New Hampshire Department of Environmental Services, written commun., 2011). Each residential water well was simulated in all model simulations (table 2). Withdrawal rates represent an average of the cyclical nature of residential water use. There are 2,220 residential water wells that have a combined withdrawal of 444 gal/min; the effect of residential withdrawals is dispersed over a large area. Returned flows from septic systems were not included but considered within the margin of error in estimated withdrawals and indirectly simulated by the recharge parameter in layer 1 of the model. Groundwater is withdrawn from layers coincident with the open interval of the well as reported in well reports. Most residential water-supply wells are open to bedrock layers 2 and 3 (model layers 4 and 5).

Low-Permeability Barrier

A low-permeability barrier constructed of bentonite clay penetrates the full thickness of the sand and gravel unit and partially to fully penetrates the till in OU1 (Harte, 2006). Where the barrier fully penetrates the till, it sits atop the bedrock. In some places, the barrier depth exceeds 100 ft from land surface. The barrier encircles an area of 169,506 square feet (fig. 1).

The barrier was installed in 1998 and is specified in December 2010 and short-pump test simulations. The model simulates the barrier with the Horizontal Flow Barrier package of MODFLOW–2005 (Hsieh and Freckleton, 1993). The k_h of the barrier was set at comparable values (0.071 to 0.000046 ft/d) as that simulated in the MSGD aquifer models (Harte, 2006). The thickness of the barrier was set at 3 ft in the sand and gravel unit (model layer 1) and at 1 ft in the till unit (model layer 2) to approximate partial penetration of the barrier into the till at some locations. During calibration, model-computed head gradients inside the barrier were compared with observed head gradients to ensure reasonable rates of simulated flow across the barrier.

Solution of Model

During convergence, some cells that represent thin till layers (less than 2 ft thick) numerically dewatered and dried up prematurely, causing the cells to go inactive. Typically, this occurred on less than 100 cells or 0.03 percent of model cells. Where the cells numerically dewatered, heads and flows in adjacent cells were checked to ensure reasonable solutions and that the premature dewatering did not cause pervasive model error; none were observed. Additionally, numerical water budgets (volumetric inflows and outflows) were balanced to zero, indicating that a good overall numerical solution was achieved.

The Geometric Multigrid Solver package of MODFLOW–2005 (Wilson and Naff, 2004) was used to solve the finite-difference equations for the MODFLOW model. The Geometric Multigrid Solver package was determined to be a more robust solver for the automated sensitivity and parameter estimation routines invoked for this study than other solvers, such as the Preconditioned Conjugate Gradient solver (Hill, 1990).

Parameterization

Selected model input was grouped into parameters (table 3) to facilitate model calibration and analysis. Parameterization helps facilitate model calibration, sensitivity analysis, and parameter estimation (Poeter and others, 2005).

Sensitivity Analysis of Parameters

Understanding the sensitivity of model parameters as related to field observation is important to gage the relative effect of parameter adjustments to model results and model fit. UCODE_2005 (Poeter and others, 2005) allows for automated sensitivity evaluation. A particularly useful sensitivity analysis is the ranking of parameter importance associated with the observed dataset based on the composite scaled sensitivities (CSS). The CSS value reflects the total amount of information provided by the observations for the estimation of one parameter (Hill and Tiedeman, 2007). A high CSS value means that the information provided by the observations gives more confidence that the parameter can be estimated than a low CSS value.

The horizontal anisotropy of the bedrock is the best-informed parameter (most sensitive) as measured by the CSS in the ability of the model to replicate results (match simulated transient heads to observed heads) from the short-pump test simulations. A CSS graph (fig. 6) from a typical example of a short-pump test simulation shows that the HANI_PAR1 parameter (table 3) is ranked the highest in its ability to affect model results and replicate transient head data from the observed data associated with the short-pump test simulation. Observed drawdowns from the short-pump test simulation show a strong preferred northeast to southwest pattern that aligns with the predominant strike direction of the rock (Burton and Harte, 2013). Recharge is the most informed parameter for the calibration of model-computed heads to observed heads from residential water wells because heads from the residential water wells are spread over the model domain and recharge controls heads, particularly in upland areas. The k_h of the middle to deep layers of the bedrock is important for model fit to the short-pump test simulation because of the preponderance of water levels from the middle to deep bedrock layers. The storage coefficient of the bedrock is ranked the seventh most informed parameter overall, but in terms of model fit to the short-pump test simulation, storage coefficient is equally as important as the k_h of the middle to deep bedrock layers. Because the residential water levels and some of the other water-level data represent one value in time, those data have no effect on the storage coefficient parameter. The least informed parameters are the k_v of the upper units (sand and gravel and till units), which indicate parameter estimation of the k_v for this model cannot be solved.

The steady-state December 2010 simulation includes observations from river seepage (Harte and others, 1997). Valley recharge (parameter RECH_MULT) has the highest CSS indicating that valley recharge is the most informed parameter by the river observations used in this simulation. Figure 7 shows an example of typical CSS results from a steady-state December 2010 simulation. Upland recharge (parameter RECH_UP) and bedrock anisotropy (parameter HANI_Par1) are the most informed by residential heads (water levels).

Calibration Results

There are three simulation sets in this report (long-term mean annual, December 2010, and short-pump test). In discussing the simulation, different parameter variations are presented other than the calibrated version in order to highlight parameter sensitivity. The calibrated version per simulation set is always referred to as the calibrated simulation or base simulation for either long-term mean annual, December 2010, and short-pump test (Harte and Watt, 2021).

The long-term mean annual simulations were preliminary in nature and focused on bulk adjustments of hydraulic properties (k_h, k_v), recharge, and river properties (stage, bed, and hydraulic conductance). Layers were kept homogeneous, horizontally isotropic, and vertically anisotropic. Results from the long-term mean annual calibrated model indicate that model-computed heads are slightly higher than observed heads at low altitudes but slightly lower than observed heads at high altitudes (fig. 8). The automated calibration yielded a residual mean head difference (observed minus model-computed heads) of –6.3 ft, a median head difference of –1.4 ft, and an absolute mean head difference of 30.6 ft. Parameter estimation computed a k_h of 0.154 ft/d for the uppermost bedrock layer and 0.0038 ft/d for the lowermost bedrock layers. Confidence intervals for both values are large (infinity), suggesting high uncertainty in parameter estimation.

The high uncertainty indicates that additional observation data are needed to constrain parameter estimation and improve the calibration. Parameters from the calibrated model of the long-term mean annual simulations were initially used as input into models of the December 2010 and short-pump test simulations and additional head data were incorporated into the December 2010 simulation. Further, to help the calibration, various levels of heterogeneity and anisotropy were evaluated in a series of simulations to minimize the bias in model-computed heads at high altitudes, as shown by the residuals. The following discussion describes the changes to the short-pump test and December 2010 simulations to evaluate different parameter complexities and improve calibration.

Parameters from the calibrated model of the long-term mean annual simulations were initially used as input into models of the December 2010 and short-pump test simulations. The initial parameter estimates from the calibrated model for the long-term mean annual simulations were used in the short-pump test simulation (called “isotropic, layered homogeneous”; table 4). Heterogeneity was incorporated into the model by interpolating hydraulic conductivities, calculated from the packer hydraulic tests (Weston Solutions, Inc., 2010) between wells.

Head residuals (observed minus model-computed heads) for horizontal isotropic conditions, either layered homogeneous or heterogeneous, are large (table 4). When horizontal anisotropy in k_h was incorporated into the model at a ratio of 20:1 (short-pump test-anisotropic heterogeneous simulation) in the row (y) to column (x) direction to approximate preferential flow along the prevalent fracture strike, head residuals greatly decreased. Head residuals changed by relatively small amounts when incorporating heterogeneity into the model from the packer hydraulic test data. Results of the short-pump test simulations indicate that horizontal anisotropy is one of the most important parameters controlling flow in the bedrock.

Anisotropy ratios of 10:1 and 20:1 (y to x direction) produce an elongated potentiometric head contour that extends to the northeast as a result of withdrawals at well BR–6. The heads are shown in figure 9A for the contour at 240 ft altitude above msl for model layer 4 (bedrock layer 2) at the end of short-pump test simulation (simulations identified in table 4 as isotropic and homogeneous, and anisotropic and heterogeneous). The 240-ft contour is a representative contour that illustrates the variations in model-computed heads when compared with results from multiple simulations. Simulations presented highlight differences between isotropic and anisotropic conditions. Model-computed heads from other model layers (layers 5, 6, and 7) show greater drawdown than the other layers and model-computed heads show that simulated heads vary based on the k_h for those layers. The observed heads, where water levels were measured, are labeled next to the monitoring well (fig. 9A).

Iterative calibration between the short-pump test and December 2010 simulations indicated that, at an anisotropy of 20:1 ratio (y to x direction), head residuals were larger than using a 10:1 ratio for the December 2010 simulations alone. The final calibrated model for the December 2010 and short-pump test simulations, which uses an anisotropy of 10:1, shows that, at the end of the short-pump test simulation, a groundwater divide forms near bedrock well BR–1. At the divide, flow diverges to the west back to the withdrawal well (BR–6) and to the east following the valley-wide direction of flow (fig. 9B). The location of the divide is beneficial in that it contains contaminants in the bedrock underlying OU1 from lateral transport to the east. The distribution of observed heads also corroborates that a divide forms in this area.

The calibrated results with a fully anisotropic heterogeneous model for the December 2010 simulation show an improvement in model fit (fig. 10) and less altitudinal bias compared with the long-term mean annual simulation (fig. 8). Model-computed heads in the December 2010 simulation match observed heads better at high altitudes (more than 700 ft above msl) unlike the long-term mean annual simulation. The automated calibration yielded a mean head difference (between observed and model-computed heads) of 14.5 ft, a median head difference of 3.7 ft, and an absolute mean head difference of 34.1 ft.

Model-computed river seepage in the December 2010 simulation model matches well (relative percent difference of –2; table 5) with the observed river seepage for the second reach of the main stem (MS2) of the Souhegan River (app. 1). Model-computed river seepages for several tributaries match poorly (table 5). Less weight was applied in matching simulated to observed river seepages for the tributaries because of smaller flows.

Parameters were estimated using a multiplicative factor to the heterogeneity distribution derived from the packer hydraulic test (Weston Solutions, Inc., 2010) and adjusted for rock type and altitude. The k_h in areas underlying land surface altitudes more than 330 ft was set less than the k_h less than 330 ft. The range in k_h specified in the final models for December 2010 (and short-pump test) simulations is provided in table 6.

The uncertainty in parameter estimates are generally large (several orders of magnitude), indicating that model reliability could be improved with additional transient head and seepage data. In particular, river seepage data in upland areas would provide considerable leverage to constrain model parameterization.

Model-computed potentiometric heads in the bedrock for the calibrated December 2010 simulation model range from approximately less than 240 ft to more than 810 ft above msl (fig. 11). In the Milford-Souhegan River Valley, the potentiometric heads in the bedrock follow a similar pattern as do the heads in the MSGD aquifer as identified by Harte and others (1997). The head residuals are visible in figure 11 as a range of colors at wells with available water-level data; spatial bias is visible as a clustering of similar colors. Some clusters of similar colors occur in several areas, indicating head residuals differ by a similar amount; additional spatial parameterization of the bedrock k_h could be incorporated into the model to help improve calibration. An area with large head residual occurs at the southeastern corner of the model (fig. 11).

Model-computed heads in the valley for the calibrated December 2010 simulation model tend to be slightly higher than observed heads by 3 to 5 ft. Given that model-computed river seepage for the main-stem Souhegan River matches observed seepage, which is indicative that simulated rates of recharge in the valley are reasonable, it is possible that the hydraulic conductivity of the bedrock and unconsolidated sediments could be further refined in the valley to improve the calibration.

Alternative Conceptualization of the Crystalline-Rock Aquifer

Fracture length is an important control on flow in the bedrock aquifer and is an unknown variable. Fractures that truncate without intersecting an adjacent fracture would result in low permeability zones with little flow. To simulate finite fracture length, a simplified adjustment was made to the model where the anisotropy in the model was decreased by a factor of 100 for every 20th model row, or approximately 1,400 ft over the area-wide model. Simulated pathlines (not shown) moved to the northeast as in the base case but in a convex pattern, coinciding with the low anisotropy row where the fracture was hypothesized to terminate. At fracture terminal locations, pathlines moved upward, which would ultimately result in additional upward discharge to the sand and gravel and till units at fracture terminal locations and a corresponding decrease in groundwater flow in the bedrock. It is likely that many fractures terminate at cross-cutting connections with other fractures, which would cause sharp angular pathlines to form instead of upflow. Nevertheless, this simplified modification helps illustrate the importance of fracture lengths on flow.

Horizontal anisotropy, as shown in simulations of short-pump test, is a major control on flow in the crystalline-rock aquifer. Simulations with and without anisotropy show a distinctly different pattern of pathlines (figs. 12 and 13, at back of report). Without anisotropy, many pathlines discharge to tributaries to the Souhegan River farther east in the valley (fig. 13, at back of report) than in the base case (fig. 12, at back of report).

Changes in Hydrologic Conditions

Adjusting recharge rates by ±50 percent from the base case has a small effect on pathlines (not shown) that were forward tracked from PCE-contaminated areas underlying OU1. The general path line trajectory and patterns look similar to the base case, particularly for the high recharge model. However, low recharge causes some pathlines to move to the west by 100 ft, which puts the projection of the pathlines for the low recharge simulation closer to some residential water wells than pathlines for the base case. Recharge from the upland areas west of the residential water wells is likely an important control on limiting the northwestern transport of contaminated PCE from OU1.

Changes in Well Withdrawals

Several simulations were completed where well withdrawals were adjusted based on potential future scenarios. Scenarios included increased residential development and withdrawals at new residential water-supply wells, cessation of withdrawals at the Milford Fish Hatchery withdrawal wells, increase in withdrawals at OU2 (OU2–EW1 well), termination of remedial operation in OU1 and OU2, and initiation of a pump and treat remedy in the crystalline-rock aquifer underlying OU1.

The addition of new residential water-supply wells north and near the river affects advective transport of PCE contaminants. A scenario where four new hypothetical residential water wells were added to the base case was completed to represent the potential expansion of housing in that area (fig. 14, at back of report). The wells were simulated to depths of 400 ft bls and withdrawals were specified at rates similar to other residential water wells. For this scenario, some of the pathlines forward-tracked from contaminated areas discharge to the hypothetical residential water wells. The remaining pathlines either discharge to the Milford Fish Hatchery withdrawal wells or discharge to the remedial wells in OU1. Depending on the proportion of noncontaminated to contaminated PCE groundwater discharging to the hypothetical residential water wells, PCE concentrations in some wells could exceed maximum contaminant levels.

The hypothetical cessation of withdrawals at the Milford Fish Hatchery withdrawal wells (as of 2013 withdrawals were approximately 1,800 gal/min) causes some pathlines to discharge farther north toward residential water wells in the bedrock. Some pathlines discharge near tributaries to the north near high-yielding residential water wells in this area, whereas others discharge to the Souhegan River to the east (fig. 15, at back of report). A few pathlines discharge to remedial extraction wells in OU2. The overall pattern of pathlines is complex and is affected by multiple hydrologic sinks (withdrawal wells and discharge to rivers).

The hypothetical termination of remedial injections and extractions at OU1 and OU2 has a small effect on pathlines (not shown). For this scenario, pathlines show similar movement as in the base case excluding the pathlines that discharge to the Milford Fish Hatchery withdrawal wells, which formerly were captured by OU1 remedial extraction wells in the base case.

Increasing the rates of extraction at remedial wells in OU2 at EW–1 from 0 to 125 gal/min causes some pathlines to discharge to OU2 extraction wells (EW–1 and EW–2). Further, it has a ripple effect and causes a decrease in flow velocities by the Milford Fish Hatchery withdrawal wells that allows the partial capture of pathlines at one of the residential water wells (fig. 16, at back of report). However, it should be noted that the interception of one particle pathline to a single residential well is viewed as not strongly supportive of a major shift in flow but rather a minor shift in flow.

The addition of well BR–6 as a long-term (to steady-state conditions) withdrawal well captures almost all pathlines from the contaminated area (fig. 17, at back of report). Pathlines track back toward well BR–6 and no pathlines discharge to the north. Few pathlines discharge to the OU1 remedial wells inside the barrier, instead pathlines are captured by extraction at well BR–6. Results of this simulation indicate that continuous withdrawals at well BR–6 could be effective in preventing or reducing the transport of PCE contamination in the bedrock.