Home Introduction Hydrogeology Simulation of GW Discharge -Equations -SEAWAT -X-Sxn. Model -Regional- Scale Model -Model Limits Conclusions References Appendix 1 Appendix 2 Appendix 3 Appendix 4 Plates PDF Version

Simulation of Ground-Water Discharge to Biscayne Bay

Regional-Scale Ground-Water Flow Model

The regional-scale model simulates the flow of ground water from January 1989 to August 1998 for most of Miami-Dade County and parts of Broward and Monroe Counties. Results from preliminary model runs suggested that 5 to 10 years would be sufficient to eliminate the effects of initial conditions. Model results early in the simulation period should be interpreted with additional caution because they may contain residual errors from inadequate initial conditions. The hydrogeologic stresses included in the model are flows from lateral boundaries, areal recharge, evapotranspiration, ground-water flow to and from canals, recharge from surface water in the Everglades, ground-water withdrawals from municipal well fields, and ground-water discharge to Biscayne Bay (fig. 15).

Estimates of ground-water discharge to Biscayne Bay were obtained through the development and calibration of the regional-scale, variable-density, ground-water flow model. A major assumption with this approach for estimating ground-water discharge to Biscayne Bay is that realistic discharge quantities will be obtained by calibrating the numerical model to observed conditions. This means that the accuracy of the discharge estimates is dependent on calibration, which is inherently limited. A sensitivity analysis with the calibrated model was used to bracket a reasonable range of probable discharge estimates.

Spatial and Temporal Discretization

To simulate ground-water flow to Biscayne Bay, a regularly spaced, finite-difference model grid was constructed and rotated so that the y-axis would roughly parallel the coast (fig. 21). Each cell is 1,000 by 1,000 m square in the horizontal plane. The grid consists of 89 rows and 71 columns, and the rotation angle from true north is clockwise 14 degrees. The purpose for rotating the grid is to align model rows with the principal direction of ground-water flow, which primarily is toward Biscayne Bay. When ground-water flow is parallel to one of the primary model axes, problems with numerical dispersion that result from solving the transport equation are minimized. Another reason for rotating the model grid is that future modifications to the model may require a higher level of discretization at the coast. A rotated model grid will allow the resolution along a flow line to be increased by dividing columns near the coast. This modification would be relatively simple to implement but would increase computer run times.

Accurate simulation of variable-density flow systems requires a higher level of vertical discretization compared to that required for simulating constant-density flow systems. Accordingly, the model grid consists of 11 layers, more layers than would be required for a typical ground-water flow model that does not simulate variable-density flow. During the initial stages of calibration, the upper layer was used to approximate overland surface-water flow in Everglades National Park and Water Conservation Areas 3A and 3B. This upper layer caused numerical instabilities and increased simulation times. As a result, the entire upper layer (layer 1) was inactivated early in the development of the model. Rather than remove this upper layer, it remains in the model in case future simulations require the representation of overland flow.

Layers 2 through 11 represent the Biscayne aquifer. The top of layer 2 is spatially variable and corresponds with land-surface elevation (fig. 22), based on a topographic contour map provided by Everglades National Park. For model cells that lie within Biscayne Bay, the top of layer 2 has a value of 0.0 m. The bottom of layer 2 is set at an elevation of 5.0 m below sea level. To minimize the effects of numerical problems that can occur in solute-transport models, the grid was designed so that most of the cells would have a uniform volume (1,000 x 1,000 x 5 m). Although the volume for each model cell in layer 2 may be slightly variable because of the variation in land-surface elevation, model cells in layers 3 to 11 have a uniform thickness of 5 m and thus a uniform cell volume. The bottom elevation for layer 11 is 50 m below sea level. The base of the Biscayne aquifer, as defined by Fish and Stewart (1991), generally slopes from west to east. This variability in aquifer thickness was included in the model by inactivating model cells with cell centers below the base of the Biscayne aquifer.

The simulation is divided into 116 monthly stress periods. For each stress period, the average hydrologic conditions for that month are assumed to remain constant. This means that the model does not simulate hourly or daily hydrologic variations, but rather seasonal and yearly variations. Further temporal discretization is introduced in the form of transport steps. Each stress period is divided into one or more transport steps, the lengths of which are determined by SEAWAT to meet specified criteria associated with solving the solute-transport equation. For the regional-scale model, about 20 transport steps were required for each stress period.

Assignment of Aquifer Parameters

