Comparison of Seepage Simulation in a Saline Environment below an Estuary Using MODFLOW and SEAWAT

This paper compares the results produced by MODFLOW, a constant-density model, to results produced by SEAWAT, a variable-density model, to investigate the feasibility of using MODFLOW in a saline environment below an estuary known as the Indian River Lagoon. The comparison was conducted over sixteen numerical simulation cases at different conditions of estuarine salinity CL, hydraulic conductivity anisotropy ratio Kr, and water table elevations on the freshwater boundaries in a two-dimensional vertical domain. The use of MODFLOW at the study site under the calibrated Kr distribution ranging from 1000-20,000 was found to accurately match the field-measured and SEAWATsimulated results with a remarkable increase in accuracy at higher groundwater elevations. The study determined a critical value of Kr of 1000 above which, MODFLOW simulations of the variable-density problem produced results that agreed well with those produced by SEAWAT. However, MODFLOW starts to produce significant errors with Kr below the critical value and hence, it should not be used for simulating variable-density environments when Kr<1000. The amount of submarine groundwater discharge (SGD) predicted by either model, and also MODFLOW accuracy in predicting the SGD are directly proportional to the head difference between the groundwater divide elevation and the lagoon water surface, but to a lower extent, are inversely proportional to CL.


Introduction
The presence of high concentrations of dissolved solids in groundwater alters the fluid density and may result in spatial density variations within the flow system. Significant fluid density gradients can substantially affect the groundwater flow patterns introducing thereby mathematical and numerical complexities for simulating such density-dependent flow systems [1]. Examples of such variable-density environments are: saltwater intrusion [1][2][3][4][5][6][7], Submarine Groundwater Discharge (SGD) [8,9], aquifer storage and recovery [10][11][12][13], brine migration [14], coastal wetland hydrology [15], injection of liquid waste in deep saline aquifers [16], and disposal of radioactive waste in salt formations [17,18]. Numerical modeling of variable-density groundwater flow and transport environments such as in saline environments, where the physics of flow and transport are densitydriven, typically relies on the use of variable-density numerical models that incorporate the relationship between fluid density and solute concentration by iteratively solving the flow and transport governing equations. These models are also termed as coupled models. An example of a variable-density coupled model is SEAWAT [19][20][21] which was developed by the U.S. Geological Survey (USGS). SEAWAT couples a modified version of MODFLOW [22][23][24], as a flow simulator for solving the flow equation, and MT3DMS [25] as a solute-transport simulator for solving the mass transport equation [8]. One of the coupling schemes in SEAWAT is termed the implicit approach which is adopted in this paper. An implicit coupling incorporates iterative solution of the flow and transport equations at each time step until the density difference is within a specified value [19]. On the other hand, in situations where spatial density variation is so small that it can be considered to be negligible, a constant-density model such as MODFLOW, also developed by the USGS, may also be used in variabledensity environments to solve the flow equation when the modeler's main objective is to simulate the flow but not the transport conditions.
While it is clearly preferable and generally accepted that a variabledensity coupled flow and transport model best simulates the physical situation in a saline groundwater environment, there may be reasons why a modeler may wish to use a constant-density model in a variabledensity environment. Reasons why modelers would prefer to use a constant-density model such as MODFLOW over a variable-density coupled model such as SEAWAT, even in a saline or variable-density environment, could be: (a) constant-density model simulations require remarkably smaller computation times compared to a variable-density coupled model, and thus, an improved computational efficiency is achieved especially if numerous simulations are required for calibration and validation, b) high level-accuracy simulation results may not always be required in some problems where dropping some important physics to get faster results that are reasonably accurate for making some groundwater management decisions is more practical, c) the ability of the constant-density simulation, like a MODFLOW simulation to be followed by multiple transport codes simulations such as MT3DMS [25] and RT3D [26], in case the model is to be used later for transport simulations, d) familiarity of groundwater modeling community with the constant-density model [19,27], and e) variabledensity coupled models require a larger number of parametric values such as values for molecular diffusion and dispersivity which leads to a higher level of uncertainty.
Comparisons of results obtained by variable-density or coupled models and constant-density or uncoupled models in saline environments have been reported in past several studies. Some of these comparison studies have been conducted in the context of the Henry problem [28] or on laboratory scale physical models [29][30][31][32][33] while other studies used large scale models [6,[34][35][36]. Researchers have found that under some saline conditions, a constant-density model is capable of producing results similar to a variable-density model. For example, Simpson et al. [29] solved the Henry's problem using both variabledensity and constant-density models and found that the predicted location of the 0.5 isochlor, under steady state conditions, were quite similar although the constant-density isochlor did not advance as far into the aquifer as the variable-density isochlor. However, they also found that the matching of the predicted steady state 0.5 isochlors by the two models was much closer when the recharge rate, q, was increased to 2q and 4q. Finally, Simpson et al. [29] noted that the variable-density solution reached the steady state condition much sooner than the constant-density solution indicating that a constantdensity model may not be suitable if transient solutions are required. Simpson et al. [30] found that the discrepancies between results obtained from variable-density and constant-density models were more significant for a modified Henry's problem where the freshwater recharge rate of the standard Henry's problem was halved. Dentz et al. [31] compared variable-density and constant-density analytical and numerical solutions of the Henry problem over a range of two dimensionless groups, the coupling parameter, defined as the ratio of buoyancy flux and viscous forces, and the Péclet number, defined as the ratio of advective and dispersive effects. Dentz et al. [31] concluded that variable-density and constant-density solutions resulted in similar flow patterns when the coupling parameter was far lower than 1 (this corresponds to higher freshwater recharge rate) at moderate Péclet numbers. Goswami et al. [32], in an experimental and numerical analysis, using variable-density and constant-density (i.e., coupled and uncoupled) SEAWAT simulations, of a laboratory-scale porous media tank, made similar conclusions to those of Simpson and Clement [29,30]. Abarca et al. [33] solved the Henry problem in its standard diffusive form and a modified dispersive form for both variabledensity and constant-density conditions. Abarca et al. [33] concluded that, unless the diffusion coefficient of the standard form is reduced by a factor of 10, or the standard form is replaced by the proposed diffusive form, constant-density solution would not result in dramatic changes in the simulated freshwater hydraulic heads or concentration distributions compared to variable-density solution.
Arlai et al. [34] found that both the density-dependent coupled model and the density-independent uncoupled model produced basically similar plume migrations in a Bangkok aquifer. Arlai et al. [34] hypothesized that this was due to groundwater pumping playing a major role in plume migration compared to the role played by density effects. However, in a subsequent paper, Arlai et al. [35] showed that the predicted steady state saline concentrations, in a 1000 by 100 m vertical plane of a saturated coastal aquifer, were much different when using a density-dependent model than a density-independent model, and that the density-dependent model results were closer to the predicted Ghyben-Herzberg interface. Motz et al. [6] conducted multiple numerical experiments on a two-dimensional vertical section of a theoretical coastal aquifer problem with freshwater recharge boundary on the upstream and a seacoast boundary on the downstream using MODFLOW and SEAWAT. Motz et al. [6] found that MODFLOW could closely match the hydraulic heads and fluxes simulated by SEAWAT on the freshwater side of the aquifer when the coastal boundary was represented by specified freshwater hydraulic heads. Subsequently, Motz et al. [36] compared coupled (densitydependent) and uncoupled (density-independent) SEAWAT solutions of saltwater intrusion and seepage circulation on the seacoast boundary of the same system described in Motz et al. [6]. In that comparison study, Motz et al. [36] observed that the density-independent solution produced similar results of saltwater intrusion and seepage circulation to that produced by the density-dependent solution when the ratio of freshwater recharge rate to the density-driven flux was increased. Thus, depending on the boundary conditions and the aquifer properties of a specific saline environment problem, it is reasonable to conclude that there may be situations when a constant-density model can replicate the results of a variable-density coupled model in simulating saline environments although that is not always the case.
The studies discussed above showed through, numerical, analytical, or non-dimensional formulation of variable-density problems, that aquifer anisotropy ratio of hydraulic conductivity K r , and the advective effect of the regional freshwater recharge are amongst the most significant factors controlling the mixed convection ratio (ratio of buoyancy forces to advective forces) on which, density effects depend substantially in saline environments. Sensitivity analyses conducted in this paper also demonstrated the importance of these parameters. Thus, these two parameters, in addition to the salt concentration which is the main source of density-variation, are therefore expected to affect the accuracy of a constant-density solution of variable-density problems where they may increase or dampen the density effects. Also, it is generally known that the density variation has limited effect on horizontal flow and is mostly observed when the flow in a saline environment is vertical, the case that occurs at high K r . Therefore, it is expected that a constant-density solution would resemble a variabledensity solution if there is a relatively small vertical flow component which is likely to happen if K r is large. However, no critical K r value, where the two solutions become similar, has been defined in the literature especially for a real-world saline environment. It is also unknown if there would still be some vertical flow component at that critical K r , i.e., if the two models would produce similar results even if there is a vertical flow component. This paper uses a real-world, two-dimensional transect in the vertical plane below the Indian River Lagoon (IRL), an estuary located on the east-central coast of Florida, to determine if a constant-density model such as MODFLOW and a variable-density model such as SEAWAT, that are using the same calibrated distribution of K r can produce similar results and if they can match 1) measured freshwater hydraulic head contours in the unconfined aquifer below the estuary, 2) groundwater flow directions in the unconfined aquifer, 3) total submarine groundwater discharge (SGD) into the estuary from the adjacent unconfined aquifer, and 4) spatial distribution of SGD below the estuary. The investigation then extends to determine: a) the critical K r below which, results of constant-density and variabledensity models become significantly different, b) if there is still a vertical flow component occur at the critical K r , c) how the IRL salinity C L can affect the amount of SGD into the lagoon and the accuracy of MODFLOW in predicting that amount, and d) how the MODFLOW accuracy improves by increased regional groundwater table elevations on the freshwater boundaries under the calibrated K r and high C L . The constant-density and variable-density modeling results are also compared under both field measured and wide ranges of C L and water table elevations at freshwater boundaries using the calibrated K r values. Thus, this study provides guidance to modelers, regarding values of the anisotropy ratio, at which constant-density models can be used for solving real-world saline environment problems and efficiently overcome the computational burden associated with variable-density models. The paper also provides further guidance by quantifying the error resulting from dropping the physics from the variable-density problem by approximating it as a constant-density problem. It is wellknown that constant-density model simulations take less computation time than variable-density model simulations. However, the authors are not aware of studies that have compared and quantified the stark contrast between the computation times used by the two models in a real-world variable density environment. Thus, the computation time difference has been identified and discussed in this study as well. The MODFLOW version used in this study was MODFLOW-2000 [23]. This study also uses SEAWAT-2000 [20], which couples MODFLOW-2000 and MT3DMS [25] using the implicit finite difference scheme.

