Received 21 February 2016; accepted 26 April 2016; published 29 April 2016

1. Introduction

Coastal aquifers are a widely relied source of fresh water supply to the population residing in the coastal areas in many parts of the world. Due to the overuse of groundwater in such areas around the world, aquifers are subjected to threats of seawater encroachment. Intrusion of seawater into coastal aquifers poses a serious challenge in the utilization of groundwater [1] . It occurs as an outcome of more groundwater being pumped from coastal aquifers in excess of that which can be replenished causing the decline of piezometric head. This results in the encroachment of seawater into the fresh water aquifer [2] . For solving seawater intrusion problem, a coupled groundwater flow and transport equation is required to be solved. The problem is nonlinear and complex leading to high computational cost and hence most authors limit their domain to vertical 2D models. For the past two decades many researchers have modeled seawater intrusion using FEM [3] - [6] and FDM [7] [8] based models.

In this study, Meshfree (MFree) method is proposed for the analysis of seawater intrusion problem. Mesh Free method has been used earlier to model the groundwater flow and solute transport problems [9] - [17] . It was found to be more effective because of its advantages over the mesh based methods. Most widely applied FDM and FEM methods use a properly predefined grid or mesh. In the MFree method, the system of algebraic equations for the entire problem domain is solved using convenient nodes thereby avoiding a predefined mesh. The problem domain is represented in MFree method by a set of nodes which are scattered within the problem domain and also on the boundaries. There is no gridding required and hence when large scale problems are being solved, computational cost and time spent on pre-processing can be saved. Kansa [18] was the first to attempt radial basis functions for solving partial differential equations. In the present study a Point Collocation Method (PCM) based MFree model is developed for a two dimensional confined aquifer seawater intrusion problem and applied for two distinct synthetic seawater intrusion problems. The presently results are compared with solutions from well-established research problems.

2. Governing Equations and Boundary Conditions

In the diffusive interface model, there is a transition zone caused due to hydrodynamic dispersion across which density of the water gradually varies from that of the fresh water to seawater unlike sharp interface models where the transition zone is approximated by a sharp interface [19] . The diffusive interface problem of saltwater intrusion into coastal aquifers can be solved by simultaneous solution of coupled groundwater flow and advective dispersive equations. The first equation describes the flow of a variable density fluid, and the second equation considers the transport of dissolved salt [20] .

The flow equation in heterogeneous anisotropic aquifer in steady state is given by Bear [21] as:

(1)

where K_xx and K_zz are components of hydraulic conductivity tensor along the principal axes direction, h is the hydraulic head and is normalized density difference between contaminant and fresh water i.e.. The solute transport equation for 2D flow is given by Bear [21] as:

(2)

where D_xx, D_zz, D_xz, and D_zx are dispersion tensor; c is the solute concentration; V_x and V_z are the seepage velocity in x and z directions.

Both these equations are coupled using the fundamental equation that relates fluid density, ρ to salt water concentration, c as follows:

(3)

where ρ and ρ_f are contaminant fluid and freshwater density respectively. For the present study, some non- dimensional variables have been defined below as found in Rastogi et al. [8] as.

(4)

(5)

(6)

where d is the depth of the confined coastal aquifer, V is the average linear velocity, c^* is seawater concentration, α_L and α_T are the longitudinal and transverse dispersivities respectively. After replacing the above non dimensional variables in Equation (1) and (2) and after dropping primes in all terms for convenience, the following equations are obtained.

(7)

where K_zz may vary from 0.0 to 1.0, and

(8)

α_T varies from 0.0 to 1.0 and the terms D_xx, D_zz, D_xz are defined in 2-D x-z coordinate system as suggested by Bear [21] ,

(9)

(10)

(11)

Darcy’s equation for principal x and z axes can be written as

(12)

(13)

Equation (7), (8), (12), (13) constitute the governing equations for seawater intrusion in coastal aquifers. In Equation (8), V_x and V_z are defined as follows.

(14)

(15)

Here; where (16)

θ is the porosity of the medium and 1/A is the head gradient of the fresh water. Dropping primes in all the above equation leads to

(17)

Similarly

(18)

(19)

Dropping primes the above equation becomes

(20)

Equation (17) and (20) give the velocity in non-dimensional form along the principal x and z directions respectively.

3. MFree PCM Formulation for Steady State Diffusive Interface Model