The approach for assigning aquifer parameters that pertain to ground-water flow and solute transport was to use the simplest distribution that would result in adequate representation of the flow system. For the model parameter that has the largest effect on ground-water flow, horizontal hydraulic conductivity, most of the model was assigned a single value of 9,000 m/d, though two other zones were assigned lower values (pl. 3). The value of 9,000 m/d is similar to the value of hydraulic conductivity used by Merritt (1996a) for most of the model domain in southern Dade County to describe the Miami Limestone, Fort Thompson Formation, and permeable zones of the Tamiami Formation. For the northeastern and southern parts of the model, a relatively low value for horizontal hydraulic conductivity (3 m/d) was assigned to the upper layer to represent the peat and marl units that comprise the upper part of the Biscayne aquifer. Another zone, located near the central part of the model along the coast, was assigned a value of 1,500 m/d. Merritt (1996a) includes this zone to match the flow patterns and head values observed in this area. A constant anisotropy ratio of 100:1 was arbitrarily selected to represent the ratio of horizontal hydraulic conductivity to vertical hydraulic conductivity. Results from sensitivity analyses indicate that the anisotropy ratio does not significantly affect model results.

SEAWAT requires the input of primary (SF1) and secondary (SF2) storage values. The primary value is used in the flow equation when the head in a cell is above the top elevation of that cell. The secondary storage value is used for the water-table case when the water table lies within the model cell. For the uppermost active layer (layer 2), SF2 is equal to 0.2, the value of specific yield estimated by Merritt (1996a). SF1 is set to 1.0 for layer 2, a value that is appropriate when the water table rises above land surface. For layers 3 through 11, SF1 is 5.9 x 10^-5 and SF2 is 0.2. The storage value of 5.9 x 10^-5 was calculated using an average thickness of 30 m for the Biscayne aquifer and a value of 2 x 10^-10 m² (square meters) per newton for aquifer compressibility of fractured rock (Domenico and Schwartz, 1990, p. 111).

Aquifer parameters required for solute transport include porosity and dispersivity. A constant porosity value of 0.2 was assigned to each layer. This value was selected because it is thought that the specific yield, which was determined by Merritt (1996a), is indicative of the effective porosity. The longitudinal and transverse dispersivity values used for the model are 5 and 0.5 m, respectively. These dispersivity values were determined from the calibration of the cross-sectional models.

Boundary Conditions

For most simulations of ground-water flow, boundaries of the model are extended to locations in the aquifer where hydrogeologic boundaries reside. Ideally, these hydrogeologic boundaries are persistent flow lines, impermeable hydraulic barriers, or areas that can be represented with a constant flux or head. For this study, adequate hydrogeologic boundaries do not exist for the inland portion of the model domain. For this reason, the northern, western, and southern boundaries were extended outside of the area of interest and represented with head-dependent boundaries (pl. 3). The eastern boundary, numerically represented as a constant-head and constant-concentration boundary, corresponds to Biscayne Bay.

Biscayne Bay

The model cells in layer 2 representing Biscayne Bay (pl. 3) were simulated with the time-varying Constant-Head (CHD) package in SEAWAT and a constant salt concentration of 35 kg/m³. By specifying a constant freshwater head boundary in layer 2, the elevation for the bottom of Biscayne Bay corresponds to the center elevation of layer 2; thus, the simulated bay bottom is flat with an elevation of -2.5 m. Average monthly values for downstream stage at structure S-123 (fig. 14) were used to set the water level in the constant-head cells. The downstream stage measurement at structure S-123 was used because this structure is located at the coast (fig. 4) and lies near the center of the model. As previously shown, the downstream stage value at structure S-123 is a reasonable approximation for sea-level elevation (fig. 14). From September 1992 to January 1993, the downstream stage at structure S-27 (fig. 4) was used to replace missing data at structure S-123. The average value for monthly stage was converted to freshwater head using the specified salt concentration of 35 kg/m³ and an elevation of 2.5 m below sea level, which is the center elevation of the model cell. Monthly averages for the downstream stage at structure S-123 ranged from less than 0.05 m to greater than 0.45 m from 1989 to 1998. These relatively large fluctuations in average monthly sea level are included in the representation of Biscayne Bay to better approximate ground-water flow near the coast.

Inland Model Domain Boundaries

To represent the lateral flow of ground water to or from the inland perimeter of the model, head-dependent boundaries were assigned to each active cell (west of the Biscayne Bay shoreline) along the northern, western, and southern perimeters of the model (pl. 3). The General-Head Boundary (GHB) package in SEAWAT is used to represent head-dependent boundaries. The reservoir heads for the boundary cells were interpolated from a triangular irregular network (TIN) data model used to represent the water-table surface. A water-table TIN was created for each month of the simulation period (1989-98) using average monthly water levels from available surface-water and ground-water monitoring stations (app. III). Boundaries on the upper model surface in wetland regions will be discussed later.

Salt concentrations for boundary reservoirs were set equal to the initial concentration of the model cell (discussed in a later section). The concentration of the boundary reservoir remains constant throughout the simulation period. The concentration is zero for most GHB reservoirs, but for boundary cells near the coast, concentrations may be greater than zero. For boundary cells with concentrations greater than zero, reservoir heads were recalculated as equivalent freshwater heads using equation 8. The calculation and assignment of equivalent freshwater heads (according to eq. 8) result in hydrostatic conditions for those boundaries. The assignment of hydrostatic conditions to the northern and southern boundaries is probably inappropriate near the coast because coastal ground-water flow patterns suggest that there can be a large vertical flow component. However, there is not sufficient quantitative data to set these coastal boundaries with a more appropriate condition. Simulated ground-water flow patterns near these boundaries, therefore, should be interpreted with caution.