Description of Study Area
The study site for this comparison is a transect in the Indian River Lagoon (IRL) which is termed as the Palm Bay transect. The IRL is a coastal estuary that extends over approximately 250 km on the east coast of Florida ( Figure 1). The width of the lagoon varies in the range of 0.8-8.0 km, while its depth ranges from 1 to 3 m. The IRL is connected to the Atlantic Ocean via the Ponce de Leon inlet, Jupiter inlet, Sebastian inlet, Fort Pierce inlet, and the St. Lucie inlet. Everywhere else, it is separated from the ocean by a coastal island strip known as the Barrier Island. The IRL is underlain by an unconfined aquifer which consists mainly of sand and shells with some silt, sandy clay, and clay. The Hawthorn formation, which is a confining impervious layer, underlies the unconfined aquifer and consists mainly of marl and clay [37]. The majority of groundwater seepage into the IRL comes from the watershed on the west side known as the Mainland although some seepage also comes from the Barrier Island.
As shown in Figure 2, the Palm Bay transect is oriented perpendicular to the IRL shoreline, and extends from the groundwater divide on the Mainland to the Atlantic Ocean passing through the IRL and the Barrier Island. As a result of the mixing of salt water from the ocean through the inlets, and freshwater received through rainfall, groundwater, non-point runoff and point runoff through canals and rivers draining an approximately 3575 km 2 [38], the IRL waters are generally brackish and have a salinity in the annual range of about 10.8 to 37.8 g/L at the study site. This corresponds to a normalized salinity range of 0.3-1.05 based on ocean salinity of 36 g/L.

