Technical Note : Groundwater flow modeling in coastal aquifers – the influence of submarine groundwater discharge on the position of the saltwater – freshwater interface

Introduction Conclusions References Tables Figures

Technical Note: Groundwater flow modeling in coastal aquifers -the influence of submarine groundwater discharge on the position of the saltwater-freshwater interface coming 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 (Post et al., 2013;Moore, 2009;Church, 1996). This kind of discharge decreases with the increase in distance offshore and is zero where the tip of the interface touches the leaky 15 confining unit (Beebe et al., 2011).
Several authors have studied seawater intrusion and the position of the saltwaterfreshwater interface owing to different factors in coastal aquifers. Strack (1976) 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 for-20 mula for the steady state flow conditions. The Strack (1976) analytic solution has been widely used by different researchers to explore seawater intrusion in coastal aquifers (e.g., Morgan et al., 2013;Beebe et al., 2011;Aharmouch and Larabi, 2001;Wriedt and Bouraoui, 2009;Mazi, 2014;Naderi et al., 2013). Different seawater intrusion assessment methods have also been developed based on the Strack (1976)  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | critical pumping flow rate in an artesian aquifer, and Bower et al. (1999), 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-5 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 (2006) 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 15 a model should be extended offshore for accurate simulations of the interface position is also an important consideration when modeling coastal groundwater flow systems. 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 20 modeling techniques could be very important. The first conceptualization is based on the Strack (1976) analytic solution, assuming that the tip of the interface lies at the shoreline; while, the second conceptualization is based on the Bakker (2006), taking the distance offshore into consideration. Madras aquifer (in the city of Madras, now called Channai, India) by Sherif and Singh in 1999. These values were then used in both conceptualisations of the analytic and numerical modellings.

Common parameters and values used
in place of h f will give the following: The discharge potential Φ is given by Q 0 X i , 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 HESSD 12,2015 where, Φ d is the discharge potential at a distance of d from the coast.