The MFree model is found to be an alternative numerical approach which eliminates the well-known drawbacks in finite difference, finite element and boundary element methods which accompany the grid based methods. In the present study, PCM based MFree model is developed and subsequently verified for the seawater intrusion problem in confined coastal aquifers based on a diffusive interface approach.

For discretizing the flow and transport equations, a PCM based meshfree strong form formulation proposed by Liu and Gu [22] is adopted. Initially a trial solution needs to be defined for and as [22]

(21)

(22)

where, h_i is the unknown head at the i^th node, and c_i is the unknown concentration at the i^th node. f_i is the shape function which is again a function of i.e. radial basis function (RBF) as given in Liu and Gu [22] and n is the total number of nodes in the support domain. Multi quadratics radial basis function (MQRBF) used presently is defined as:

(23)

where,; C_s and q are shape parameters of the radial basis function. where α_c is the dimensionless size of support domain, and d_c is the characteristic length that can be related to the nodal spacing in the local support domain. The value of shape parameter α_c has a very profound impact on the accuracy of results. In the absence of an explicit method to find out the shape parameters, trial and error approach is adopted for its determination, which poses a certain limitation to the MQRBF method. The shape parameter q can be 0.5, 0.98 or 1.03 as suggested in Liu and Gu [22] . It is observed that standard Multi quadratic RBF has the q value of ± 0.5. In the present study the value of q is taken as 0.98 from Thomas et al. [15] . The derivatives of the unknown function at the collocation node x_i can be expressed as

(24)

For present study, the shape of the support domain is taken as circular where the radius can be varied according to accuracy required. Kovarik and Muzik [16] also used the concept of prescribing a radius for determining the supporting points. With a larger radius, more number of nodes will be present in the support domain but it will lead to an increase in the time of computation. Support domain can also be rectangular or square as adopted by Meenal and Eldho [13] . By implementing both Equation (21) and Equation (22) and also by using PCM, Equation (7) can be written as:

(25)

where, , and values are to be calculated for each support domain

which are subsequently incorporated in the global matrix for the whole problem domain with consideration of Kronecker delta property.

Further Equation (25) can be written in matrix form as:

(26)

Similarly the transport Equation (8) can be approximated using PCM as

(27)

Equation (27) can be simplified and written as,

(28)

where [G₂] and [G₃] are global matrices of single derivative of f_i with respect to x and z respectively [G₄], [G₅], [G₆] and [G₇] are global matrices of double derivative of f_i with respect to x, z, xz, and zx respectively.

4. Model Development

Based on PCM based MFree formulation, a model in 2D was developed for diffusive interface seawater intrusion problem for a confined coastal aquifer in MATLAB. The specific steps in the modeling involve: 1) Discretization of the problem domain with nodes placed at a regular interval; 2) Specification of hydraulic conductivity, porosity, piezometric head gradient, depth of the aquifer, longitudinal and transverse dispersivity as input parameters; Variables are initialized and the initial head and concentration values are assigned; 3) Construction of support domain for all the nodes in the flow region. For every support domain single derivative of f_i with respect to x and z and also double derivative of f_i with respect to x, z, xz, and zx are found; 4) Assemblage of support domain matrix is assembled after applying boundary conditions into the global matrix; 5) Solution of the systems of equations for the state variables defining the head and solute concentration distribution in the flow domain. Initially for an assumed value of c, h is found which is used to compute c and the coupled equation is updated with the new c; 6) Progression of iterative solution till a user defined convergence is realized. Figure 1 shows the flowchart for the presently developed PCM based MFree model which was subsequently verified for its accuracy by examining its solution with two different well established research problems.

5. Case Studies

To investigate the effectiveness of the developed PCM model, the first case study considered is the well-recog- nized Henry’s problem [23] and the second work is a confined aquifer seawater intrusion problem reported by Pandit and Anand [3] .

5.1. Henry’s Problem

The domain of solution for this study is rectangular two-dimensional vertical cross-section of a confined coastal aquifer of length L and thickness d as shown schematically in Figure 2. The length-to-depth ratio for the problem domain is 2:1. The domain is bounded on top and bottom by impervious boundaries. A seawater body is there on right side of the problem domain and hence right boundary is referred as the seaward side and the left as landward side with freshwater inflow to the coastal aquifer.

The boundary conditions applied for the problem are relative concentration 0 at the fresh water inflow boundary and equal to 1 on the seaward boundary. Initial concentration in the aquifer domain is assumed to be 0. The mathematical representation of the applicable boundary conditions is listed in Table 1. A is the piezometric head gradient value substituted to take into account the inflow from the landward side which is taken as 151 to simulate the general inflow taken for the Henry’s problem.