Data Collection
Initially, the horizontal domain of the transect was determined by locating the groundwater divide ( Figure 2) using multiple monitoring wells installed on the Mainland. Relevant field data is measured to determine appropriate boundary conditions, and for the calibration of the numerical models. Water table elevations were measured in several wells installed on the Mainland and Barrier Island to obtain water surface profiles which were used as boundary conditions in the numerical models. The lagoon bed profile at the transect location shown in Figure 3 was determined by traversing the lagoon cross section by a boat and measuring the bed depth under the water surface at approximately 50 locations. Piezometric heads and groundwater salinity are simultaneously measured at eight stations situated across the IRL transect as shown in Figures 2 and 3. Three of these stations are located on the lagoon shores and are termed WS1, WS2, and ES, while five stations are somewhat uniformly spaced along the lagoon and are labelled S1 to S5. Nested shallow and deep piezometers were installed at each of the eight lagoon stations as depicted in Figure 3. The depth of the shallow locations ranges from 0.5 to 1.5 m while that of the deep locations ranges from 1.5 to 6.1 m below the IRL bed. Piezometers were made of 1.9 cm diameter PVC pipes. A 30 cm long, 1.9 cm diameter well screen made of slotted PVC pipe with #10 slot size, was attached to the end of each piezometer. Mid-screen locations range from (0.5-4.31 m) and (1.72-6.82 m) below the National Geodetic Vertical Datum of 1929 (NGVD 29) for the shallow and deep piezometers, respectively.
Piezometers were driven into the aquifer using a jetting technique in which, a 3.175 cm diameter PVC casing was jetted into the sediment to the desired depth by a 1.49 Kw centrifugal pump. Once the desired depth is reached, a piezometer was inserted into the casing. The casing was then pulled out of the sediment slowly and sand was then packed around the screen and then the annular space around the piezometer was sealed by bentonite clay chips. The depth of the clay seal from the lagoon bed was about 180 cm and 30 cm for the deep and the shallow piezometers, respectively. The water that was introduced into the wells during the jetting phase was removed by wells development. In order to protect the piezometers from potential sabotage, shallow piezometers were terminated approximately 5 cm above the lagoon bed while the deep piezometers were terminated within 60 cm from the lagoon water surface. The top of each piezometer was threaded and fitted with an O-ring seal where a cap was screwed to create a leak-proof seal. Whenever measurements were taken, the cap was removed and a PVC extension pipe was screwed into the top of the piezometer so that it extends above the water surface.
Groundwater salinity was measured by extracting samples from the deep and shallow piezometers, and reading the salinity with a YSI salinity meter, Model 85. The salinity meter converts conductivity readings to salinity based on ASTM algorithms found in ASTM Designation D1125-82, and has a measurement range of 0 to 80 parts per thousands (ppt), a resolution of 0.1 ppt, and an accuracy of 0.5 percent. Lagoon water salinity was also measured at each of eight lagoon stations at the same time as the other measurements. The lagoon   water surface elevation from the NGVD 29 was established during each sampling event based on three bench marks installed in the study area.
The horizontal hydraulic conductivity K h at the stations S1 through S5 in addition to WS1 and ES, at both shallow and deep depths, was estimated in the laboratory by obtaining grain size distribution curves and using the Hazen's equation, K=C(d 10 ) 2 , where d10 is the effective diameter, and C is an empirical constant=0.01. Soil samples used for the K h analyses were extracted at these stations from shallow locations ranging in depth of 0 to 0.5 m and deep locations ranging in depth of 1.4 to 1.9 m. An average K h of about 30 m/day is obtained for both deep and shallow locations. A sampling event on any given day, comprised of measuring the water table elevations in the Mainland and Barrier Island wells, the piezometric heads and groundwater salinity in the offshore and onshore piezometers, and the IRL water surface elevation and salinity. This study uses field data of three sampling events conducted on May, August, and September. Table 1 gives the measured elevations of water table at the locations of the groundwater divides on the mainland and the Barrier Island as well as the measured elevations of the IRL water surface for the three sampling events. Values of the head difference ΔH between the groundwater divide and the IRL surface are also shown in Table 1.