Bakker (2006) analytic solution
The vertically integrated freshwater discharge of the confined aquifer in the horizontal-5 direction is given in the Bakker (2006) analytic solution by: where, h f is the thickness of the freshwater zone. From this, the following procedures were followed to derive simplified equations, based on the Bakker (2006), for the discharge potential at the toe (where the thickness of the freshwater zone is equal to the 10 thickness of the aquifer, i.e., h f = H), 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., h f = H) at the toe of the interface. Substituting h f by H, both sides of Eq. (3) were integrated with respect to x and h, respectively and yielded the following 15 equation.
The discharge potential at any distance X = i 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 X = d from the shoreline, then the discharge potential at the toe will be given by Φ d = Q 0 d . Thus, because of Q 0 d = 1 2 δK H 2 , 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 (2006) 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 X = 0. Therefore, the point where the discharge potential becomes zero lies where the tip of the interface touches the confining unit. So, the procedure to 10 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 Z 0 below the sea level or (Hl − Z 0 ) below the bottom of the confining unit, where Hl is the depth of the confining unit below the sea level. Therefore, substituting (Hl − Z 0 ) 2 in place of H 2 in Eq. (4) above will give 15 the following equation for the discharge potential at the shoreline.
The discharge potential (Φ 0 ) is zero and Z 0 = Hl at X = 0 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 Eq. (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: The equations for g c and λ d are also given as Q 0 K H and λ H , respectively. Therefore, Eq. (8) has been re-written as follows to develop an expression for the 5 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., According to Bakker (2006), ϕ = 1 and u = 0 when the toe lies at the intersection point of the confined and semi-confined sections. Substituting these values will then give us: 15 Therefore, the interface head at the toe is calculated as follows: (thickness of the freshwater head). This can also be written as: Thus, the depth of the interface below sea level (Z 0 ) is calculated as follows: Hence, the plotting plane can be divided into two zones, the offshore and onshore 5 zones. The saltwater-freshwater interface heads along the two zones can be plotted against the offshore and onshore distances. For example, let "X 1" represents the list of equally spaced 100 numbers from 0 to −d and "X 2" 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, 10 plotting can be completed within two steps.
Step 1 is to plot X 1 against (Z1 + 30) and step 2 is to plot X 2 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 Eq. (4) Therefore, to derive an equation for the depth of the interface below the confining unit at any distance X 1 (h fX 1 ), we need to calculate the total discharge potentials at any 3761 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | distance X 1. The value of Φ at a distance X 1 is given by Q 0 X 1. Moreover, the total discharge potential at a distance of X 1 is (Q 0 X 1 + Φ 0 ). But because of that X 1 represents a list of 100 numbers from 0 to −d , a negative sign should be included within the equation Q 0 X 1. Therefore, substituting (−Q 0 X 1 + Φ 0 ) in place of Φ d in Eq. (12), the depth of the interface (h fX 1 ) below the confining unit at any distance of X 1 is given by: The depth of the interface below sea level is, therefore, calculated as Z1 = Hl − h fX 1 .
Further, the depth of the interface (h fX 2 ) at any distance offshore (X 2) can also be calculated from Eq. (10). The point where the interface tip touches the leaky confining 10 unit is when Hl is equal with Z2. This point is located at a distance where X 2 = L. Therefore, because of that X 2 has 100 list of values, substituting L 2 by (X 2−L) 2 would help to calculate the depth of the interface below the confining unit at the 100 different X 2 values. Hence, the depth of the interface below the leaky confining unit (h fX 2 ) will be given as: The depth of the interface below sea level can then be calculated as: 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 X 1 and X 2 (from Eqs. 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 den- to simulate three-dimensional variable-density groundwater flow coupled with multispecies solute and heat transport. While, SUTRA is a general-purpose, densitydependent, 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 (Werner et al., 2013). In this case, the model 20 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 25 modeling section, numerical modeling was conducted based on the Strack (1976) and Bakker (2006) analytic solutions for the interface problems. The other consideration, in this simulation was that the confining unit offshore (the sea floor) is assumed to 3763 Introduction 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, 5 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 m, based on the Strack (1976)  while, the right side boundary was set to a constant flux freshwater with an inflow rate of 1 m 3 and density of 0 kg m −3 .

15
The model was initially run for 50 000 days 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 30 m. The initial concentrations were based on seawater concentrations, 35 kg m −3 for columns 2 to 129 in all layers; while, column 1 (seawa-20 ter 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 day −1 , the porosity was set to 0.35, and the values for longitudinal and transverse dispersivities were set to 0.1 m each. The specific storage was also set to 0.0001.

25
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 3764 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of 1 × 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 5 model with an inland and offshore distances of 2600 and 400 m, respectively, based on the Bakker (2006) 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 10 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 day −1 and density of 35 kg m −3 . This indicates that the model takes the SGD into consideration over the offshore distance specified above. 15 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 simu-20 lating 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. 12,2015  The Strack (1976) analytic solution is based on the use of the single potential which is defined throughout the multiple zones of an aquifer. Refer to Strack (1976) to look at 5 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.

HESSD
Therefore, taking the above assumptions into consideration and using the common parameter values listed in Table 1, and Eqs. (1) and (2), a python script of this conceptualization was developed, which has yielded the graph shown in Fig. 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 15 found at a distance of 1300 m inland.

Case-2
As indicated in the methodology part above, a python script was developed based on the Bakker (2006) analytic solution. Common parameter values were also used with the Strack (1976) analytic solution for comparison purposes. Similar to case-1, the 20 interface heads vs. distance graph was plotted using the python package (Fig. 3).
Bakker's analytic solution was developed for steady interface flow in aquifers consisting of confined and semi-confined sections. According to Bakker (2006), 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 HESSD 12,2015  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 (1976) and Bakker (2006), is the distance offshore. Bakker (2006) has taken the distance offshore into consideration to include the influence of SGD on the seawater-freshwater interface position; while, Strack (1976) has not. This difference 10 has created a difference on the position and shape of the interfaces developed using the two solutions. As shown in Fig. 4, the toe of the interface was found at a distance of 1222.97 m from the coastline using the Bakker (2006) analytic solution, while it was found at a distance of 1300 m using the Strack (1976) analytic solution. In this case, the SGD has created 15 (1300-1222.97 m), 77.03 m, gap in the location of the toe of the interface. In other words, the Strack (1976) analytic solution has overestimated the location of the toe of the interface, for it does not 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 20 to hundreds and/or thousands of kilometers offshore. There are plenty of scientific evidences for the availability of SGD (Post et al., 2013;Moore, 2009Moore, ,1996Johannes, 1980;Fanning et al., 1981;Church, 1996;Li et al., 1999;Burnett et al., 2006). Therefore, it is important to consider the model extent offshore for accurate simulations of the position of saltwater-freshwater interface.

25
As shown in Fig. 4, the length of the outflow zone was found to be 308.11 m, using the Bakker (2006) analytic solution; while, there is no outflow zone considered in the Strack (1976) analytic solution. In fact, the gap between the simulated interface locations using 3767 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 (2006) analytic solution. The Bakker (2006) 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. 5 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 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 (Kooi and Groen, 2001). However, regard-10 less of whether conditions are confined or unconfined, error occurs in Strack's solution because it ignores the SGD.

Case-1
This conceptualization was based on the Strack (1976) analytic solution for the inter- 15 face 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 seawater-freshwater interface was found at the shoreline; while, the toe was found at an inland distance of 880 m (Fig. 6). 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 20 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 Fig. 5, 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 (2006) 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 5 300 m; while, the toe was found at an inland distance of 700 m (Figs. 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-10 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 semiconfined 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 (Werner et al., 2012). However, sharp interface approximations can only be used in areas where the mixing 20 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 (Shalev et al., 2009). 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 (Dausman and Langevin, 2005 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, 5 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 [1 − (α T /b) 1/6 ] developed by Pool and Carrera (2011) to include the influence of mixing on the interface location, where α T is transverse dispersivity and b is aquifer thickness. While, simulation results for the saltwater-freshwater interface position from the density-dependent solute transport numerical model (SEAWAT), which includes advection and dispersion, are believed to be accurate. 15 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: As shown in Figs. 9 and 10, the analytic solution results for both case-1 and case-20 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 (Fig. 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 5 from the numerical modeling. However, after correcting the analytic solution result, the toe was found at an inland distance of 717.25 m (Fig. 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 when modeling 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 (1976) analytic solution for the interface problems have overestimated the interface location, compared to the results based on Bakker (2006) analytic solution. The reason for this situation is that the Strack (1976) 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 20 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.