Lower Model Boundary

The lower model boundary represents the base of the Biscayne aquifer. This boundary is represented in the model as an impermeable barrier, an approach typically used in Miami-Dade County (Merritt, 1996a; Swain and others, 1996). Although there may be upward vertical ground-water flow from deeper aquifers into the Biscayne aquifer, there is no evidence to suggest that this is necessary to include in the model.

Internal Hydrologic Stresses

Hydrologic stresses that are internal to the model domain are represented with internal boundary conditions. Internal hydrologic stresses include: canals, surface water in the Everglades, recharge, evapotranspiration, and ground-water withdrawals from municipal well fields.

Canals

The dynamic exchange of water between the Biscayne aquifer and major canals is simulated with the River (RIV) package in SEAWAT. The RIV package assigns a head-dependent flux to each model cell that is intersected by a major canal. The flux value is a function of the difference between canal stage and aquifer head and the level to which the canal is hydraulically connected to the aquifer. In the model, canal stages vary spatially and temporally, depending on the stage values that were observed in the field. Simulated flows to and from the canal cells are compared with the measured or computed flows in the canals to ensure that the model is accurately representing the system. This process is described later. Although surface-water flow in the canals is not explicitly simulated, it is thought that this approach is reasonable for representing the interaction between canals and the Biscayne aquifer.

The canals included in the model primarily consist of major water-management canals (pl. 3). Minor canals, sometimes referred to as secondary or tertiary canals, generally were not included in the model. To facilitate model development, the major canal network was divided into individual canal segments based on hydrologic characteristics, network geometry, structure locations, and the hydrologic basin delineations by Cooper and Lane (1987). Each canal segment was assigned an identifying number, a canal bottom elevation, a canal width, a representative value for hydraulic conductivity between the aquifer and the canal, a representative value for flow length, and the upstream and downstream monitoring stations used to assign canal stages. Values for hydraulic conductivity and flow length were formulated using the conceptual model described by Swain and others (1996, p. 66). These values were modified during calibration so that simulated values of canal baseflow and aquifer heads would match those observed in the field.

Of the stress and boundary packages used in the regional-scale model, the RIV package was the most complex to develop. The RIV input file read by SEAWAT was created with the RIVGRID program (Leake and Claar, 1999). RIVGRID geometrically intersects a river (or canal) network with a model grid and calculates the hydraulic conductance for those model cells that contain a river. RIVGRID also calculates a stage for each RIV cell by linearly interpolating a value along the river segment.

For tidal canal segments to the east of coastal control structures, the stages calculated by RIVGRID are corrected to equivalent freshwater heads. In the freshwater head calculation, it is assumed that the surface water in the tidal canal is pure seawater with a salt concentration of 35 kg/m³. As required by SEAWAT, the elevation of the river bottom is used in equation 8 as the elevation term for the freshwater head calculation. A concentration of 35 kg/m³ also is used to represent the salt concentration of the water that enters the aquifer from the tidal canals.

Ponded Surface Water

To represent the hydrologic effects of surface water in the Everglades and coastal lowlands, a GHB was included in model cells of layer 2 that had standing water above land surface. To determine whether or not there was standing water in a model cell, a TIN of the water table was created for each month using average monthly stage data from canals and average water-table elevations from monitoring wells. A water-level elevation was then interpolated from the TIN for each model cell. A GHB that represents surface water is active in the model when the water level interpolated from the TIN is higher than land surface and a RIV boundary does not exist in the same model cell. Plate 3 shows the model cells in layer 2 that have an active GHB for at least one stress period of the simulation. The conductance value for each GHB, determined through calibration, is 1 x 10⁵ m²/d. The salt concentration of the GHB fluid is specified as 0 kg/m³.

Recharge and Runoff

Figure 23. Land use of southern Florida for 1995. [larger image]

The Recharge (RCH) package in SEAWAT is used to apply areal recharge to the model, with a specified salt concentration of 0 kg/m³. The general procedure for estimating recharge values was to subtract from rainfall a quantity thought to be representative of runoff. The remaining quantity is applied to the model as recharge.