Numerical Modeling
The details regarding the MODFLOW and SEAWAT numerical modeling set up including models domain and discretization, initial and boundary conditions, and calibration and aquifer parameters are presented subsequently.

Models domain and discretization
The model domain for both MODFLOW and SEAWAT is a twodimensional vertical cross-section extending horizontally over a distance of 9,740 m from the groundwater divide at the Mainland to the Atlantic Ocean ( Figure 4). The left-hand boundary of the domain is the groundwater divide at the Mainland, while the right-hand boundary is the Atlantic Ocean. The width of the IRL at this transect is 2.61 km. The groundwater divide on the Mainland is located at 6.23 km from the west shore of the IRL, while the Atlantic Ocean is 0.9 km from the east shore. The bottom boundary is the top of the Hawthorn Formation, which is assumed to be horizontal, impermeable, and is at a depth of approximately 33.5 m at the Palm Bay transect location. The model domain was discretized into one row, 76 columns and 22 layers ( Figure 4). The single row of the finite difference mesh was arbitrarily set with a width of 1 m. The columns spacing ranged from 15 to 300 m, while layers spacing ranged from 0.3 to 6 m. The exact finite difference mesh shown in Figure 4 is used identically in both MODFLOW and SEAWAT models.

Initial and boundary conditions
Initial condition for the hydraulic heads is 33.5 m (i.e., the total depth of the model domain) at all internal nodes for both MODFLOW and SEAWAT, while the saltwater concentrations are set to zero at all internal nodes for SEAWAT. The boundary conditions are described in the following section.