To solve the Henry’s problem, after initial trials number of divisions considered in the x and z directions are 40 and 20 respectively. Consequently the number of nodes in the longitudinal and transverse direction is 41 and

21 respectively which makes the total number of nodes in the domain 861. The nodal distribution in the domain is as shown in Figure 3.

For the PCM model, q is presently considered as 0.98 adopted from Thomas et al. [15] after initial tuning of this parameter within a recommended range in literature varying from 0.5 to 1.03. Similarly α_c is considered as 8 for both the case studies. A radial support domain was found appropriate with a radius equal to 3 times the nodal distance d_c. For a d_c value of 0.05, the radius of the support domain works out to 0.15. Figure 4 shows the solute mass concentration variation represented as isochlors for the Henry’s problem. The “S” shape of the isochlors indicates the existence of seawater circulation at the bottom of the aquifer and a reduction in the landward extent of seawater penetration due to dispersion as demonstrated by Henry [23] . In Figure 5, a comparison of 0.5 isochlor positions obtained by the present model with the solutions of [20] [24] [25] and [26] is presented.

The results obtained using the PCM model are in near alignment with the known solutions. A marginal shift in the position of 0.5 isochlor towards the landward side is noticed when compared to the existing solutions. The PCM produced shape of the interface is “S” type trending similar to the shapes given by other researches. The “S” shape of the interface suggests that in the top aquifer region, the intrusion is much less as when compared to the bottom of the aquifer. The shift of seawater wedge towards landward side is attributed to the assumed values of α_T = 0.1 and d_l = 0.8, (i.e. the velocity dependent dispersion effects) while other parameters are kept same. It is to be noted that previous models assumes dispersion coefficient to be independent of the velocity which is a concept lacking in simulation of real flow conditions. It can also be understood from the comparison that position of the toe (considering 0.5 isochlor as interface) is better predicted in the MFree method than by FEM as the proximity of the MFree solution is nearer to the semi analytical solutions [23] when compared to the FEM solutions. From the velocity vectors plotted in Figure 6, it is observed that there is a clockwise circulation of seawater in the zone of transition along the interface which has earlier been demonstrated by other researchers [3] [16] .

5.1.1. Sensitivity Study of Meshfree Parameters

For the Henry’s problem, sensitivity of the model to number of nodes is performed to understand the effect of nodal distribution. The analysis is done by considering 210,861 and 3321 nodes in the domain. The 0.5 isochlors obtained for the Henry’s problem by varying the numbers of nodes 210, 861 and 3321 is shown in Figure 7.

It is encouraging that a larger nodal density involving 861 and 3321 nodes leads to near similar solutions. Presently simulation is performed on a Core i5 processor with speed of 3.1 GHz and 4 GB RAM. For a nodal specification of 210, 861 and 3321 nodes a computational time of 14.1322, 40.760, 248.944 seconds respectively was required for this configuration. This has suggested that after reaching a particular density of nodes in the domain, the accuracy of the result is not affected by increasing the number of nodes which is a positive outcome expected of a good numerical model. As expected when number of nodes is increased, the computational time also grows up considerably. It is thus understood from Figure 7 that increasing the number of nodes beyond certain limit does not improve the accuracy of the result further. This indicates a plausible strength of MFree simulation which is not much influenced by consideration of additional nodes in the flow domain and assumes a saving in simulation run time. FEM and FD solution scheme are known for nodal density dependence limiting their solution accuracy.

For the Henry’s problem, similar analysis was also carried out by varying the size of the support domain. Support domain is described as a function of number of times of the nodal distance. The simulation time required to model 861 nodal Henry’s problem for a support domain of 2 times the nodal distance was 33.03 seconds, where as three and four fold increase led to computational time increase to 40.76 seconds and 53.52 seconds respectively. The accuracy of the results however, remained almost the same in all the cases suggesting to keep the size of the support domain as minimum for the present problem without violating the Peclet number condition [22] . This analysis suggests robustness of the MFree simulation and recommends adopting a smaller support domain without sacrificing the accuracy of the solution.

5.2. Case Study 2