For each month of the simulation period, Thiessen polygons were constructed using the rainfall monitoring stations that had continuous data for that month. Because few rainfall stations had continuous data for the entire simulation period, the geometry and distribution of Thiessen polygons varies by month. A monthly rainfall total was then estimated for each model cell by intersecting the model grid with the Thiessen polygons and weighting the rainfall totals by the areas of the Thiessen polygons within the model cell. To calculate a recharge value for each cell and for each month, the rainfall totals were multiplied by runoff coefficients. The calculation of runoff coefficients was based on the assumption that runoff quantities are dependent on land use, an approach used by Restrepo and others (1992) to construct a ground-water flow model for Broward County. For the present study, runoff coefficients were estimated for each land-use category through model calibration and results provided by Restrepo and others (1992). Table 5 contains the calibrated runoff coefficients that were assigned to each land-use category. These runoff coefficients may include the effects of other processes, such as pumpage from the Biscayne aquifer for agricultural irrigation, and may not be appropriate for other areas or spatial scales. An average runoff coefficient was calculated for each model cell by intersecting the model grid with the 1995 land-use coverage (fig. 23) and calculating the runoff coefficient based on a weighted-area average.

Evapotranspiration

The Evapotranspiration (EVT) package in SEAWAT, which withdraws ground water as a function of depth to the water table, is used to simulate evapotranspiration from the water table. The EVT package works by "pumping" ground water from the uppermost active model cell based on three, user-specified parameters (described in McDonald and Harbaugh, 1988): maximum evapotranspiration rate (QMAX), evapotranspiration surface (SURF), and extinction depth (EXTD). When the elevation of the water table is above the SURF elevation, the evapotranspiration rate is equal to QMAX. When the depth to the water table is greater than EXTD, the evapotranspiration rate is zero. When the water table is between SURF and EXTD, the evapotranspiration rate is linearly calculated between QMAX and zero. The SURF value assigned to each model cell is equal to the land-surface value. It is assumed that the extinction depth is related to land use (table 5; Restrepo and others 1992). The maximum evapotranspiration rates are assigned by month based on a study by Merritt (1996a). The monthly values for QMAX are listed in the table below. Evapotranspiration is not active for a model cell when a GHB (representing surface water) is active in that same cell. In that case, the effects of evapotranspiration are indirectly incorporated in the model through the specified stage of the GHB.

	January	February	March	April	May	June - October	November	December
Maximum evapotranspiration rate (centimeters per day)	0.20	0.28	0.36	0.43	0.46	0.53	0.30	0.28

Municipal Well Fields

The effects of ground-water withdrawals from municipal well fields are included in the numerical model using the Well (WEL) package in SEAWAT. Pumpage was assigned to the model by well field rather than by pumping well. Each well field was assigned to a single cell in the model based on well-field location and the open-hole interval of the pumping wells (pl. 3). In general, this simplification is reasonable because most of the pumping wells for a well field are contained within one model cell. Additionally, the model is not designed to characterize flow patterns at the well-field scale. Average values of monthly pumping rates for each well field were obtained from South Florida Water Management District (Water Use Division, written commun., 1999) and from Bertha Goldberg (Miami-Dade Water and Sewer Department, written commun., 1999). These pumping rates were assigned as monthly averages to the appropriate cells in the model.

Initial Conditions

Initial water levels were specified for the model by creating a TIN of the average water-table surface for January 1989. The TIN was generated using stage data from canals and water-level elevations from monitoring wells. Initial heads were assigned to each cell and layer in the model by interpolating values to the cell centers from the TIN of the water table.

Early attempts to specify initial concentrations consisted of specifying seawater concentrations (35 kg/m³) east of the 1995 saltwater intrusion line (Sonenshein, 1997) and freshwater concentrations west of the saltwater intrusion line. This procedure was used for each layer and results in a vertical wall of seawater at the 1995 saltwater intrusion line. Results from early simulations suggested that a better estimate for initial salinity concentrations would reduce the length of simulation time required to achieve a stable position for the saltwater interface, thus increasing the length of time model results would be unaffected by the initial salinity distribution. The initial salinity distribution was improved by linearly interpolating concentrations between the 1995 saltwater intrusion line and Biscayne Bay. The initial salinity distribution was further improved by incorporating the results from the airborne geophysical survey (fig. 11; Fitterman and Deszcz-Pan, 1998) for the southern part of the model area.

Calibration Procedure and Model Results

Calibration is a procedure of adjusting input parameters until the model is a reasonable representation of the physical system. For this study, calibration was achieved by adjusting the input parameters for the model within a reasonable range until simulated values of head, canal baseflow, net recharge, and position of the saltwater interface approximated values observed in the field. In addition to these quantitative comparisons, the general flow patterns simulated by the model were qualitatively inspected to ensure that the model reasonably represents ground-water discharge to Biscayne Bay. To ensure that the variable-density component of the model was working properly, special attention was focused on the cyclic flow patterns that result within the saltwater interface.

Comparison of Simulated and Observed Heads

Ground-water models traditionally are calibrated by matching simulated heads with water levels recorded in ground-water monitoring wells. Although this approach was used in this study, caution is exercised in evaluating the comparison because simulated and observed heads tend to match even when inappropriate aquifer parameters are used in the simulation.