MODFLOW boundary conditions:
When setting up the boundary conditions for MODFLOW, the mainland and Barrier Island freshwater boundaries (AB and CD) shown in Figure 4 were modeled as constant head boundaries (Dirichlet boundaries). The constant head values of AB and CD boundaries were calculated from the functions f 1 (x) and f 2 (x) obtained from the field measured water table elevations in the Mainland and Barrier Island, respectively. Each of these functions is a polynomial equation of a trend line obtained from statistical regression between the head values measured at each monitoring well and its distance from the lagoon shoreline. Both of the Hawthorn formation (EF) and the groundwater divide (FA) boundaries were modeled as no flow boundaries (Neumann boundaries). Motz et al. [6] found that specifying equivalent freshwater heads at the brackish water boundaries in a constant-density model yielded results closest to those obtained by a variable-density coupled model. Representing the boundary conditions at the saline boundaries in constant-density models by specified freshwater hydraulic heads was also used by Simpson and Clement [29][30][31]36]. Therefore, in MODFLOW, the hydraulic heads at the lagoon bed boundary (BC) and the ocean boundary (DE) were converted to freshwater hydraulic heads using equations 2 and 4, respectively, and these boundaries were also treated as constant head boundaries. This is equivalent to coupling the flow and transport only at the saltwater boundaries (i.e., the lagoon and the ocean) while ignoring the coupling everywhere else inside the model domain. The boundary conditions used in MODFLOW can be mathematically described in equations 1 to 6 below: Boundary BC: ( ) Boundary CD: h=f 2 (x) Boundary DE: ( ) Boundary EF: ∂h/∂z=0 Boundary FA: ∂h/∂x=0 (6) where h is the hydraulic head specified in the boundary conditions, h f is the equivalent freshwater hydraulic head, f 1 (x) and f 2 (x) are functions obtained by respectively measuring water table elevations on the Mainland and Barrier Island, z is the elevation of the nodes from the top of the Hawthorn Formation, C L and C S are the measured concentration of the lagoon and seawater in the form of normalized salinity, respectively, ρ s is the saltwater density, ρ f is the freshwater density, and h S is the height of the piezometric surface above the Hawthorn formation.
SEAWAT boundary conditions: SEAWAT model requires the boundary conditions for both hydraulic heads and saltwater concentrations. The hydraulic head boundary conditions used for SEAWAT are identical to those used for MODFLOW in equations 1-6. However, the conversion into freshwater hydraulic heads at the IRL and Ocean boundaries BC and DE in equations 2 and 4, respectively, is done internally in SEAWAT using the specified concentration at each boundary [21]. The saltwater concentration boundary conditions for SEAWAT can be mathematically described as: Boundary BC: C=C L Boundary CD: C=0 Boundary DE: C=C S Boundary EF: ∂C/∂z=0 (11) Boundary FA: ∂C/∂x=0 (12)

Model calibration and aquifer parameters
The vertical hydraulic conductivity of the aquifer below the estuary at the study transect was calibrated using both statistical and visual methods. In the statistical method, the model-predicted and fieldmeasured freshwater hydraulic heads nodal values at the shallow and deep locations shown in Figure 3 were compared by estimating the Root-Men Squared Error (RMSE), Nash-Sutcliffe Efficiency (NSE) index, and testing the null hypothesis using the two-sided test. In the visual comparison method, the calibration was conducted by visually matching the predicted and measured freshwater head equipotential lines below the study transect. Model calibration resulted in a RMSE of 0.05 m, a NSE of 0.96, a two-sided t-test not rejecting the null hypothesis that the means of the measured and predicted values are equal, and a good visual comparison of the measured and predicted freshwater head equipotential lines.
The calibrated values of the aquifer vertical hydraulic conductivity ranged from 0.0015 to 0.03 m/day with a predominant value of   0.015 m/day. The value of the horizontal hydraulic conductivity was predominantly 30 m/day. This range of hydraulic conductivity values leads to K r ranging from 1000-20,000 with a predominant value of 2000. The K r is defined in this paper as the ratio of horizontal and vertical hydraulic conductivities (K h /K v ). The input data and aquifer properties used in the MODFLOW and SEAWAT models in this study are listed in Table 2. The K h was estimated from Hazen's equation as described previously. The K v and lagoon salinity are variable as discussed in the study. Freshwater and seawater normalized salinities and densities are constants. Multiple dispersivity values were investigated and found not to affect the predicted seepage values.

Results and Discussion
Results produced by MODFLOW and SEAWAT models were compared for sixteen cases described in Table 3. In presenting the SGD results of these cases, the term, relative error, implies the difference between the two sets of model results with the SEAWAT results assumed to be the "true" results. The results of Cases 1 to 3 compare the outputs of both models using the calibrated distribution of K r ranging from 1000-20,000 at different boundary conditions measured during field sampling events conducted on May, August, and September. The respective lagoon normalized salinity C L in the three cases was 0.844, 0.583 and 0.306 and the measured groundwater elevations  Table 1. In the Cases 4 through 8 (Table 3), C L was arbitrarily increased from 0.9 to 1.05 to determine how MODFLOW results compare with SEAWAT results under higher IRL salinity conditions. Cases 4 through 8 utilize the same Mainland and Barrier Island boundary conditions and the same IRL water surface elevation measured on May and used in Case 1. Cases 4 through 8 also utilize the calibrated K r distribution, hence the aquifer K r for these cases also varies from 1000-20,000 as in Cases 1 to 3. Cases 9 through 13 tested the accuracy of MODFLOW at decreasing K r under relatively high saline condition as the lagoon salinity was kept at 0.9. In these cases, the K r was changed from 100,000 to 10 while the boundary conditions and lagoon salinity were kept identical to those used for Case 4. Case 14 is identical to Case 12 with the exception that the C L was reduced from 0.9 to 0.3. This numerical experiment was conducted to determine the accuracy of MODFLOW at low anisotropy ratio and at low lagoon salinity. Cases 15 and 16 are also identical to Case 1 with the exception that the water table elevations used for the boundary conditions of the Mainland and the Barrier Island cells were increased by 5 percent in Case 15 and 10 percent in Case 16 to examine the effect of higher water table elevation. In Case 15, the water table elevations at the water table divides on the Mainland and the Barrier Island were 2.05 m and 1.7 m, respectively higher than in Case 1. Also, in Case 16, the respective elevations on the Mainland and the Barrier Island were 4.1 m and 3.4 m higher than in Case 1.

