Groundwater Flow Modeling in Coastal Aquifers: The Influence of Submarine Groundwater Discharge on the Position of the Saltwater– Freshwater Interface

An investigation of the impact of submarine groundwater discharge on the position of saltwater-freshwater interface is presented in this manuscript. Two conceptualizations were considered and analyzed using both analytic and numerical techniques, for comparison purposes. The first conceptualization assumes that the tip of the saltwater- freshwater interface occurs at the shoreline, and the second conceptualization allows for the tip to extend off-shore. Analytic solutions exist for both conceptualizations. Results from both analytic and numeric analysis for the two conceptualizations are presented. Results from the first conceptualization were found to overestimate the inland distance to the interface toe, compared to the second conceptualization, for it ignores the influence of submarine groundwater discharge on the interface location. Moreover, results from the analytic solutions as a whole were found to overestimate the interface location compared to the numerical modeling results, for analytic solutions are based on the sharp interface approximations. Therefore, an empirically derived dispersion factor should be used to correct the analytic solution results so as to compare them with the numerically simulated values. Furthermore, offshore model extents should be incorporated when modeling coastal aquifer systems to include the influence of submarine groundwater discharge on the saltwater-freshwater interface position.


Introduction
Seawater intrusion is the migration of saline water into freshwater in coastal aquifers. Saline water is denser than freshwater, for it has higher mineral contents. Consequently, it forms a wedge beneath freshwater in coastal aquifers ( Figure 1). Seawater intrusion can occur naturally owing to the connectivity between seawater and groundwater, and due to certain human activities. Therefore, modeling the coastal groundwater flow system enables the evaluation of the potential for seawater intrusion into aquifer systems as a result of different factors. However, modeling saltwater intrusion is considered difficult. Factors such as heterogeneity of the aquifer hydraulic properties, the complicated aquifer geometries and the temporal and spatial variability in groundwater density make it difficult to model seawater intrusion [1]. The accuracies of model outputs are strongly based on the assumptions made on the model input parameters. Seawater intrusion model results for the saltwater-freshwater interface position, for example, are strongly affected by different factors such as boundary conditions, initial (head and concentration) conditions and aquifer hydraulic properties. Furthermore, nowadays, Submarine Groundwater Discharge (SGD) is also becoming an important issue to be considered in modeling coastal groundwater systems. Owing to seawater intrusion, the land driven fresh groundwater can discharge to the seafloor through the leaky confining unit and the process is called SGD [2][3][4]. This kind of discharge decreases with the increase in distance offshore and is zero where the tip of the interface touches the leaky confining unit [5].
Several authors have studied seawater intrusion and the position of the saltwater-freshwater interface owing to different factors in coastal aquifers. Strack [6] developed an analytic solution for the regional interface problems in coastal aquifers based on the single-valued potentials, the Dupuit assumption and the Ghyben-Herzberg formula for the steady state flow conditions. The Strack [6] analytic solution has *Corresponding author: Haile Arefayne Shishaye, Haramaya University, Dire Dawa, Ethiopia, Tel: +251942124463; E-mail: haile.4.hiwot@gmail.com been widely used by different researchers to explore seawater intrusion in coastal aquifers [1,5,[7][8][9][10]. Different seawater intrusion assessment methods have also been developed based on the Strack [6] analytic solution [11,12].
Other authors like Huyakorn et al. [13], who presented a numerical model based on the sharp interface approach and taking into account the flow dynamics of saltwater and freshwater, Motz [14], who proposed an analytic solution for calculating the critical pumping flow rate in an artesian aquifer, and Bower et al. [15], who modified the critical interface rise based on the analytic solution which allows the critical pumping rate to be increased are also some of the well-known studies conducted on seawater intrusion in costal aquifers.
However, none of the above papers consider the influence of SGD on the seawater-freshwater interface position. There is no possibility to simulate the offshore groundwater discharges using the above analytic solutions, unless modifications are made to include the offshore outflow zone of the land driven fresh groundwater through the seafloor by taking model extent offshore into consideration. Recently, Bakker [16] has modified the Strack solution so as to include the offshore freshwater outflow zone. It is a solution for a steady interface flow in confined coastal aquifers discharging to a semi-confined section below the ocean. Bakker has shown that the tip of the saltwater-freshwater interface can perhaps touch the leaky confining unit at some distance offshore, and this depends on the head of the land driven fresh groundwater and the leakage factor of the seafloor. Hence, the decision on how long a model should be extended offshore for accurate simulations of the interface position is also an important consideration when modeling coastal groundwater flow systems [17]. In general, to identify the knowledge gap in our current advances in seawater intrusion and the influence of submarine groundwater discharge on the position of the saltwaterfreshwater interface, a comparison among the previously developed and widely in use analytical and numerical models on the area is very crucial. Such a comparison can show the gaps within the developed models and can also give directions to fill the knowledge gaps in the area.
The objective of this manuscript is therefore to investigate the influence of SGD on the position of saltwater-freshwater interface. To do so, comparing the steady state interface location when two conceptualizations are used in both analytic and numerical modeling techniques could be very important. The first conceptualization is based on the Strack [6] analytic solution, assuming that the tip of the interface lies at the shoreline; while, the second conceptualization is based on the Bakker [16], taking the distance offshore into consideration. [ * 2

Common parameters and values used
The discharge potential Φ is given by i X Q 0 , where Q 0 is the inflow to the confined aquifer and X i is the onshore distance. If the toe of the interface is assumed to be located at where, Ф d is the discharge potential at a distance ofd from the coast. Bakker [16] analytic solution: The vertically integrated freshwater discharge of the confined aquifer in the horizontal-direction is given in the Bakker (2006) analytic solution by: where, f h is the thickness of the freshwater zone. From this, the following procedures were followed to derive simplified equations, based on the Bakker [16], for the discharge potential at the toe (where the thickness of the freshwater zone is equal to the thickness of the aquifer, i.e., H h f = ), lengths onshore and offshore and the interface heads.
As described above, the thickness of the freshwater zone is equal to the thickness of the aquifer i.e., f h H = at the toe of the interface. Substituting f h by , H both sides of equation (3) were integrated with respect to x and , h respectively and yielded the following equation.
The discharge potential at any distance i X = is given by the product of the vertically integrated discharge in the confined aquifer and the distance " "i from the shoreline. The value of X is 0 at the shoreline. Therefore, if it is assumed that the toe of the interface is at a distance d X = from the shoreline, then the discharge potential at the toe will be given by In this case, if it is assumed that the tip of the interface lays at the coastline, the distance of the toe from the shoreline will be given by: However, Bakker [16] has taken the distance offshore into consideration. Thus, the place where the discharge potential values become zero will not necessarily be at the coordinate where 0 = X . Therefore, the point where the discharge potential becomes zero lies where the tip of the interface touches the confining unit. So, the procedure to develop an equation for the discharge potential at the shoreline based on the Bakker (2006) analytic solution is similar to that of followed above, but with different head value. The interface head at the shoreline lies at a depth 0 Z below the sea level or (  (4) above will give the following equation for the discharge potential at the shoreline.
The discharge potential ( 0 Φ ) is zero and 0 Z Hl = at 0 = X when the tip of the interface lays at the shoreline. From here, the equation for the distance of the interface toe from the shoreline might be different from equation (5), if the tip of the interface lays at some distance offshore. i.e, The model extent offshore (L) is given in the Bakker (2006) analytic solution as: Therefore, equation (8) has been re-written as follows to develop an expression for the distance offshore (L) in terms of the parameters given in (Table 1). So, Bakker has also developed an equation relating the distance offshore (L) and the dimensionless head (φ ), i.e., Substituting these values will then give us: Therefore, the interface head at the toe is calculated as follows: . This can also be written as: Thus, the depth of the interface below sea level ( 0 Z ) is calculated as follows: Hence, the plotting plane can be divided into two zones, the offshore and onshore zones. The saltwater-freshwater interface heads along the two zones can be plotted against the offshore and onshore distances. For example, let " X1 " represents the list of equally spaced 100 numbers from 0 to -d and " X2 " represents the list of equally spaced 100 numbers from 0 to L. It is also possible to increase or decrease the list of numbers within the range to increase the plotting accuracy. Therefore, in this case, plotting can be completed within two steps.
Step 1 is to plot X1 against (Z1 +30) and step 2 is to plot X2 against (Z2 + 30), where, " 30 " is the depth of the bottom of the aquifer below sea level and Z1 and Z2 are calculated as follows.
From equation (4), Therefore, to derive an equation for the depth of the interface below the confining unit at any distance 1 X ( 1 fx h ), we need to calculate the total discharge potentials at any distance . 1 X The value of Ф at a distanceX1 is given by 1 0 X Q . Moreover, the total discharge potential at a distance of X1 is 0 0 ( 1 ) Q X + Φ . But because of that X1 represents a list of 100 numbers from 0 to -d, a negative sign should be included within the equation Q 0 X1.Therefore, substituting (12), the depth of the interface ( 1 fx h )below the confining unit at any distance of X1 is given by: The depth of the interface below sea level is, therefore, calculated as Further, the depth of the interface ( 2 fx h ) at any distance offshore (X2) can also be calculated from eqn. (10). The point where the interface tip touches the leaky confining unit is when Hl is equal with Z2. This point is located at a distance where Therefore, because of that X2 has 100 list of values, substituting L 2 by (X2-L) 2 would help to calculate the depth of the interface below the confining unit at the 100 different X2 values. Hence, the depth of the interface below the leaky confining unit ( 2 fx h ) will be given as: The depth of the interface below sea level can then be calculated as: The depth from the sea level to the bottom of the aquifer is 30 (i.e., H b = -30). Furthermore, the depth to the interface locations offshore is represented by Z2 and to that of onshore is represented by Z1. Therefore, Z1 + 30 will result the hydraulic heads of the interface all along 0 to d, and Z2 + 30 will give the interface hydraulic heads all along 0 through -L. Remember that, the values for Z1 and Z2 are based on X1 and X2 (from equations 13 and 15), respectively. Therefore, Z1 and Z2 will represent a list of 100 numbers each. This implies that we have 100 points onshore and 100 points offshore to plot.
Numerical modeling: In addition to their use as planning tools for improving water supply and management, numerical models of groundwater systems are also useful for understanding groundwater flow processes. In terms of the use of modeling packages, groundwater flow systems can be divided into two, i.e., groundwater flow processes with constant density and the one with variable density. Simulating the groundwater systems in coastal aquifers which include saltwater and freshwater requires the use of a numerical modeling code that solves the variable-density flow equation (examples and perhaps widely used packages are SEAWAT and SUTRA).
SEAWAT is a generic MODFLOW/MT3DMS based computer program designed to simulate three-dimensional variable-density groundwater flow coupled with multi-species solute and heat transport. While, SUTRA is a general-purpose, density-dependent, fluid flow and mass-transport numerical model that applies a finite element and integrated finite-difference hybrid method, which is mainly used to model both the coastal surficial and confined aquifers [17]. In this case, the model used to investigate the impact of SGD on the position of saltwater-freshwater interface is the three dimensional SEAWAT model. SEAWAT has been used widely for groundwater studies including saltwater intrusion in coastal aquifers.
The type of aquifer considered in this paper is a confined coastal aquifer which is hydrogeologically connected with the seawater. Similar to what was done in the analytic modeling section, numerical modeling was conducted based on the Strack [6] and Bakker [16] analytic solutions for the interface problems. The other consideration, in this simulation was that the confining unit offshore (the sea floor) is assumed to be a leaky confining unit. The common parameters used in both simulation cases are listed in Table 1.
Case 1-Modeling with no distance offshore: To obtain a steady state solution, the simulation run was divided into 10 stress periods, which in turn are divided into 10,000 time steps and 70,000 days of period length each, which corresponds to a total simulation period of 700,000 days. Modeling was conducted for case-1 by constructing a three dimensional SEAWAT model with an inland distance of 2600 meters, based on the Strack [6] analytic solution. The model was initially run for 50,000days in a steady state simulation flow type. Then, the initial and prescribed head was taken from the steady state simulation as an input for the transient simulation flow type. Therefore, the initial and prescribed head used in this simulation was 30m. The initial concentrations were based on seawater concentrations, 35 kg/m 3 for columns 2 to 129 in all layers; while, column 1 (seawater column) was fixed to a concentration of 35 kg/m 3 . The concentration for column 130 (freshwater column) was also set to 0. A uniform and isotropic value for hydraulic conductivity was set to 260 m/d, the porosity was set to 0.35, and the values for longitudinal and transverse dispersivities were set to 0.1m each. The specific storage was also set to 0.0001.
The SIP package of MODFLOW and the GCG package with the finite-difference option of MT3DMS/SEAWAT were used to solve systems of the flow and transport equations, respectively. The SIP solver was used with a head convergence criterion of 1 x 10 -4 m to solve the flow equation. Furthermore, the GCG solver was used with a courant number (number of cells or fraction of a cell that a parcel of water can advect during one time-step) of 0.75 to solve the solute-transport equation.
Case 2-Modeling with distance offshore: Modeling was conducted for case-2 by constructing a three dimensional SEAWAT model with an inland and offshore distances of 2600m and 400m, respectively, based on the Bakker [16] analytic solution.
The same dimensions, initial conditions, solver packages and aquifer properties to that of used in case-1 were used in case-2. However, the boundary conditions were different between the two cases. The left, bottom, and from column 21 to 149 of the top layer were set to no flow boundaries. The right hand side (column 150) was set to constant flux with a constant inflow rate of 1 m 3 and density of 0 kg/m 3 . From column 1 to 20 of the top boundary, a general head boundary was applied with an external head of 10 m, conductance of 2 m 2 /d and density of 35 kg/ m 3 . This indicates that the model takes the SGD into consideration over the offshore distance specified above.
In this case, the impact of the seawater is through the general head boundary (the seafloor leakage). The vertical seawater column simulated in case-1 is now neglected. However, when considering a long model extent offshore, up to the end of the continental shelf, the vertical constant head and density seawater column should be taken into account. The only reason for neglecting the vertical seawater column and simulating only the impact of the seawater through the leaky confining unit (General head boundary) in our case is because of that the offshore model extent considered is very short (400 m). Otherwise, both the vertical constant head and density seawater column and the leakage through the seafloor (in a semi confined scenario) should be taken into account during actual simulations.

Result and Discussion
Analytic solution Case 1: The Strack [6] analytic solution is based on the use of the single potential which is defined throughout the multiple zones of an aquifer. Refer to Strack [6] to look at the assumptions based on which he developed his analytic solution.
However, this conceptualization only works in coastal aquifers with no SGD. Furthermore, the solution neglects the mixing factor, for it is based on the sharp interface approximation. Nevertheless, Strack's solution has been widely used in the area of coastal groundwater flow modeling. Therefore, taking the above assumptions into consideration and using the common parameter values listed in table 1, and equations 1 and 2, a python script of this conceptualization was developed, which has yielded the graph shown in Figure 2. In this case, the graph is showing that there is no offshore out flow zone included in the Strack's solution. The tip of the interface is located exactly at the shoreline; while, the toe was found at a distance of 1300m inland.

Case 2:
As indicated in the methodology part above, a python script was developed based on the Bakker [16] analytic solution. Common parameter values were also used with the Strack [6] analytic solution for comparison purposes. Similar to case-1, the interface heads vs. distance graph was plotted using the python package ( Figure 3).
Bakker's analytic solution was developed for steady interface flow in aquifers consisting of confined and semi-confined sections. According to Bakker [16], the integrated discharge from all layers is constant in the confined section and is directed towards the semi-confined section, which is bounded on top by a leaky seafloor. The land driven fresh groundwater flows up through the leaky confining unit (seafloor). This discharge is, therefore, totally based on the freshwater-saltwater head differences, hydraulic properties of the aquifer materials and the leakage factor of the seafloor. Using the common values (Table 1), the toe of the interface in this case was found at an inland distance of 1222.97 m; while, the tip was found at an offshore distance of 308.11 m.
Plotting the two cases on one graph: The main difference between the two well known analytic solutions for the interface problems, Strack [6] and Bakker [16], is the distance offshore. Bakker [16] has taken the distance offshore into consideration to include the influence of SGD on the seawater-freshwater interface position; while, Strack [6] has not. This difference has created a difference on the position and shape of the interfaces developed using the two solutions.
As shown in Figure 4, the toe of the interface was found at a distance of 1222.97 m from the coastline using the Bakker [16] analytic solution, while it was found at a distance of 1300m using the Strack [6] analytic solution. In this case, the SGD has created (1300 m-1222.97 m), 77.03 m, gap in the location of the toe of the interface. In other words, the Strack [6] analytic solution has overestimated the location of the toe of the interface, for it doesn't consider the model extent offshore.
It is highly unlikely to find the saltwater column vertically fixed at the shoreline. Even though it varies from place to place, the continental shelves can perhaps be extended to hundreds and/or thousands of kilometers offshore. There are plenty of scientific evidences for the availability of SGD [2][3][4][19][20][21][22][23]. Therefore, it is important to consider the model extent offshore for accurate simulations of the position of saltwater-freshwater interface.
As shown in Figure 4, the length of the outflow zone was found to be 308.11 m, using the Bakker [16] analytic solution; while, there is no outflow zone considered in the Strack [6] analytic solution. In fact, the gap between the simulated interface locations using the Bakker and Strack analytic solutions is getting wider from the toe to the tip. And this is because of the leakage factor included in the Bakker [16] analytic solution. The Bakker [16] analytic solution assumes two sections within the aquifer system, the confined and semi-confined conditions. Therefore, the appearance of the interface can possibly vary based on the values used for the leakage factor.
Thus, according to the simulation results, errors can be caused by ignoring the SGD during simulations. The extent of errors in confined aquifers may be greater than that of unconfined aquifers, for the flow

Numerical modeling
Case 1: This conceptualization was based on the Strack [6] analytic solution for the interface problems. In this case, the concentration contour line at the 50% of the maximum concentration was chosen to represent the interface location. Accordingly, the tip of the seawaterfreshwater interface was found at the shoreline; while, the toe was found at an inland distance of 880m ( Figure 5). The distance to the toe location found in this case is shorter than the distance found in case-1 of the analytic solution section. Both the analytic and numerical solutions of this case (case-1) were applied for the same conceptualization. However, the analytic solution has overestimated the interface location, for it ignores the mixing factor. As it is shown in Figure 6, the saline water (dark brown color) is overlaid by the mixing zone which in turn is overlaid by the land driven fresh groundwater discharge. However, this structure is missed on the analytic solution.

Case 2:
This conceptualization was based on the Bakker [16] analytic solution for the interface problems. Similar to case-1 above, the concentration contour line at the 50% of the maximum concentration was chosen to represent the interface location. Accordingly, the tip of the seawater-freshwater interface was found at an offshore distance of 300m; while, the toe was found at an inland distance of 700m ( Figures  7 and 8). Similar to case-1, the distance to the toe location found in this case is shorter than the distance found in case-2 of the analytic solution section.
Comparison between case-1 and case-2 results of the numerical modeling: Similar to the analytic solution, the numerically simulated results for case-1 and case-2 are also different. Furthermore, the locations and shapes of the interfaces in these two cases are different. The main reason for this situation is, therefore, the SGD incorporated in case-2. Therefore, ignoring the influence of SGD on the interface position when modeling coastal groundwater systems, especially those with confined and semi-confined sections, overestimates the interface location.

Comparison between analytic and numerical modeling results:
Analytic solutions for the position of the saltwater-freshwater interface are based on the sharp interface approximations. The advantages of analytic solutions are that they are considerably less computationally intensive and require less data [11]. However, sharp interface approximations can only be used in areas where the mixing zone can be ignored; while, in reality, mixing between freshwater and saltwater occurs in all coastal aquifers as a result of dispersion, changes in head, changes in hydraulic conductivity and tides [25]. In fact, the extent of mixing varies from a few tens of centimeters in tight clays or sandstone to several tens of meters in karstic limestone [26]. But mixing is always there, regardless of its extents. Therefore, analytic solutions obviously overestimate the extent of saltwater penetration further inland.
In reality, freshwater overlies the mixing zone which in turn overlies saline water. It is not, therefore, possible to provide stable and accurate results using analytical solutions, unless a correction factor is incorporated to include the influence of mixing. However, the complex density-dependent groundwater flow and solute transport models provide stable and convincing results when employed with proper spatial and temporal discretisations.
In this case, simulation results for the interface heads from the two analytic solutions were corrected by the empirically derived dispersion factor developed by Pool and Carrera [12] to include the influence of mixing on the interface location, where α T is transverse dispersivity and b is aquifer thickness. While, simulation results for    systems in equilibrium in confined aquifers may discharge hundreds of meters or even kilometers offshore depending on the hydraulic gradient and the geology of the formation [24]. However, regardless of whether conditions are confined or unconfined, error occurs in Strack's solution because it ignores the SGD.   the saltwater-freshwater interface position from the density-dependent Solute Transport Numerical Model (SEAWAT), which includes advection and dispersion, are believed to be accurate.
The transverse dispersivity used in both of the two cases of the numerical modeling section was 0.1. Therefore, the correction factor (f) can be calculated as follows:  Figures 9 and 10, the analytic solution results for both case-1 and case-2 were corrected by the correction factor (f). Finally, the corrected analytic solutions became very close to the numerically simulated values. Initially, the location of the toe for case-1 was found at an inland distance of 1300 m; while, it was found at a distance of 880 m from the numerical modeling. However, after correcting the analytic solution result, the toe was found at an inland distance of 762.43 m (Figure 9), which is very close to the numerically simulated value than the initial one.
Similarly, the analytically calculated toe location for case-2 was also corrected by the empirically derived dispersion factor (0.5865). Initially, the location of the toe for case-2, was found at an inland distance of 1222.97 m; while, it was found at a distance of 700 m from the numerical modeling. However, after correcting the analytic solution result, the toe was found at an inland distance of 717.25 m (Figure 10), which is also very close to the numerically simulated value than the initial one.

Conclusion
In conclusion, two inferences can be derived from this work. Firstly, SGD has an impact on the seawater-freshwater interface position. Hence, model extents offshore should be taken into account whenmodeling groundwater flow in coastal aquifer systems to incorporate the influence of SGD on the interface position. In this case, simulation results based on the Strack [6] analytic solution for the interface problems have overestimated the interface location, compared to the results based on Bakker [16] analytic solution. The reason for this situation is that the Strack [6] solution neglects model extents offshore. Thus, it needs to be modified to incorporate the influence of SGD on the interface position.
Secondly, analytic solutions overestimate the location of the seawater-freshwater interface, for they are based on the sharp interface approximations. Therefore, results from analytic solutions need to be corrected by the empirically derived dispersion factor (f) to include the mixing factor in the solution so as to compare results from numerical modelling with those from the analytical solutions.