To further validate the applicability of the PCM model and to understand the pumping and recharge effects on seawater intrusion, a second case study from Pandit and Anand [3] for L/d ratio of 10:1 is considered. The rectangular two dimensional confined coastal aquifer with vertical cross section of length L and depth d is shown in Figure 8. The L/d ratio for the domain is 10:1 and data considered is as follows. A = 1000; K_ZZ = 1 porosity = 0.3, d_l = 1.0 and α_T = 0.1.

Relative concentration in the flow domain is 0 at the fresh water inflow boundary and 1 on the seaward boundary. Initial concentration in the aquifer domain is assumed 0. Boundary conditions for the problem are same as given in Table 1. For solving this problem, number of divisions considered after initial trials in the longitudinal x and transverse z directions were chosen 49 and 8 respectively. Consequently the number of nodes in x and z are 50 and 9 respectively. Thus the total number of nodes considered in the entire flow domain is 450. The nodal distribution in the problem domain is given in Figure 9. As stated earlier, the value of shape parameter q is considered as 0.98 and α_c is taken equal to 8. Based upon the experience of earlier case study, a radial support domain was considered with a radius equal to 4 times the nodal distance. There is a fresh water flow from left to

right which is horizontal at the fresh water face. The piezometric head gradient in the study is expressed in terms of A and for the present problem is taken as 1000. The domain is bounded by no flow boundary condition on top and bottom. Figure 10 shows the solute mass concentration variation measured as isochlors. The 0.5 isochlors obtained in present model is compared with that of Pandit and Anand [3] and Ukrande [26] which are plotted in Figure 11. Comparison of 0.5 isochlors shows reasonable agreement. MeshFree results are closer to the Pandit and Anand [3] and validate the developed meshfree model.

Figure 12 shows the velocity vector for this case study. As seen in the figure, in the seaward side the hydraulic head gradient is found to be oriented vertically upward almost against the gravitational force points vertically downward. There is a fresh water inflow due to which solute density is less at the bottom portion of the domain. Velocity vectors also indicate the circulation of seawater inside the domain as shown in Figure 12.

6. Sensitivity Study of Parameters

Saltwater intrusion primarily depends on the four non dimensional parameters presently such as: A―Gradient of peizometric surface; α_T―Transverse dispersivity; d_l―depth of aquifer and K_zz―the ratio of principle hydraulic conductivity tensors. Here an attempt is made to study the influence of these parameters in the salt water intrusion using the developed PCM model. The ranges of the parameters studied are listed in Table 2.

6.1. Effect of Gradient of Peizometric Surface

The seawater intrusion to the costal aquifer depends on fresh water gradient or the inflow rate from the landward side to the seaward side. The gradient of the piezometric surface varies from 0.01% to 0.5% in most cases as observed by Pandit and Anand [3] . For present study, piezometric surface is expressed in terms of A as can be seen in Equation (16). The bound for the parameter Ais therefore established as 10,000 to 200. To understand the effect of piezometric head gradient on the intrusion process, plots of 0.5 isochlors were compared for d_l = 1, K_zz = 1, α_T = 0.5 while varying the values of A from 200 to 10,000. Figure 13 shows the variation of 0.5 isochlor for a set of chosen values of A. It was found from the study that with the larger values of A, the lesser is the piezometric head gradient at the landward side and hence the seawater encroaches further into the aquifer. It is found that this effect is more prominent at the bottom of the aquifer as observed in Figure 13.

6.2. Effect of Longitudinal Dispersivity and Transverse Dispersivity

The longitudinal dispersivity values are found much larger in the field than obtained in laboratory experiments. After compiling data inferred from various field tests Gelhar [27] recommended that longitudinal dispersivity would be in the range of tens of meters to few hundreds of meters. Here in order to consider the problem in the non-dimensional approach both longitudinal and transverse dispersivity is divided by longitudinal dipersivity as expressed in Equation (6). To study the effect of dispersivity, the range of α_T which is presently the ratio of transverse to longitudinal dispersivity and therefore lower bound of 0.1 to an upper limit of 1.0 was adopted which is considered reasonable by Pandit and Anand [3] . Keeping the other parameters constant i.e. d_l = 1, K_zz = 1, A = 1000, α_T is varied from 0 to 1. Figure 14 shows the variation of 0.5 isochlor for the variation of α_T. It was observed that with an increase in the value of α_T, the seawater intrusion decreases as can be observed in Figure 14.

6.3. Effect of Depth of Aquifer