Cases No. 1 to 3
The following results are compared for Cases 1 to 3: a) freshwater hydraulic head distributions below the IRL, b) flow directions in the form of velocity vectors, c) total SGD into the IRL from the underlying aquifer, and d) spatial distribution of SGD below the IRL. Figures 5-7 compare the calibrated MODFLOW and SEAWAT model-simulated freshwater hydraulic head distributions in the entire modeling domain for Cases 1, 2, and 3, respectively. Even though the shape and magnitude of the equipotential freshwater hydraulic heads are fairly different for the three cases, both MODFLOW and SEAWAT predicted almost identical distributions. Inspecting the comparisons at the two saline boundaries of the models, i.e., the lagoon and Ocean boundaries, (Figures 5-7), it can be seen that the equipotential lines predicted by the two models show some discrepancy on the bottom half portion of the sea boundary where salt intrusion is predominant, while they seem to be fairly similar on the top portion. On the lagoon boundary, the two distributions are very similar for all three cases except that MODFLOW is always predicting the contours locations a little higher than their actual locations predicted by SEAWAT. However, this discrepancy in the contour locations below the estuary seems to become very minor for Case 3 (Figure 7), where the IRL salinity is the lowest and the groundwater elevation on the Mainland (Table 1) is the highest, compared to the first two cases. In general, the variable-density effects predominating the sea side of the model, do not seem to have significant effect on the accuracy of predicting the equipotential lines below the estuary especially in Case 3.
A comparison of the model predicted velocity vectors by MODFLOW and SEAWAT in Case 1 is shown in Figure 8. In general, both models predicted very similar flow patterns in that: a) the groundwater flows upward into the IRL, and b) there is recirculation of the ocean water from the lower part of the aquifer adjacent to the ocean boundary back into the ocean at the upper region of the aquifer. The meteoric groundwater discharge (MGWD) originating from the Barrier Island splits into two directions with a portion going to the

Density of seawater, ρ s (kg/m 3 ) 1025
Density of freshwater, ρ f (kg/m 3 ) 1000    Note: SGD in m 3 /day/m is the total SGD in m 3 /day per m of transect width; (a) and (b) imply that the groundwater elevations collected on May were increased by 5% and 10% in Cases 15 and 16, respectively.  IRL and the rest going to the Atlantic Ocean. Both models predict that the flow direction below the lagoon is mostly oriented upward towards the lagoon even though the vertical hydraulic conductivity K v is significantly small. The simulated results of Cases 2 and 3 are not shown as they were very similar to the patterns of groundwater flow shown in Figure 8.
The values of the predicted total SGD into the IRL transect, over the three cases, are shown in Table 3, and the relative error in the MODFLOW results ranged from 3.2 to 9.4 percent. In all these cases, MODFLOW predicted slightly higher SGD with the highest error (9.4 percent) occurring in Case 1 when the IRL salinity was relatively higher and the groundwater elevation on the Mainland (Table 1) was relatively lower. Regardless of the difference in the predicted SGD, both models predicted that the SGD into the IRL is the highest on August (Case 2) followed by September (Case 3) then May (Case 1). However, referring to the measured groundwater elevations on the Mainland, it can be noticed that those elevations were the highest on September followed by August and then on May. Therefore, one may expect that the highest SGD into the IRL would follow that sequence. However, the measured IRL water surface elevations shown in Table 1, seem to also increase with increasing Mainland groundwater elevation. Therefore, higher mounts of rainfall on the boundaries do not necessarily produce a significantly higher SGD since the rainfall not only increases the water table elevation but also elevates the water level in the estuary thereby maintaining the hydraulic gradient. Infact, the predicted SGD values seem to be proportional to the head differences ΔH between the groundwater divide on the mainland and the IRL water surface (Table 1).
A comparison of the predicted spatial distribution of the SGD, from the west shore to the east shore of the IRL transect for the first three cases, is illustrated in Figures 9-11, respectively. It can be seen that the patterns predicted by MODFLOW and SEAWAT are nearly identical. Both models predict higher seepage on the west shore side as most groundwater seepage comes from of the Mainland. Figures 9-11 don't only show that both constant-density and variable-density models agree very well, but also show that the SGD is not just a near-shore phenomenon but extends all the way up to 1600 m from the shoreline. This result indicates that the SGD into the IRL can occur at much greater distances away from the shoreline than predicted by previous studies (Martin et al.). Also, there can be high and low SGD producing zones within a few meters of each other due to a sudden variation in the vertical hydraulic conductivity K v values. For example, there is a sharp decline, followed by a sudden rise, in the SGD flux within 400 m from the west shore. This sharp variation in SGD into the IRL can be missed in seepage meter studies unless seepage meters are placed within 10 m (or even closer) of each other. Although Cases 1 to 3 share the same K r distribution, the IRL salinity and the groundwater elevations on the boundaries and the IRL water surface elevations were different. Therefore, it cannot be decided at this point if the results of these cases presented in Table 3 and Figures 5-11 were mostly affected by the changing C L , the boundary conditions, or both. Therefore, in the next phase (Cases 4-8) of this investigation, comparisons were accomplished at a changing C L while keeping the K r and the boundary conditions constant as in Case 1.