Figure 24. Location of ground-water monitoring wells in southern Florida used to calibrate the regional-scale ground-water flow model and average differences between observed and simulated values of head. [larger image]

This potential problem results from two related conditions unique to the study area: (1) most of the model is overlain by hydrologic boundaries (canals and standing water in the Everglades), and (2) the hydraulic conductivity of the Biscayne aquifer is so high that the simulated heads are primarily controlled by boundary heads. Nevertheless, a reasonable match between observed and simulated heads is required for calibration, though a reasonable match does not ensure that the model is adequately calibrated.

Three different error statistics are commonly used to evaluate the comparison between observed and simulated values of head or other model variable (Anderson and Woessner, 1992). The mean error (ME) is simply the average of the differences between observed and simulated heads. Although the ME is a useful statistic, an ME value of zero does not indicate a perfect match between observed and simulated because positive and negative errors can average to zero. The mean absolute error (MAE) is calculated by taking the average of the absolute values of the differences between observed and simulated heads. The MAE may be the most useful statistic because it provides a true measure of the average difference between observed and simulated heads. The root-mean square error (RMSE) is calculated by taking the square root of the average of the squared differences. The RMSE is useful because it evaluates the spread of the errors by approximating the standard deviation of the errors.

The error statistics can be calculated in a variety of ways. For example, error statistics can be calculated for the differences by well and stress period, the differences for each well, or the differences for each stress period. The ME in head for all stress periods and wells (6,525 values in total) is 0.08 m. The positive value for ME indicates that simulated values of head generally are higher than observed values of head. The MAE is 0.15 m, which is an acceptable value considering that observed heads range from -2.23 to 2.51 m. This means that the MAE relative to the range in observed heads is about 3 percent. The RMSE is 0.27 m, indicating that about 68 percent of the head differences are within plus or minus 0.27 m of the ME. A histogram constructed from the differences between observed and simulated heads suggests that the differences are normally distributed with a mode of about 0.05 m.

Figure 26. Calibration statistics for errors between observed heads and heads simulated by the regional-scale ground-water flow model. [larger image]

The locations of monitoring wells used in the comparison between observed and simulated heads are shown in figure 24. Before head differences were calculated, the simulated heads were converted from freshwater head. In general, the ME is randomly distributed throughout the model domain (fig. 24). Wells near the coast tend to have an ME value that is relatively small, suggesting that simulated flow patterns near the coast are probably reasonable. To illustrate the level of calibration between observed and simulated values of head, the comparison is also presented as time-series plots for selected monitoring wells along the coast (fig. 25). These time-series plots, organized from north to south, illustrate the capability of the model to reasonably simulate head and fluctuations in head. A time-series plot of the error statistics also provides insight into the temporal distribution of error between observed and simulated heads (fig. 26). Errors tend to be large during the wet season and wet years (fig. 26).

Comparison of Simulated and Observed Canal Baseflow

Calibration of the ground-water flow model requires a reasonable match between observed and simulated values of canal baseflow, which is defined as the discharge of ground water to a canal. Canal baseflow can be positive or negative depending on location and the time of year. Matching canal baseflow is a required part of the calibration process because simulated and observed volumetric flow rates, rather than values of head, are compared. In the Biscayne aquifer where values of hydraulic conductivity are high, heads are relatively easy to match with a ground-water flow model. Matching baseflows, however, is much more difficult to achieve, but doing so results in a more reliable model.

Hydrologists divide study areas into smaller areas called surface-water basins. Runoff within a surface-water basin typically flows into a major river or canal and then exits the basin at a discrete outflow point.

Figure 27. Location of surface-water basins in southern Florida and the mean absolute error between observed and simulated canal baseflow. [larger image]

The surface-water basins shown in figure 27 were delineated for southeastern Florida by Cooper and Lane (1987).

Calibration of the regional-scale model to canal baseflow was performed for selected surface-water basins. For each of these basins, simulated baseflow was combined with the total runoff value estimated from the land-use method described earlier. This combined value was then compared with the net canal loss or gain measured for the basin. Net losses and gains were calculated for the selected surface-water basins by subtracting canal inflows from canal outflows. It is assumed that the difference between inflow and outflow is equal to baseflow and runoff. The mathematical expression is:

Baseflow + Runoff = Canal outflow - Canal inflow

where the left-hand-side is based on the simulation and the right-hand-side is based on measured quantities. Terms on the right side of the equation were obtained from the database of the South Florida Water Management District. Canal outflow is the sum of all surface-water discharge that flows out of the surface-water basin. Canal inflow is the sum of all surface-water discharge that flows into the surface-water basin. This approach assumes that direct rainfall to a canal, direct evaporation from a canal, and storage within a canal are negligible.

