Comparison of Seepage Simulation in a Saline Environment below an Estuary Using MODFLOW and SEAWAT
Received Date: Apr 03, 2017 / Accepted Date: Apr 12, 2017 / Published Date: Apr 18, 2017
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.
Keywords: Coastal aquifer; MODFLOW; SEAWAT; Indian River lagoon; Estuary; Submarine groundwater discharge
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 . Examples of such variable- density environments are: saltwater intrusion [1-7], Submarine Groundwater Discharge (SGD) [8,9], aquifer storage and recovery [10-13], brine migration , coastal wetland hydrology , injection of liquid waste in deep saline aquifers , 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-21] which was developed by the U.S. Geological Survey (USGS). SEAWAT couples a modified version of MODFLOW [22-24], as a flow simulator for solving the flow equation, and MT3DMS  as a solute-transport simulator for solving the mass transport equation . 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 . 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  and RT3D , 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  or on laboratory scale physical models [29-33] while other studies used large scale models [6,34-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.  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.  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.  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.  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.  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. , 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.  solved the Henry problem in its standard diffusive form and a modified dispersive form for both variabledensity and constant-density conditions. Abarca et al.  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.  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.  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.  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.  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.  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.  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. . In that comparison study, Motz et al.  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 Kr, 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 Kr. 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 Kr is large. However, no critical Kr 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 Kr, 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 Kr 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 Kr 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 Kr, c) how the IRL salinity CL 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 Kr and high CL. The constant-density and variable-density modeling results are also compared under both field measured and wide ranges of CL and water table elevations at freshwater boundaries using the calibrated Kr 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 . This study also uses SEAWAT-2000 , which couples MODFLOW-2000 and MT3DMS  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 . 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 km2 , 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.
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 Kh 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(d10)2, where d10 is the effective diameter, and C is an empirical constant=0.01. Soil samples used for the Kh 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 Kh 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.
|Sampling event||Groundwater divide elevation (m above NGVD29)||IRL water surface elevation (m above NGVD29)||ΔH (m)|
Table 1: Measured Elevations of Groundwater Divides on Mainland and Barrier Island and Water Surface of IRL at Different Sampling Events.
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 f1(x) and f2(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.  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-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 AB: h=f1(x) (1)
Boundary CD: h=f2(x) (3)
Boundary EF: ∂h/∂z=0 (5)
Boundary FA: ∂h/∂x=0 (6)
where h is the hydraulic head specified in the boundary conditions, hf is the equivalent freshwater hydraulic head, f1(x) and f2(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, CL and CS 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 hS 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 . The saltwater concentration boundary conditions for SEAWAT can be mathematically described as:
Boundary AB: C=0 (7)
Boundary BC: C=CL (8)
Boundary CD: C=0 (9)
Boundary DE: C=CS (10)
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 Kr ranging from 1000-20,000 with a predominant value of 2000. The Kr is defined in this paper as the ratio of horizontal and vertical hydraulic conductivities (Kh/Kv). The input data and aquifer properties used in the MODFLOW and SEAWAT models in this study are listed in Table 2. The Kh was estimated from Hazen’s equation as described previously. The Kv 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.
|Data or Aquifer Property||Value|
|Horizontal hydraulic conductivity, Kh (m/day)||30|
|Vertical hydraulic conductivity, KV (m/day)||variable (see Table 3)|
|Specific Storage, Ss (m-1)||0.00001|
|Specific yield, Sy||0.01|
|Longitudinal dispersivity, αL (m)||30|
|Transverse dispersivity, αT (m)||3|
|Diffusion coefficient, Dm (m2/d)||0|
|Seawater normalized salinity, Cs||1|
|Freshwater normalized salinity, Cf||0|
|Lagoon water normalized salinity, CL||variable (see Table 3)|
|Density of seawater, ρs (kg/m3)||1025|
|Density of freshwater, ρf (kg/m3)||1000|
Table 2: Input Data for Models.
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 Kr 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 CL in the three cases was 0.844, 0.583 and 0.306 and the measured groundwater elevations on the boundaries in addition to the IRL water surface elevations are presented in Table 1. In the Cases 4 through 8 (Table 3), CL 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 Kr distribution, hence the aquifer Kr 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 Kr under relatively high saline condition as the lagoon salinity was kept at 0.9. In these cases, the Kr 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 CL 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.
|Case no.||Anisotropy ratio Kr||IRL salinity CL||Boundary conditions||Total SGD into IRL|
|SEAWAT (m3/day/m)||MODFLOW (m3/day/m)||Relative error (%)|
|1||1000-20,000||0.844||May||1.790 × 10-4||1.960 × 10-4||9.4|
|2||1000-20,000||0.583||August||2.170 × 10-4||2.250 × 10-4||3.9|
|3||1000-20,000||0.306||September||1.820 × 10-4||1.880 × 10-4||3.2|
|4||1000-20,000||0.9||May||1.777 × 10-4||1.957 × 10-4||10.1|
|5||1000-20,000||0.95||May||1.764 × 10-4||1.953 × 10-4||10.7|
|6||1000-20,000||1||May||1.751 × 10-4||1.949 × 10-4||11.3|
|7||1000-20,000||1.025||May||1.747 × 10-4||1.947 × 10-4||11.5|
|8||1000-20,000||1.05||May||1.744 × 10-4||1.946 × 10-4||11.6|
|9||1,00,000||0.9||May||4.120 × 10-5||4.340 × 10-5||5.3|
|10||10,000||0.9||May||1.070 × 10-4||1.170 × 10-4||9.1|
|11||1,000||0.9||May||2.080 × 10-4||2.380 × 10-4||14.6|
|12||100||0.9||May||2.590 × 10-4||3.570 × 10-4||38.2|
|13||10||0.9||May||2.190 × 10-4||5.020 × 10-4||129.4|
|14||100||0.3||May||3.170 × 10-4||3.530 × 10-4||11.2|
|15||1000-20,000||0.844||(a)||4.730 × 10-4||4.950 × 10-4||4.5|
|16||1000-20,000||0.844||(b)||7.680 × 10-4||7.930 × 10-4||3.3|
Note: SGD in m3/day/m is the total SGD in m3/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.
Table 3: Description and Results of Numerical Modeling Cases.
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 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 Kv 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 Kv 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 Kr 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 CL, the boundary conditions, or both. Therefore, in the next phase (Cases 4-8) of this investigation, comparisons were accomplished at a changing CL while keeping the Kr 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 Kr 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 Kr is 10,000 or higher (Cases 9 and 10) and less than 15% where Kr is 1000 or higher. However, when Kr is lowered below 1000, the relative error becomes significantly higher and is as high as 129.4 percent when Kr 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 Kr 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 Kr 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 Kr 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 CL is 0.3 instead of 0.9 while Kr 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 Kr (Figure 15). It is obvious that modeling this system based on a constant- density assumption under low Kr conditions regardless of CL 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 Kr and CL 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 CL, hydraulic conductivity anisotropy ratio Kr, 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 Kr 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 Kr 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 Kr 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 Kr 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 CL, 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 Kv 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.
We greatly appreciate the support provided by the St. Johns River Water Management District (SJRWMD) for funding the work of this research.
- Paniconi C, Khlaifi I, Lecca G, Giacomelli A, Tarhouni J (2001) A modelling study of seawater intrusion in the Korba Coastal Plain, Tunisia. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere 26: 345-351.
- Shoemaker WB, Edwards KM (2003) Potential for saltwater intrusion into the lower Tamiami aquifer near Bonita Springs, southwestern Florida. Water Resources Investigations Report 03- 4262, U.S. Geological Survey, Tallahassee, FL.
- Dausman A, Langevin CD (2005) Movement of the saltwater interface in the surficial aquifer system in response to hydrologic stresses and water-management practices, Broward County, Florida. Scientific Investigations Report 2004-5256, U.S. Geological Survey, Reston, Virginia.
- Lin J, Snodsmith JB, Zheng C, Wu J (2009) A modeling study of seawater intrusion in Alabama Gulf Coast, USA. Environ Geol 57: 119-130.
- El-Bihery M (2009) Groundwater flow modeling of quaternary aquifer Ras Sudr, Egypt. Environ Geol 58: 1095-1105.
- Motz L, Sedighi A (2009) Representing the coastal boundary condition in regional groundwater flow models. J Hydrol Eng 14: 821-831.
- Ding F, Yamashita T, Lee HS, Pan J (2014) A modelling study of seawater intrusion in the liao dong bay coastal plain, china. Journal of Marine Science and Technology (Taiwan) 22: 103-115.
- Langevin CD (2003) Simulation of submarine ground water discharge to a marine estuary: Biscayne Bay, Florida. Ground Water 41: 758-771.
- Li X, Hu BX, Burnett WC, Santos IR, Chanton JP (2009) Submarine ground water discharge driven by tidal pumping in a heterogeneous aquifer. Ground Water 47: 558-568.
- Ward JD, Simmons CT, Dillon PJ (2007) A theoretical analysis of mixed convection in aquifer storage and recovery: How important are density effects. J Hydrol 343: 169-186.
- Ward JD, Simmons CT, Dillon PJ, Pavelic P (2009) Integrated assessment of lateral flow, density effects and dispersion in aquifer storage and recovery. J Hydrol 370: 83-99.
- Minsley BJ, Ajo-Franklin J, Mukhopadhyay A, Morgan FD (2011) Hydrogeophysical methods for analyzing aquifer storage and recovery systems. Ground Water 49: 250-269.
- Zuurbier KG, Zaadnoordijk WJ, Stuyfzand PJ (2014) How multiple partially penetrating wells improve the freshwater recovery of coastal aquifer storage and recovery (ASR) systems: A field and modeling study. J Hydrol 509: 430-441.
- Shafer JM, Brantley DT, Waddell MG (2010) Variable-density flow and transport simulation of wellbore brine displacement. Ground Water 48: 122-130.
- Langevin C, Swain E, Wolfert M (2005) Simulation of integrated surface-water/ground-water flow and salinity for a coastal wetland and adjacent estuary. J Hydrol 314: 212-234.
- Baidariko EA, Pozdniakov SP (2011) Simulation of liquid waste buoyancy in a deep heterogeneous aquifer. Water Resources 38: 972-981.
- Rumynin VG, Mironenko VA, Sindalovsky LN, Boronina AV, Konosavsky PK, et al. (2000) Conceptual and numerical modelling of density induced migration of radioactive contaminants at the Lake Karachai waste disposal site. Proceedings of the ModelCARE'99 Conference, September 20, 1999 - September 23, IAHS Press, Zurich, Switz, pp: 405-411.
- Fraser Harris AP, McDermott CI, Kolditz O, Haszeldine RS (2015) Modelling groundwater flow changes due to thermal effects of radioactive waste disposal at a hypothetical repository site near Sellafield, UK. Environmental Earth Sciences 74: 1589-1602.
- Guo W, Langevin CD (2002) User's Guide to SEAWAT: A Computer Program for Simulation of Three-Dimensional Variable-Density Ground-Water Flow. Techniques of Water-Resources Investigations 6-A7, U.S. Geological Survey, Tallahassee, FL.
- Langevin CD, Shoemaker WB, Guo W (2003) MODFLOW-2000, the U.S. Geological Survey modular ground-water model: Documentation of the SEAWAT-2000 version with the variable- density flow processes (VDF) and the integrated MT3DMS Transport Processes (IMT). USGS Open-File Rep. 03-426, U.S. Geological Survey, Tallahassee, FL.
- Langevin CD, Guo W (2006) MODFLOW/MT3DMS-based simulation of variable density ground water flow and transport. Ground Water 44: 339–351.
- McDonald MG, Harbaugh AW (1988) A Modular Three-Dimensional Finite-Difference Ground- Water Flow Model. Techniques of Water-Resources Investigations 6-A1, U.S. Geological Survey, Reston, Virginia.
- Harbaugh AW, Banta ER, Hill MC, McDonald MG (2000) MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model: User Guide to Modularization Concepts and the Ground-Water Flow Process. USGS Open-File Rep. 00-92, U.S. Geological Survey, Reston, Virginia.
- Harbaugh AW (2005) MODFLOW-2005, the U.S. Geological Survey Modular Ground-Water Model-The Ground-Water Flow Process. USGS Techniques and Methods 6-A16, U.S. Geological Survey, Reston, Virginia.
- Zheng C, Wang P (1999) MT3DMS-a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reaction of contaminants in ground-water systems: Documentation and user’s guide. Jacksonville, Florida. Contact Report SERDP-99-1, U.S. Army Corps of Engineers.
- Clement TP (1997) RT3D-A Modular Computer Code for Simulating Reactive Multi-Species Transport in 3-Dimensional Groundwater Aquifers. Pacific Northwest National Laboratory, Richland, Washington.
- Hill MC, Poeter E, Zheng C, Doherty J (2003) MODFLOW 2001 and other modeling Odysseys. Ground Water 41: 113-113.
- Henry HR (1964) Effects of dispersion on salt encroachment in coastal aquifers. Sea water in coastal aquifers. Geological Survey Water Supply Paper 1613-C, U.S. Geological Survey, Washington, DC, pp: 70-84.
- Simpson MJ, Clement TP (2003) Theoretical analysis of the worthiness of Henry and Elder problems as benchmarks of density-dependent groundwater flow models. Adv Water Resour 26: 7-31.
- Simpson MJ, Clement TP (2004) Improving the worthiness of the Henry problem as a benchmark for density-dependent groundwater flow models. Water Resour Res 40: 1-11.
- Dentz M, Tartakovsky DM, Abarca E, Guadagnini A, Sánchez-Vila X, et al. (2006) Variable- density flow in porous media. J Fluid Mech 561: 209-235.
- Goswami RR, Clement TP (2007) Laboratory-scale investigation of saltwater intrusion dynamics. Water Resour Res 43: W04418.
- Abarca E, Carrera J, Sánchez-Vila X, Dentz M (2007) Anisotropic dispersive Henry problem. Adv Water Resour 30: 913-926.
- Arlai P, Koch M (2007) Need for density-dependent flow and transport modeling of horizontal seawater and vertical saltwater intrusion in the Bangkok multilayer-aquifer system. Proceedings of the 12th National Convention on Civil Engineering, Phitsanulok, Thailand, May 2-4.
- Arlai P, Koch M (2009) The importance of density-dependent flow and solute transport modeling to simulate seawater intrusion into a coastal aquifer system. Proceedings of the International Symposium on Efficient Groundwater Resources Management (IGS-TH 2009), Bangkok, Thailand, February 16-21.
- Motz L, Sedighi A (2013) Saltwater intrusion and recirculation of seawater at a coastal boundary. J Hydrol Eng 18: 10-18.
- Brown DW, Kenner WE, Crooks JW, Foster JB (1962) Water resources of Brevard County, Florida. Report of Investigations 28, U.S. Geological Survey, Tallahassee, FL.
- Martin JB, Cable JE, Smith C, Roy M, Cherrier J (2007) Magnitudes of submarine groundwater discharge from marine and terrestrial sources: Indian River Lagoon, Florida. Water Resour Res 43: W05440.
Citation: Al-Taliby W, Pandit A, Heck H (2017) Comparison of Seepage Simulation in a Saline Environment below an Estuary Using MODFLOW and SEAWAT. Hydrol Current Res 8: 270. Doi: 10.4172/2157-7587.1000270
Copyright: © 2017 Al-Taliby W, et at. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Select your language of interest to view the total content in your interested language
Share This Article
- Total views: 1510
- [From(publication date): 0-2017 - Aug 21, 2019]
- Breakdown by view type
- HTML page views: 1403
- PDF downloads: 107