Cases No. 4 to 8
The SGD values predicted by both MODFLOW and SEAWAT for Cases 4 through 8 are shown in Table 3 and these results in addition to the results of Case 1 indicate that the accuracy of MODFLOW in predicting the SGD decreases slightly with increasing IRL salinity.
As a result, MODFLOW accuracy in predicting the SGD spatial distribution, which depends directly on the SGD value, would also decrease slightly with increasing salinity. The results also indicate that the highest relative error (11.6%) occurs when the normalized lagoon salinity has its highest value of 1.05 (Case 8). The amount of SGD into the IRL predicted by both models seems to be inversely proportional to the lagoonal water salinity. However, changing the IRL salinity for these cases seems not to change the amount of seepage as much as noticed between Cases 1 to 3. It can be concluded that the amount of SGD is primarily affected by the hydraulic gradient present between the IRL and the groundwater levels on the boundaries and is affected, to a lot lower level, by the IRL salinity. Although not presented here, the equipotential freshwater hydraulic head distributions predicted by MOFLOW and SEAWAT for Cases 4 to 8 did not change significantly with increasing IRL salinity and were very similar to those shown in Figure 5 for Case 1. These results also indicate that MODFLOWpredicted freshwater hydraulic head distributions stay similar to those predicted by SEAWAT even if the IRL salinity increases significantly although the error in predicting SGD increases slightly. Therefore, it can be concluded that the shape and magnitude of the freshwater hydraulic head distributions on one hand, and the agreement of those distributions in shape and magnitude between the two models depend primarily on the Mainland groundwater divide elevation and not the IRL salinity for the cases having the same K r distribution. This explains why MODLOW-predicted freshwater hydraulic head distribution for Case 3 (Figure 7), which has the highest groundwater divide elevation at the Mainland (Table 1), was more similar to those of SEAWAT compared to Figures 5 and 6.

Cases No. 9 to 13
It can be seen from the SGD values shown in Table 3 that the relative error is less than 10% when K r is 10,000 or higher (Cases 9 and 10) and less than 15% where K r is 1000 or higher. However, when K r is lowered below 1000, the relative error becomes significantly higher and is as high as 129.4 percent when K r is equal to 10 (Case 13). The large difference in the predicted SGD values is also reflected in the freshwater hydraulic head contours predicted by the two models as shown in Figure 12 which compares the results for Case 12. The comparison of the freshwater hydraulic head contours of Case 13, which has the highest relative error of 129.2 percent, although not shown here, was even worse than that shown in Figure 12. At those low K r values the predicted flow directions below the lagoon were mostly vertical either upward or downward. Downward vertical flow gives rise to the saline estuarine water overlaying the aquifer with a salinity as high as 0.9 to penetrate deeper into the aquifer while the upward vertical flow drifts more seawater up into the aquifer. Those vertical flow patterns make the variable-density effects below the estuary to be predominant over the advective effects causing the constant-density MODFLOW model to fail at this point. Inspecting Figure 13 that presents the SGD results of Cases 9 to 13 shows that a K r of 1000 (Case 11) looks to be the critical value below which, the use of the constant-density MODFLOW model results in significant loss of accuracy. A comparison of flow directions in a selected portion of the model domain extending horizontally from the west shore to the east shore of the lagoon and vertically from the lagoon surface down to an arbitrarily selected depth of -9 m NGVD 29 is shown in Figure 14. It can be seen from Figure 14 that MODFLOW is still predicting very similar flow directions to those produced by SEAWAT at the critical K r of 1000 although there is still a vertical flow component below the lagoon.