The calibration to canal baseflow was achieved by adjusting conductance values for the RIV cells. Plots showing a comparison between observed and simulated values of canal baseflow (fig. 28) suggest the model reasonably simulates the dynamic exchange of water between canals and the Biscayne aquifer. The plots in figure 28, organized from north to south, suggest that the model is reasonably capable of simulating the complex transfer of water between canals and the aquifer. Error statistics for the comparisons are included in figure 27.

Net Recharge

Net recharge is generally described as the percentage of rainfall that reaches the water table after runoff and evapotranspiration are withdrawn. Since evapotranspiration is a function of the water-table elevation in the model, the rate of net recharge actually is a result of the simulation. The simulated rates of recharge, therefore, must be inspected to ensure that values are within a reasonable range. This is particularly true for the regional-scale model because net recharge is a large component of the water budget.

A subroutine was added to the SEAWAT code to save net recharge rates by transport step. A post-processing code was then used to calculate the average value of net recharge for each cell based on the entire simulation period. Net recharge was calculated by combining rainfall, evapotranspiration, and discharge to or from a GHB representing surface water in the Everglades. The spatial distribution of average annual net recharge for the simulation period is illustrated in figure 29. Net recharge rates are expressed in centimeters per year to facilitate a direct comparison with the average annual rainfall rate of 141 cm/yr. Results suggest that the regional-scale model reasonably simulates the spatial distribution of net recharge. High recharge values are located along the Atlantic Coastal Ridge (fig. 2) where the water table is the deepest. In the Everglades and water-conservation areas, the model simulates moderate recharge values (fig. 29). Negative values of net recharge indicative of high evapotranspiration rates and ground-water discharge areas are simulated in the low-lying areas along the southeastern part of the model domain.

Saltwater Interface

For the regional-scale model to accurately simulate the discharge of ground water to Biscayne Bay, the model must accurately simulate the three-dimensional distribution of ground-water salinity. Unfortunately, data are lacking to adequately characterize the three-dimensional distribution of ground-water salinity because most monitoring wells are installed near the toe of the saltwater interface rather than within the transition zone. The observed data provide a two-dimensional description of the position of the saltwater interface. Rather than compare the few point measurements of ground-water salinity with simulated values of salinity, model results are reduced to two dimensions to facilitate a comparison with the 1995 position of the saltwater intrusion line (Sonenshein, 1997). In figure 30, the simulated salt concentrations in the lowermost active model cell are shaded according to their corresponding salinity values. Salinities in the lowermost active model cells are the most appropriate model results to use in the comparison with the observed 1995 saltwater intrusion line because the line was drawn for the base of the Biscayne aquifer.

In general, the model reasonably simulates the position of the saltwater interface, particularly south of the Miami Canal. North of the Miami Canal, the simulated position of the saltwater interface extends too far inland by 5 to 10 km. Three possible explanations for the discrepancy between observed and simulated positions of the saltwater interface for the area north of Miami Canal include: (1) the aquifer parameters and hydrologic stresses do not accurately represent the physical system, (2) numerical dispersion has caused the transition zone to spread too far inland, and (3) the observed position of the 1996 saltwater intrusion line is in disequilibrium with the current hydrologic stresses, but with time will migrate inland to the position simulated by the model.

Figure 31. Total salt mass in the Biscayne aquifer as simulated by the regional-scale ground-water flow model. [larger image]

Aquifer parameters and hydrologic stresses that did not accurately represent the physical system seem to be the most likely explanation for the discrepancy. Simulated heads at monitoring wells G-852 and F-45 generally are too low (fig. 25), suggesting that the hydrologic stresses or aquifer parameters used in the model may be in error. By improving the model in this area, the simulated position of the saltwater interface probably would more closely approximate the observed position of the saltwater interface.

To test the second possible explanation for the discrepancy (numerical dispersion), hypothetical, cross-sectional models with different levels of discretization were run with SEAWAT. The cross-sectional model with resolution similar to that of the regional-scale model (1,000 m horizontal by 5 m vertical) gave results identical to those from a much more detailed model (200 m horizontal by 1 m vertical), suggesting that numerical dispersion probably does not significantly affect the regional-scale model.

The third possible explanation for the discrepancy (disequilibrium) also is unlikely. Chloride concentrations in monitoring wells north of the Miami Canal do not show a long-term increasing trend, which suggests that the current position of the saltwater interface is in equilibrium with the current hydrologic stresses.

A plot of total salt mass in the aquifer relative to time suggests that the initial salt concentrations in the model are too low (fig. 31). However, as the simulation progresses through time, the total salt mass in the aquifer generally increases, with slight decreases during the wet seasons. Toward the end of the simulation, it appears that the model approaches steady state with respect to salt mass in the model domain. The plot of total salt mass suggests that the ground-water salinities simulated by the model are approaching equilibrium with the hydrologic stresses and aquifer parameters used in the model. Results from early in the simulation period should be evaluated with caution because initial ground-water salinities probably are too low.

Ground-Water Flow to Biscayne Bay