The depth of the aquifer is also another parameter that influences the encroachment of seawater into the fresh water. The depth of the aquifer normally varies in the range of 10 m to 500 m in the field problems. Here the non-dimensionalized parameter d_l which is the ratio of aquifer depth to its longitudinal dispersivity as defined in Equation (6) is adopted for present sensitivity study. Larger aquifers generally have larger dispersivities and hence even though the depth of the aquifers varies from tens of meters to few hundreds of meters, a value ranging between 0.01 to 20 is found very reasonable for d_l [3] . Figure 15 shows the 0.5 isochlors variation for different d_l values. The study regarding the depth of the aquifer proved interesting as it was found that initially when the d_l values are less, the 0.5 isochlors are found to be nearly vertical however with increase in the value of d_l the isochlor tends to take more like S shape with more encroachment at the bottom of the aquifer as can be observed in Figure 15.

It is due to the fact that for higher values of d_l the hydraulic effect causing cyclic flow becomes prominent as suggested by Cooper [28] whereas dispersion effects become negligible.

6.4. Effect of Hydraulic Conductivity

The property of hydraulic conductivity facilitates the movement of water through the porous media. The larger the pore size, the larger the hydraulic conductivity and vice versa generally holds true. To simplify the model, often a measured hydraulic conductivity value found by certain approach is used as uniform input value either for the entire model or for a particular grid. This method is suited for simulating the average response of the modeled system. In certain other approaches, spatial conductivity distributions from measured values are assigned to the model which is more tedious. The non-dimensional parameter K_zz is ratio of the hydraulic conductivity in transverse to the longitudinal direction as explained in Equation (5). It is commonly found that transverse hydraulic conductivity is lesser in magnitude than the longitudinal hydraulic conductivity. Presently, it is considered that K_zz vary from 0.1 to 1.0. To study the effect of anisotropy of the porous medium, all other parameter were kept constant and K_zz was assigned the values 0.1 and 1 in two respective cases. Figure 16 shows the 0.5 isochlor variation for different hydraulic conductivities. It was implicit that varying K_zz only produced a very small change in the location or shape of the isochlor as shown in Figure 16.

6.5. Effect of Pumping and Recharge

Further the simulation is carried for recharge effects. The recharge wells cause an increase of groundwater head in the flow domain. To simulate a steady state recharge well, the head imposed on the well node is kept slightly higher than the nearby area nodes where it is installed. The equivalent recharge rate can be calculated from the relevant parameters of the well-known drawdown (or built up for this case) equation. Presently head ratio at node 218 (8.775, 0.5) is kept 1.2. Figure 18 shows the 0.5 isochlors with and without recharge effects. It can be seen from the Figure 18 that 0.5 isochlor has moved towards seaward side i.e. toe of the 0.5 isochlor is shifted from 7.8 to 9.4 because of injection of fresh water from these wells when compared to the isochlors of the domain without a recharge well.

It was also seen by the sensitivity analysis that the movement of the dispersive interface also depends on the location of the recharge or pumping wells. Location of pumping wells and well head was also found to play an important role in the sustainable freshwater yield from the costal aquifers. The advantage of this technique of well simulation is that the influence of well location and its penetration on the movement of interface can be adequately analyzed for practical considerations.

7. Concluding Remarks

In the present study, a simulation model based on PCM MFree method is developed for a confined coastal aquifer. The model is based on the diffusive interface approach using point collocation based MFree method. Adequacy of the model was established by comparing the results with well-known problems from the previous research works. Velocity vectors show the presence of circulation of seawater in the transition zone of the flow region. Further a sensitivity study of various parameters was carried out. The sensitivity study shows that a decreasing piezometric head gradient and increasing depth of the aquifer has more pronounced effect on the encroachment of seawater intrusion into the landward side when compared to the variation in the transverse dispersivity or the ratio of principle hydraulic conductivity tensors. An examination of the effects of injection and pumping wells revealed that proper control measures such as involving installation of recharge wells can be effectively used to seaward pushing of the salt-fresh water interface. Study suggests that location of wells can be an important determining factor in the movement of the dispersive interface and it is possible to analyze the interface push back effect for various locations of injection wells. Such simulation model can be coupled with optimization models to suggest optimal location of wells for adequate management of coastal aquifers. A comparison of the PCM model results with the available FEM model solutions showed that PCM results are marginally better when both (PCM and FEM) were compared with the available semi analytic concentration distributions. PCM based MFree model has advantages in preprocessing of meshing/remeshing and when appropriate PCM parameters are identified, the model is found to be computationally more effective. The proposed method may be a better choice for seawater intrusion studies in coastal aquifers.