Case No. 14
Case 14 is a repetition of Case 12 with the exception that the C L is 0.3 instead of 0.9 while K r was kept at 100. Results of this case show that although reducing the IRL salinity from 0.9 to a value closer to freshwater salinity resulted in reducing the relative error in SGD prediction from 38.2 to 11.2 percent (Table 3), MODFLOW fails to produce accurate freshwater hydraulic head distribution below the lagoon at such a low K r (Figure 15). It is obvious that modeling this system based on a constant-density assumption under low K r conditions regardless of C L value is not accurate at all.

Cases No. 15 and 16
These cases assess the degree of improvement in the MODFLOW accuracy in predicting the SGD values and the freshwater hydraulic head equipotential lines if Case 1 was re-run with the same K r and C L but with higher groundwater elevations on the Mainland and Barrier Island boundaries. The freshwater hydraulic head equipotential lines predicted by the two models for Cases 15 and 16 are shown in Figures  16 and 17, respectively. These figures as well as the SGD values shown in Table 3 indicated that the accuracy of MODFLOW improves substantially if the water table elevations are higher. Increasing the water table elevations by 5% in Case 15 and 10% in Case 16 resulted in reducing the error in estimating the SGD from 9.4% in Case 1 to 4.5% and 3.3% in Cases 15 and 16, respectively. The drop of the relative error from 9.4% into 4.5% and 3.3% is equivalent to an increase in MODFLOW accuracy by 52% and 65%, respectively.

Model computation time
The computation times for all cases described in Table 3 were approximately 0.07 seconds for MODFLOW and approximately 564 seconds for SEAWAT. Analysis conducted by refining the finite difference mesh to 258 columns and 110 layers compared to 76 columns and 22 layers in the original model showed that the discrepancy between the computation times for the two models increased even more than that. MODFLOW required approximately 1.5 seconds while SEAWAT required approximately 12,600 seconds (3.5 hours) for simulating the results of Case 1. The computation times reported here are clock times of PC with an Intel Core i3-2310M 2.10 GHz CPU. These computation times indicate that, depending on the mesh size, a constant-density model may be faster by a factor of 8000 or more than a coupled model in solving identical problems. This difference in computation times is likely to become even higher for complex three-dimensional (3-D) models.

Summary and Conclusions
This paper investigates the accuracy of results produced by MODFLOW, a constant-density model, in a saline environment below an estuary known as the Indian River Lagoon (IRL), by comparing its results with those of SEAWAT, which is a variable-density coupled model. The results that were compared included: 1) measured freshwater hydraulic head contours in the unconfined aquifer below the estuary, 2) groundwater flow directions in the unconfined aquifer, 3) total submarine groundwater discharge (SGD) into the estuary from the adjacent unconfined aquifer, and 4) spatial distribution of SGD below the estuary. The comparison was conducted over sixteen numerical experiments at different conditions of estuarine salinity C L , hydraulic conductivity anisotropy ratio K r , and water table elevations on the freshwater boundaries.
The results presented in this paper showed that: a) the use of MODFLOW for modeling the IRL at the study site under its calibrated K r range of 1000-20,000 was satisfactory and accurate to within approximately 3 to 9% regardless of the IRL salinity and groundwater elevations on the boundaries with an increase in its accuracy by about 52% and 65% by increasing the measured groundwater elevations by 5% and 10%, respectively, b) results produced by MODFLOW can be in close agreement with those obtained by SEAWAT if K r is greater than a critical value of 1000 regardless of the lagoon salinity, the conditions   under which, MODFLOW produced results within less than 15% to those predicted by SEAWAT, c) MODFLOW should probably not be used in saline environments if K r is less than 1000 under any conditions even when lagoon salinity is low, d) there is still vertical flow component predominating below the lagoon even at the critical K r of 1000, e) the amount of SGD predicted by either model and also the MODFLOW accuracy in predicting the SGD are directly proportional to the head difference between the groundwater divide elevation and the lagoon water surface, but, to a lower extent, are inversely proportional to C L , f) both MODFLOW and SEAWAT predicted that the SGD can occur at much greater distances away from the shoreline than predicted by previous studies, g) both MODFLOW and SEAWAT agreed very well in showing high and low SGD producing zones across the transect and both models showed that depending on the K v value, those zones can occur within a few meters of each other and may be missed by seepage meter studies, and h) for the two meshes used in the analyses, MODFLOW was faster than SEAWAT by a factor of greater than 8000 and this discrepancy in computation times becomes even more significant as the mesh is refined.
In summary, under certain conditions, constant-density models such as MODFLOW could be a viable option for modeling saline environments, particularly estuarine environments, if the primary reason for the modeling is to determine flow rates and not saltwater transport, as there is a vast discrepancy in the computation times needed by the two models to solve identical problems. Of course, it is advisable that modelers compare results for their boundary value problem with both variable-density and constant-density models prior to deciding if the constant-density model is feasible or not.