The simulated ground-water flow from the active model cells into the coastal constant-head cells is assumed to represent the discharge of ground water to Biscayne Bay. In addition, to simulating the volumetric flow rate, the model also simulates the salt concentration of the ground-water flow into the constant-head cells. This salt concentration is assumed to represent the salinity of the ground water that discharges to Biscayne Bay. The simulated salinity of ground-water discharge to Biscayne Bay ranges from nearly fresh at the shoreline to nearly seawater some distance offshore. To simplify the results, simulated discharge estimates are presented as the freshwater portion of the total discharge. The freshwater portion of the ground-water discharge is calculated from the total discharge using the following equation

,

(13)     (D)

where:

Qf	is fresh ground-water discharge, in cubic meters per day,
QT	is total ground-water discharge, in cubic meters per day,
S	is fluid density of pure seawater, in kilograms per cubic meter,
	is simulated fluid density of ground water discharging to Biscayne Bay, in kilograms per cubic meter, and
f	is fluid density of pure freshwater, in kilograms per cubic meter.

When the fluid density of the ground-water discharge is equal to 1,000 kg/m³, the fresh ground-water discharge is equal to the total ground-water discharge. When the fluid density of the ground-water discharge is equal to the fluid density of seawater (1,025 kg/m³), the fresh ground-water discharge is equal to zero. During certain times, simulated flow is from the bay into the aquifer. When this condition occurs, Q_T is negative, but Q_f is zero. A subroutine was added to SEAWAT to calculate and sum the Q_T and Q_f terms between each constant-head cell and the adjacent active cells. These flow terms are written to a file following each timestep. At the end of a model run, a post-processing routine calculates the average flows for each constant-head cell by stress period.

Figure 32. (above) Simulated fresh ground-water discharge to Biscayne Bay. [larger image]

Figure 33. (above) Simulated fresh ground-water discharge compared with measured surface-water discharge to Biscayne Bay. [larger image]

One potential problem with calculating fresh ground-water discharge estimates from simulated density is that the resulting freshwater discharge quantities are directly dependent on salt concentrations. (Density is a linear function of salt concentration.) As previously mentioned, ground-water salinities are considered calibrated when the simulated position of the saltwater interface matches the observed position of the saltwater interface. This does not ensure that the simulated salt concentrations are calibrated. The average salt concentration of simulated discharge to Biscayne Bay is about half that of seawater, and thus half of the total discharge is freshwater. This suggests that if salt concentrations are in error by 100 percent (17.5 kg/m³), estimates of fresh ground-water discharge would be in error by about a factor of two.

Results from the regional-scale model indicate that fresh ground-water discharge to Biscayne Bay occurs along the coastline and into the tidal portions of the Miami, Coral Gables, and Snapper Creek Canals (fig. 32). The model suggests that fresh ground-water discharge at the coastline of Biscayne Bay is similar in magnitude to the fresh ground-water discharge to the tidal canals. The average rate of fresh ground-water discharge is about 2.2 x 10⁵ m³/d for the coastline of Biscayne Bay, about 1.6 x 10⁵ m³/d for the tidal portion of the Miami Canal, about 1.4 x 10⁵ m³/d for the tidal portion of the Coral Gables Canal, and about 3.2 x 10⁴ m³/d for the tidal portion of the Snapper Creek Canal. The annual fluctuation in fresh ground-water discharge is about 1 x 10⁵ m³/d for the coastline of Biscayne Bay and the tidal portions of the Miami and Coral Gables Canals and is about 3 x 10⁴ m³/d for the tidal portion of the Snapper Creek Canal. Fluctuations in ground-water discharge appear to be damped because sea level was highest during the wet season when the water table was relatively high and lowest during the dry season when the water table was relatively low.

A comparison was made between simulated fresh ground-water discharge and measured surface-water discharge from the coastal control structures (figure 4 and figure 5, from S-29 to S-197). Based on the results for the simulation period (1989-98), fresh ground-water discharge seems to be an important mechanism of freshwater delivery to Biscayne Bay during the dry season (fig. 33). For the dry seasons of 1989, 1990, and 1991, model results suggest that fresh ground-water discharge exceeded the surface-water discharge to Biscayne Bay. During the wet season, however, fresh ground-water discharge is about one order of magnitude less than the surface-water discharge. For the total simulation period, ground-water discharge is about 6 percent of surface-water discharge.

Nearly 100 percent of the fresh ground-water discharge to Biscayne Bay is to the northern half of Biscayne Bay, north of structure S-123 (figure 4 and figure 5). South of structure S-123, where land-surface elevations are less than 1 m above sea level, the water table was unable to develop enough head to drive ground-water discharge into the bay.

Sensitivity Analysis

A sensitivity analysis was performed with the regional-scale model to determine the effects of various boundary conditions and aquifer parameters on the simulated estimates of fresh ground-water discharge to Biscayne Bay. Boundary conditions that were evaluated in the sensitivity analysis include stage in Biscayne Bay and the type of solute-transport boundary used to represent Biscayne Bay. Aquifer parameters that were evaluated with the sensitivity analysis include horizontal hydraulic conductivity, vertical conductance, dispersivity, and river conductance. Simulated values of freshwater discharge to Biscayne Bay do not appear to be sensitive to vertical conductance.

For the first sensitivity run, a different type of concentration boundary was used to represent Biscayne Bay. In addition to the constant-concentration boundary type, MT3D and SEAWAT have another solute-transport boundary type that may be appropriate for representing concentrations in Biscayne Bay. This other boundary type is referred to as specified concentration for inflow, which means that any fluid entering the active model domain has a specified concentration. With the constant-concentration boundary, salt mass is transferred from the boundary into the active model domain by advection and dispersion. With the specified concentration for inflow, however, salt mass is transferred only by advection. By specifying the concentration for inflow, the simulated concentration of the constant-head cell may be less than the actual concentration specified for inflow. Compared with the calibrated, regional-scale model (referred to in this section as the base case), using a specified concentration for inflow as the boundary type for Biscayne Bay increases the fresh ground-water discharge simulated by the model (fig. 34) by about a factor of two. The average discharge rate for the sensitivity run is 4.9 x 10⁵ m³/d as compared with the average base case value of 2.2 x 10⁵ m³/d.

The effect of Biscayne Bay salinity was evaluated by running a simulation in which the salinity values assigned to the constant-concentration cells in layer 2 were specified using the results from a circulation model of Biscayne Bay, which represents the four years from 1995 to 1998 (John Wang, University of Miami, written commun., 2000). For this sensitivity run, the salinity values from the circulation model were repeated so that a 10-year run could be performed with the regional-scale model and that the time period between the two models would align from 1995 to 1998. The circulation model does not cover the northern part of Biscayne Bay. For this reason, the salinity value assigned to the northern part of Biscayne Bay remains unchanged from the base case. By using salinity values from the circulation model, the simulated value of fresh ground-water discharge to Biscayne Bay increases to about 3.7 x 10⁵ m³/d. A similar run was also performed in which a specified concentration for inflow boundary was used with salinity values from the circulation model. For this run, the average rate of fresh ground-water discharge to Biscayne Bay is 5.8 x 10⁵ m³/d. In both of these runs with the bay salinities set from results of the circulation model, the salinity values assigned to the constant-concentration cells in layer 2 may be less than the salinity of seawater, which means that Biscayne Bay contains a freshwater component. This freshwater component in Biscayne Bay, which is most likely from canal discharge, can be circulated downward into the aquifer. The rate of fresh submarine ground-water discharge determined by this sensitivity run, therefore, contains a freshwater component that results from salinity values in the bay that are less than the salinity of seawater.

The effect of stage in Biscayne Bay was evaluated by running two simulations in which the stage value for each month was increased or decreased by 0.25 m. When the average monthly stage in Biscayne Bay increases by 0.25 m, the simulated ground-water discharge decreases to an average value of 0.7 x 10⁵ m³/d. A decrease in the average monthly stage by 0.25 m substantially increases the ground-water discharge to an average value of 5.4 x 10⁵ m³/d. Results from this sensitivity analysis may be used to infer the effects of sea-level rise on ground-water discharge to Biscayne Bay.

A sensitivity analysis was completed to evaluate the effect of hydraulic conductivity on ground-water discharge to Biscayne Bay by performing two simulations. Hydraulic conductivity was increased by a factor of two for the first run; the simulated fresh ground-water discharge to Biscayne Bay increases to an average value of 3.7 x 10⁵ m³/d. By decreasing the hydraulic conductivity by a factor of two, the simulated ground-water discharge decreases to an average value of 1.3 x 10⁵ m³/d. These data suggest that for the values tested, an error in estimating hydraulic conductivity results in about a one-to-one error in estimates of fresh ground-water discharge to Biscayne Bay.

Results from a sensitivity analysis of aquifer dispersivity suggest that simulated rates of ground-water discharge are sensitive to the value of dispersivity used in the model. When the longitudinal and transverse values for dispersivity are increased by a factor of 10, simulated rates of fresh ground-water discharge decrease (fig. 34). Conversely, when the longitudinal and transverse values for dispersivity are decreased by a factor of 10, simulated rates of fresh ground-water discharge increase (fig. 34).

Two sensitivity runs were performed to evaluate the effect of canal conductance on fresh ground-water discharge to Biscayne Bay. An increase in river conductance by a factor of 10 does not greatly affect the simulated discharge rates (fig. 34), but a decrease in canal conductance does affect the simulated ground-water discharge rates (fig. 34). A decrease in canal conductance increases the freshwater discharge to Biscayne Bay probably because ground water that is intercepted by canals in the base case is allowed to discharge to the bay in the sensitivity run.