Three-Dimensional Modelling of Heterogeneous Coastal Aquifer : Upscaling from Local Scale

The aquifer heterogeneity is often simplified while conceptualizing numerical model due to lack of field data. Conducting field measurements to estimate all the parameters at the aquifer scale may not be feasible. Therefore, it is essential to determine the most significant parameters which require field characterization. For this purpose, the sensitivity analysis is performed on aquifer parameters, viz., anisotropic hydraulic conductivity, effective porosity and longitudinal dispersivity. The results of the sensitivity index and root mean square deviation indicated, that the longitudinal dispersivity and anisotropic hydraulic conductivity are the sensitive aquifer parameters to evaluate seawater intrusion in the study area. The sensitive parameters are further characterized at discrete points or at local scale by using regression analysis. The longitudinal dispersivity is estimated at discrete well points based on Xu and Eckstein regression formula. The anisotropic hydraulic conductivity is estimated based on established regression relationship between hydraulic conductivity and electrical resistivity with R2 of 0.924. The estimated hydraulic conductivity in x and y-direction are upscaled by considering the heterogeneous medium as statistically homogeneous at each layer. The upscaled model output is compared with the transversely isotropic model output. The bias error and root mean square error indicated that the upscaled model performed better than the transversely isotropic model. Thus, this investigation demonstrates the necessity of considering spatial heterogeneous parameters for effective modelling of the seawater intrusion in a layered coastal aquifer.


Introduction
The coastal aquifers provide fresh groundwater for more than 2 billion people worldwide [1].Groundwater stored in the coastal aquifers is susceptible to degradation due to its proximity to seawater, in combination with the intensive water demands.SeaWater Intrusion (SWI) is caused by prolonged changes in the coastal groundwater levels due to pumping, land-use change, climate variations/and sea-level fluctuations.The groundwater models provide a scientific and predictive tool for determining the appropriate solutions for SWI problems.Substantial research effort spanning a period of about 50 years has been devoted in understanding the groundwater flow and SWI at aquifer scale [2][3][4][5][6][7][8][9][10].The recent studies pointed out the heterogeneity, anisotropy and layering are often neglected or simplified while conceptualizing numerical model at both aquifer and global scales [11,12].This indicates that there is a wide gap in the knowledge about hydrogeology in modelling the SWI.The present study bridges the gap between hydrogeological characterization and SWI modelling by considering heterogeneity, anisotropy and layering.The development of the numerical models consists of several logical steps, one of which is assigning the appropriate aquifer parameters.Thus, the first step in modelling an aquifer is to determine the input parameters, which affects the model output significantly.In the present study, the sensitivity analysis is carried out to determine the relative sensitive parameter for which heterogeneity must be considered.A Sensitivity Analysis (SA) is the process of varying model input parameters over a reasonable range and observing the relative change in the model response.SA is a useful tool for model building and evaluation.There are several prior studies that have broadly reviewed the existing SA methods [13][14][15][16][17][18][19].Considerable research has been carried out to evaluate the parameter's impact on the groundwater model output [20][21][22][23].In the present study, SA is carried out for a range of aquifer parameters such as anisotropic hydraulic conductivity (K xx and K yy ), porosity (ε) and longitudinal dispersivity (α l ).The parameters which are relatively sensitive to the model output are considered for further studies on heterogeneity.
The heterogeneous aquifer parameters can be determined by conducting field experiments at discrete points (e.g., borehole logging, pumping test, Vertical Electrical Soundings/VES) or at local scale (e.g., electrical resistivity tomography/ERT).The pumping test or slug test is widely used to estimate the hydraulic parameters at discrete points.The other field measurements based on electrical resistivity are also used to estimate aquifer parameters [24][25][26].Several published field studies have proved that there exists a strong link between electrical resistivity and hydraulic parameters of an aquifer since the flow of fluid and electric current are governed by same physical and lithological attributes.The classical regression technique can be used to determine the hydraulic parameters [27][28][29][30][31][32][33].The earlier studies commonly used one dimensional (1D) resistivity test (e.g., VES) to determine the aquifer parameters.But these VES measures gives only vertical changes in the subsurface and cannot detect lateral changes [33].Therefore, in the present study to understand the vertical layering and lateral heterogeneity, ERT measurements are used.
The estimated aquifer parameters at the local scale (e.g., ERT) or discrete point measurements (e.g., VES) can be upscaled to aquifer scale to simulate the groundwater flow and solute transport.The aquifer parameters can be upscaled without considering the surrounding flow/transport field; such techniques are referred to as intrinsic upscaling or local upscaling [34,35].In this study, an intrinsic upscaling technique with stochastic field theory is adopted, in which parameter fields are generated considering the heterogeneous medium as statistically homogeneous [36][37][38][39].The earlier studies focused on determining the effective hydraulic parameter but in the present study anisotropic parameters are estimated for a layered coastal aquifer.The parameters at the aquifer scale are estimated based on the spatial correlation, where statistical parameters such as mean, variance and correlation length do not change over the scale.
The aim of the present investigation is to consider varying aquifer thickness and anisotropic heterogeneous aquifer parameters in modelling a three-dimensional (3D) coastal phreatic aquifer.The objectives of the paper are as follows 1.
To determine the relative sensitive aquifer parameters in modelling a layered coastal phreatic aquifer.2.
To estimate the anisotropic heterogeneous aquifer parameters by using 2D resistivity data.

3.
To upscale the anisotropic hydraulic parameter from local measurements to aquifer scale.4.
To simulate transient groundwater flow and solute transport for a 3D variable density conceptual model constrained with heterogeneity, anisotropy and layering.
The analysis gives an insight into the importance of considering the anisotropy and heterogeneity to develop a 3D transient groundwater flow and solute transport model.

Study Area
The area under investigation lies in the Dakshin Kannada region on the West coast of India, with an area extent of about 8 km 2 .The area is surrounded by the Arabian Sea on the west; River Pavanje on the north and north-eastern part as shown in Figure 1.The Pavanje River bends perpendicularly Water 2019, 11, 421 3 of 17 (13 • 1 54.49"N and 74 • 47 40.89"E) and flows north before joining the Arabian Sea.The Pavanje River is tidal in nature; thus, the adjoining aquifer gets contaminated by seawater for a considerable distance upstream during the non-monsoon period.Even though the fresh water requirement is partially met by surface water supply, there exists a greater dependency on groundwater resources.The educational institutes namely National Institute of Technology, Karnataka (NITK) and Srinivas group of Institutions are the major groundwater consumers from this aquifer.
Water 2018, 10, x FOR PEER REVIEW 3 of 18 is tidal in nature; thus, the adjoining aquifer gets contaminated by seawater for a considerable distance upstream during the non-monsoon period.Even though the fresh water requirement is partially met by surface water supply, there exists a greater dependency on groundwater resources.The educational institutes namely National Institute of Technology, Karnataka (NITK) and Srinivas group of Institutions are the major groundwater consumers from this aquifer.The area under consideration is moderately populated and cultivated.This region falls under tropical climate and each year is classified into four seasons according to India Meteorological Department (IMD).But this region experiences two well-marked seasons: the non-monsoon period from October-May and the monsoon period from June-September, since about 87% of annual average rainfall is experienced during the monsoon period.
The area under investigation has low-level laterite which is heterogeneous and discordant contact with the substratum.The lateritic formation underlaid by a bed of granitic gneiss bedrock of Archean age.The basin is predominantly an unconfined aquifer with depth ranging from 12-30 m.The present study area was chosen due to the availability of hydrological, hydrogeological and geophysical data.The IRD (Institut de Recherche pour le Développement, France) and NITK (National Institute of Technology, Karnataka) have conducted discrete and local field experiments in this area.The 2D electrical resistivity data are available at eighteen locations [40,41] pumping test data are available at fifteen locations [42][43][44] and VES test are available at twelve locations [42,43] as shown in Figure 1.

Methodology
The first step of the study is to determine the sensitive parameters for understanding the dynamics between groundwater flow and solute transport.For this analysis, a conceptual model is developed using FEFLOW by DHI-WASY GmbH, Berlin (Germany).Then the significant aquifer The area under consideration is moderately populated and cultivated.This region falls under tropical climate and each year is classified into four seasons according to India Meteorological Department (IMD).But this region experiences two well-marked seasons: the non-monsoon period from October-May and the monsoon period from June-September, since about 87% of annual average rainfall is experienced during the monsoon period.
The area under investigation has low-level laterite which is heterogeneous and discordant contact with the substratum.The lateritic formation underlaid by a bed of granitic gneiss bedrock of Archean age.The basin is predominantly an unconfined aquifer with depth ranging from 12-30 m.The present study area was chosen due to the availability of hydrological, hydrogeological and geophysical data.The IRD (Institut de Recherche pour le Développement, France) and NITK (National Institute of Technology, Karnataka) have conducted discrete and local field experiments in this area.The 2D electrical resistivity data are available at eighteen locations [40,41] pumping test data are available at fifteen locations [42][43][44] and VES test are available at twelve locations [42,43] as shown in Figure 1.

Methodology
The first step of the study is to determine the sensitive parameters for understanding the dynamics between groundwater flow and solute transport.For this analysis, a conceptual model is developed using FEFLOW by DHI-WASY GmbH, Berlin (Germany).Then the significant aquifer parameters are characterized at discrete locations and at local scale using geophysical data.The parameters estimated at the local scale is upscaled without considering the flow/or solute transport field.The upscaled model output is evaluated with respect to the measured hydraulic head and solute concentration.And the results are also compared with transversely isotropic model output to highlight the importance of considering heterogeneity in modelling a coastal aquifer.

Numerical Model
For the simulation of hydraulic head (h) and solute concentration (C), a numerical model FEFLOW (Finite Element subsurface FLOW) is used.The governing equations for the groundwater flow and solute transport are derived from the basic conservation principles for the mass of fluid, contaminants and linear momentum.For the present analysis, the Galerkin finite element method without up-winding is used.The Picard iteration method is used to treat the nonlinearities.The matrices are solved by Preconditioned Conjugate Gradient (PCG) method for symmetric matrix and Bi-Conjugate Gradient Square Stabilized (BiCGSTAB) for the unsymmetrical matrix.The other settings are set to default in FEFLOW [45].

Discretization
The domain is discretized with six nodal triangular prisms and with fine discretization at the boundaries, streams and wells.The 3D triangular meshes are generated by using the triangulation code built in FEFLOW.The phreatic conceptual model is vertically stratified into three layers.The Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) of 30 m × 30 m data is used to define the top elevation of the model.The bottom topography and thickness of each layer are based on the depth of observation wells, 1D and 2D resistivity data.

Boundary Conditions
The coastal boundary is defined as a Dirichlet boundary condition with the equivalent head and constant solute concentration.The river boundary is defined as a transfer boundary condition and the data required for this boundary condition are river stage, river water solute concentration, hydraulic conductivity and thickness of the river bed.The remaining part of the domain is defined as specified head and specified concentration based on groundwater level and solute concentration contours.The top boundary is defined as a specified flux condition and phreatic.The bottom of the model is defined with Neumann boundary condition.
For the developed conceptual model, the source and sink are assigned with uniform recharge rate and transient groundwater draft, respectively.The parameters tabulated above are used throughout the model (Table 1).

Model Inputs for SA
For SA, aquifer parameters are considered as layered heterogeneous (but spatially homogeneous) based on the soil type and earlier field investigations.The respective aquifer parameter ranges are tabulated in Table 2.The hydraulic conductivity (K) is considered as anisotropic and vertical hydraulic conductivity (K zz ) is assumed to be 10% of K xx .The porosity (ε) ranges are considered based on earlier investigations and literature [46,47].The longitudinal dispersivity (α L ) is taken from the literature [48].The transverse dispersivity (α T ) is assigned as 1/10th of α L [21].
Table 2. Layer wise flow and solute transport aquifer parameter ranges to perform sensitivity analysis (SA).The parameter values in the bracket ( ) indicates initial values.

SA Method
In the present analysis, sensitivity is expressed by a dimensionless index [49], which is calculated as the ratio between the relative change in the model response to the relative change of each aquifer parameter (K xx , K yy , ε and α L ).The mathematical expression to calculate the Sensitivity Index (SI) is given by Lenhart et al. [49] as: where Y is the vector of output state variables (hydraulic head and solute concentration) and x is an aquifer parameter (K xx , K yy , ε and α L ).Y 0 is the model output vector calculated with respect to the initial aquifer parameters value x 0 .Each parameter is varied ±∆x and the other parameters remain at the initial value.To assess the calculated sensitivity indices, relative ranks are assigned to each parameter (Figure 2).The deviation of model output (Y n ) from the initial output (Y 0 ) is estimated using Root Mean Square Deviation (RMSD) which is given as: where t is the number of points at which state variables are estimated.RMSD is used to understand how the model responds to the variation in each parameter over a valid range.

Regression Analysis
The SA results indicate that the anisotropic hydraulic conductivity and longitudinal dispersivity are the significant parameters in 3D modelling of the coastal phreatic aquifer.Thus, heterogeneity analysis is carried out for these specific parameters.The field measurements on dispersivity values were not available.Therefore, regression formula (Equation ( 3)) derived by Xu and Eckstein [50] is used to estimated heterogeneous longitudinal dispersivity.
The available field measurements such as VES, pumping test and ERT are used to estimate anisotropic K.The ERT profiles are considered such that these profiles are located approximately 500 m away from the coastal line and not affected by salinity.The VES survey and pumping test carried out at same discrete locations are used to establish the regression equation between K and electrical conductivity (EC).From the linear regression model, K values are estimated locally (at electrical resistivity profile) in x and y-axis.The deviation of model output (Yn) from the initial output (Y0) is estimated using Root Mean Square Deviation (RMSD) which is given as: where t is the number of points at which state variables are estimated.RMSD is used to understand how the model responds to the variation in each parameter over a valid range.

Regression Analysis
The SA results indicate that the anisotropic hydraulic conductivity and longitudinal dispersivity are the significant parameters in 3D modelling of the coastal phreatic aquifer.Thus, heterogeneity analysis is carried out for these specific parameters.The field measurements on dispersivity values were not available.Therefore, regression formula (Equation (3)) derived by Xu and Eckstein [50] is used to estimated heterogeneous longitudinal dispersivity.The available field measurements such as VES, pumping test and ERT are used to estimate anisotropic K.The ERT profiles are considered such that these profiles are located approximately 500 m away from the coastal line and not affected by salinity.The VES survey and pumping test carried out at same discrete locations are used to establish the regression equation between K and electrical conductivity (EC).From the linear regression model, K values are estimated locally (at electrical resistivity profile) in x and y-axis.

Intrinsic Upscaling of Aquifer Parameters
Ordinary Kriging is used to estimate heterogeneous α L from discrete well locations to aquifer scale.From the established linear regression model, K values are estimated locally (at electrical resistivity profile) in x and y-axis.The direction perpendicular to the coastal line is considered as x-axis and direction parallel to the coastal line is considered as the y-axis.
The hydraulic conductivity is a random space function and is characterized statistically by their spatial moments.A spatially correlated random field is generated by considering Gaussian covariance function with grid dimension in x and y directions, arithmetic mean, variance and correlation length of each layer.For further details on a spatially correlated random field generator, the reader can refer Bellin and Rubin [38].
The generated K fields in both x and y-axis are assigned as input to the conceptual model.The transient groundwater flow and solute transport simulation results are compared with the field measured state variables.The upscaled anisotropic heterogeneous K model results are also compared with transversely isotropic (pumping test data) results.

Sensitivity Analysis
The steady-state conceptual phreatic 3D model is simulated by varying only one parameter with respect to the initial aquifer parameter for every simulation.To calculate SI based on hydraulic head Water 2019, 11, 421 7 of 17 output, h values at all the nodes are considered but in case of the C, nodes at which C values ≥ 500 mg/L are considered, because values > 500 mg/L is most important (According to the World Health Organization).Table 3 shows the ranks assigned relatively based on SI values to the aquifer parameter with respect to h and C. The aquifer parameter with rank 1 can be identified as a highly sensitive parameter.
From Table 3, it is observed that α L and anisotropic K are the important parameters for modelling SWI.The hydraulic conductivity in the y-axis (i.e., K parallel to the coast) is more significant than K xx (i.e., K normal to the coast) in the present study as the Pavanje River which flows normally to the coast contributes to the value of h as well as for C. As also indicated in Table 3 that the h and C outputs are insensitive to effective porosity.Thus, this parameter does not significantly affect groundwater flow and SWI.
To quantify the amount of output deviation from the initial model, RMSD is used.To compute RMSD, h values and C values at all the nodes are considered.Figure 3a,b show the RMSD for h and C outputs respectively for the various aquifer parameters.
From the Figure 3a, it can be observed that K yy for lower limit has higher RMSD.The RMSD was not computed for the lower limit of K xx because a lower limit of K xx = 0.35 m/d (K zz = 0.035 m/d), leads to a build-up of groundwater above ground level.The RMSD for hydraulic conductivity in the x-axis (with K zz ) and the y-axis are almost identical as the range considered for the two parameters are same.
From the Figure 3b, it can be observed that α L (considering α T ) show significant deviation from 2.5-3.5 mg/L over its range.This indicates α L is the most critical parameter in 3D modelling coupled SWI problem.It can also be noted that unlike Figure 3a the RMSD of K yy shows greater deviation than K xx .This may be because of solute transport from the tidal Pavanje River (especially for the first quartile).From the SA, it can be concluded that the effective porosity is an insensitive parameter.The longitudinal dispersivity and anisotropic hydraulic conductivity are the significant parameters to be considered for heterogeneity studies.

Determination of Aquifer Parameters at the Local Scale
The longitudinal dispersivity at well points is determined by Xu and Eckstein [50] regression formula (Equation ( 3)).The distance of the well point from the sea or river (source) is considered as L, to estimate α L in Equation ( 3).The values of estimated α L are tabulated in Table 4.The points located beyond 1 km from the source are rounded up to 1 km values because Equation ( 3) is valid up to 1 km only.
The hydraulic conductivity measured from pumping test and respective aquifer resistivity from VES is used to establish a relationship (Table 5).Figure 4 shows positive correlation between K and EC and thus K is negatively correlated with aquifer resistivity (ρ).The linear regression analysis of EC and K gave the following relationship (Figure 4 and Equation ( 4)) with a coefficient of determination (R  From the Figure 3a, it can be observed that Kyy for lower limit has higher RMSD.The RMSD was not computed for the lower limit of Kxx because a lower limit of Kxx = 0.35 m/d (Kzz = 0.035 m/d), leads to a build-up of groundwater above ground level.The RMSD for hydraulic conductivity in the x-axis (with Kzz) and the y-axis are almost identical as the range considered for the two parameters are same.
From the Figure 3b, it can be observed that αL (considering αT) show significant deviation from 2.5-3.5 mg/L over its range.This indicates αL is the most critical parameter in 3D modelling coupled SWI problem.It can also be noted that unlike Figure 3a the RMSD of Kyy shows greater deviation than Kxx.This may be because of solute transport from the tidal Pavanje River (especially for the first quartile).From the SA, it can be concluded that the effective porosity is an insensitive parameter.The longitudinal dispersivity and anisotropic hydraulic conductivity are the significant parameters to be considered for heterogeneity studies.

Determination of Aquifer Parameters at the Local Scale
The longitudinal dispersivity at well points is determined by Xu and Eckstein [50] regression formula (Equation ( 3)).The distance of the well point from the sea or river (source) is considered as L, to estimate αL in Equation ( 3).The values of estimated αL are tabulated in Table 4.The points located

Intrinsic Upscaling of Aquifer Parameters
The αL values from discrete points are upscaled to aquifer scale by using ordinary Kriging.The spatially varying longitudinal dispersivity is shown in Figure 5 and it is seen that the αL values vary between 4.96-11.9m.The VES Nos. 1, 4 and 10 which are located near the tidal river have αL < 6.2 m.

Intrinsic Upscaling of Aquifer Parameters
The α L values from discrete points are upscaled to aquifer scale by using ordinary Kriging.The spatially varying longitudinal dispersivity is shown in Figure 5 and it is seen that the α L values vary between 4.96-11.9m.The VES Nos. 1, 4 and 10 which are located near the tidal river have α L < 6.2 m.
The K xx (perpendicular to coast) and K yy (parallel to coast) values are estimated at the ERT profiles based on the regression relationship.The variogram and data required for generating layer-wise log n K xx and log n K yy fields are shown and tabulated in Figure 6 and Table 6, respectively.The correlation length of K in x-axis in all the layers range from 23.5-40.1 m.And the correlation length in y-axis in layer 2 is greater than 355 m and in layer 3 it is less than 13 m.This indicates the necessity of considering layering and spatial heterogeneity.It can also be noticed that the correlation length for log n K yy (λ y ) in layer 2 is relatively high compared to other layers.This indicates that there is no change in K yy value, for correlation length up to 357 m.The distribution of single realization layer-wise K xx and K yy fields are interpolated to get the anisotropic K field at the aquifer scale (Figure 7).The layer-wise K zz field is 10% of K xx field.

Intrinsic Upscaling of Aquifer Parameters
The αL values from discrete points are upscaled to aquifer scale by using ordinary Kriging.The spatially varying longitudinal dispersivity is shown in Figure 5 and it is seen that the αL values vary between 4.96-11.9m.The VES Nos. 1, 4 and 10 which are located near the tidal river have αL < 6.2 m.The Kxx (perpendicular to coast) and Kyy (parallel to coast) values are estimated at the ERT profiles based on the regression relationship.The variogram and data required for generating layerwise logn Kxx and logn Kyy fields are shown and tabulated in Figure 6 and Table 6, respectively.The correlation length of K in x-axis in all the layers range from 23.5-40.1 m.And the correlation length in y-axis in layer 2 is greater than 355 m and in layer 3 it is less than 13 m.This indicates the necessity of considering layering and spatial heterogeneity.It can also be noticed that the correlation length for logn Kyy (λy) in layer 2 is relatively high compared to other layers.This indicates that there is no change in Kyy value, for correlation length up to 357 m.The distribution of single realization layer-wise Kxx and Kyy fields are interpolated to get the anisotropic K field at the aquifer scale (Figure 7).The layerwise Kzz field is 10% of Kxx field.

Numerical Model Simulation Results
The upscaled anisotropic K field and longitudinal dispersivity are used as input to the conceptual model to simulate hydraulic head (h) and solute concentration (C) for a period of 1095 day with stress period of 1 day.The K values obtained from the pumping test is interpolated for the entire area and this transversely isotropic (i.e., K xx = K yy = K zz ) model results are compared with upscaled model numerical simulation results.The performance of both the models are evaluated by bias error (b) and root mean square error (RMSE).
The mean temporal b h and mean temporal RMSE h over 14 observation wells for simulation period of 1095 days are tabulated in Table 7.The mean temporal b h in upscaled model output is 23.3% less than transversely isotropic model output and mean temporal RMSE h is 16.6% lesser.Figure 8 illustrates the spatial b h comparison between upscaled and transversely isotropic model output over a period of 1095 day.The spatial bias error of hydraulic head varies from 0-3.9 m, where the high error is during the monsoon period (Figure 8).The b h for upscaled model output is less than transversely isotropic model output throughout the transient simulation period.It can be concluded with these results that the h output from the upscaled conceptual model is better than the transversely isotropic model.

Numerical Model Simulation Results
The upscaled anisotropic K field and longitudinal dispersivity are used as input to the conceptual model to simulate hydraulic head (h) and solute concentration (C) for a period of 1095 day with stress period of 1 day.The K values obtained from the pumping test is interpolated for the entire area and this transversely isotropic (i.e., Kxx = Kyy ≠ Kzz) model results are compared with upscaled model numerical simulation results.The performance of both the models are evaluated by bias error (b) and root mean square error (RMSE).
The mean temporal bh and mean temporal RMSEh over 14 observation wells for simulation period of 1095 days are tabulated in Table 7.The mean temporal bh in upscaled model output is 23.3% less than transversely isotropic model output and mean temporal RMSEh is 16.6% lesser.Figure 8 illustrates the spatial bh comparison between upscaled and transversely isotropic model output over a period of 1095 day.The spatial bias error of hydraulic head varies from 0-3.9 m, where the high error is during the monsoon period (Figure 8).The bh for upscaled model output is less than transversely isotropic model output throughout the transient simulation period.It can be concluded with these results that the h output from the upscaled conceptual model is better than the transversely isotropic model.7. The temporal bC in upscaled model output at well Nos. 2 and 10 are 2.33% and 31.165%less than transversely isotropic model output, respectively.The temporal RMSEC at well Nos. 2 and 10 are 2.44% and 28.9% lesser, respectively.Figure 9 illustrates the spatial bc at well No. 10 between upscaled and transversely isotropic model output.The spatial bh value varies from 0-0.75 kg/m 3 , which is lesser than the transversely isotropic model output.Therefore, it can be concluded that the C output from the upscaled conceptual model

Discussion
One of the main conceptual difficulties to model a 3D SWI process is assigning heterogeneous aquifer parameters.The SA is carried out to determine the sensitive parameters which require further heterogeneous investigations.The SA result on this coupled flow and transport problem showed that the effective porosity does not affect the model output significantly (Table 3 and Figure 3).The analysis highlighted that the groundwater flow and solute transport to and from the tidal river is more significant (based on the relative ranking).From the analysis, it can be concluded that the longitudinal dispersivity and anisotropic hydraulic conductivity are the sensitive parameters.The result of this analysis is reasoned as other similar studies carried out in coastal aquifers also indicate these parameters as significant [51,52].
The sensitive parameters are characterized at discrete well locations based on the regression analysis.The regression analysis between hydraulic conductivity and aquifer resistivity showed inverse correlation (Figure 4).This inverse correlation indicates the presence of the geological formation of the same sediment group [53].One of the important factors which influence the relationship between K and aquifer resistivity is anisotropy caused by layering [53,54].Therefore, in the present study three layers are considered based on their resistivity range.
An innovative approach based on the directions of ERT profiles is used to estimate anisotropic hydraulic conductivity.The intrinsic upscaling technique is used to estimate the anisotropic hydraulic conductivity at the aquifer scale because of lacking data.The upscaling technique generates the spatially correlated random field at each layer of Kxx and Kyy.The range of correlation length between the layer (i.e., 23-40 m for Kxx and 12-358 m for Kyy) and small values of correlation length (e.g., 12.8 m at layer 3 of Kyy) indicates the necessity of considering layering and spatial heterogeneity, respectively.
The upscaled model output is compared with a transversely isotropic model output which is developed from pumping test data (Table 7; Figures 8 and 9).The upscaled model performed better than the transversely isotropic model.This is because the spatial heterogeneity in the transversely isotropic model is restricted, due to the lack of available pumping test data.This comparison also highlights the importance of considering heterogeneities in modelling a coastal aquifer.However, the upscaled model output compared with observed data show high mean temporal bias error and RMSE of the hydraulic head, with bh and RMSE of −1.81 m and 2.26 m, respectively.The bias error of solute concentration at wells 2 and 10 are −0.713kg/m 3 and −0.381 kg/m 3 , respectively.The negative sign in bias error indicates that the upscaled model underestimates the values.The RMSE of solute concentration at wells 2 and 10 are 0.716 kg/m 3 and 0.433 kg/m 3 , respectively.
The upscaled model was unable to simulate accurate state variables, due to single K field generation or could be due to phenomena that are neglected.The phenomena such as tidal influence [55]/land-use change [56] which were not considered can be addressed in future studies.The main drawback of the study is considering the single K field as input to the numerical model.Multiple

Discussion
One of the main conceptual difficulties to model a 3D SWI process is assigning heterogeneous aquifer parameters.The SA is carried out to determine the sensitive parameters which require further heterogeneous investigations.The SA result on this coupled flow and transport problem showed that the effective porosity does not affect the model output significantly (Table 3 and Figure 3).The analysis highlighted that the groundwater flow and solute transport to and from the tidal river is more significant (based on the relative ranking).From the analysis, it can be concluded that the longitudinal dispersivity and anisotropic hydraulic conductivity are the sensitive parameters.The result of this analysis is reasoned as other similar studies carried out in coastal aquifers also indicate these parameters as significant [51,52].
The sensitive parameters are characterized at discrete well locations based on the regression analysis.The regression analysis between hydraulic conductivity and aquifer resistivity showed inverse correlation (Figure 4).This inverse correlation indicates the presence of the geological formation of the same sediment group [53].One of the important factors which influence the relationship between K and aquifer resistivity is anisotropy caused by layering [53,54].Therefore, in the present study three layers are considered based on their resistivity range.
An innovative approach based on the directions of ERT profiles is used to estimate anisotropic hydraulic conductivity.The intrinsic upscaling technique is used to estimate the anisotropic hydraulic conductivity at the aquifer scale because of lacking data.The upscaling technique generates the spatially correlated random field at each layer of K xx and K yy .The range of correlation length between the layer (i.e., 23-40 m for K xx and 12-358 m for K yy ) and small values of correlation length (e.g., 12.8 m at layer 3 of K yy ) indicates the necessity of considering layering and spatial heterogeneity, respectively.
The upscaled model output is compared with a transversely isotropic model output which is developed from pumping test data (Table 7; Figures 8 and 9).The upscaled model performed better than the transversely isotropic model.This is because the spatial heterogeneity in the transversely isotropic model is restricted, due to the lack of available pumping test data.This comparison also highlights the importance of considering heterogeneities in modelling a coastal aquifer.However, the upscaled model output compared with observed data show high mean temporal bias error and RMSE of the hydraulic head, with b h and RMSE of −1.81 m and 2.26 m, respectively.The bias error of solute concentration at wells 2 and 10 are −0.713kg/m 3 and −0.381 kg/m 3 , respectively.The negative sign in bias error indicates that the upscaled model underestimates the values.The RMSE of solute concentration at wells 2 and 10 are 0.716 kg/m 3 and 0.433 kg/m 3 , respectively.
The upscaled model was unable to simulate accurate state variables, due to single K field generation or could be due to phenomena that are neglected.The phenomena such as tidal influence [55]/land-use change [56] which were not considered can be addressed in future studies.
The main drawback of the study is considering the single K field as input to the numerical model.Multiple realizations can be performed to estimate anisotropic K field.And based on the performance evaluation the best suitable anisotropic K field can be used for further predictive studies.

Conclusions
In the present investigation, hydrogeological characterization is carried out by using field data and intrinsic upscaling technique.The sensitivity analysis is performed to determine the significant aquifer parameters in modelling SWI.The SA results indicate that the hydraulic conductivity in the y-direction and longitudinal dispersivity are the significant parameter in modelling groundwater flow and SWI in a coastal phreatic aquifer, respectively.The hydraulic conductivity in the y-direction (i.e., K parallel to the coast) is more significant in the present case due to the presence of the river.This illustrates the necessity of considering anisotropic hydraulic conductivity in numerical simulations.The SA results also demonstrate that for the present SWI problem, model outputs are insensitive to effective porosity.Thus, sensitive parameters, that is, anisotropic hydraulic conductivity and longitudinal dispersivity are investigated further at the aquifer scale.
The VES and pumping test data are used to establish a relationship between hydraulic conductivity and electrical conductivity.The inverse relationship between hydraulic conductivity and electrical resistivity with R 2 of 0.9244 is used to determine the local hydraulic conductivity in x and y-directions.Due to the absence of field measurement on dispersivity, Xu and Eckstein [50] regression formula is used to estimate longitudinal dispersivity at well points.The locally estimated longitudinal dispersivity are upscaled at aquifer scale by using ordinary Kriging.The layer-wise anisotropic hydraulic conductivities are upscaled based on the stochastic field theory in x and y-direction.
The upscaled anisotropic heterogeneous aquifer parameters are used as input to develop a transient 3D conceptual model.The upscaled 3D model output for both state variables (h and C) are compared with a transversely isotropic model output which is developed from pumping test data.The mean temporal and spatial bias error and RMSE of the transversely isotropic model is greater than the upscaled model.Therefore, it can be concluded that the upscaled conceptual 3D model is better than the transversely isotropic model.This study is a building block towards aquifer scale 3D modelling in a coastal anisotropic heterogeneous porous media.

Figure 1 .
Figure 1.Location of study area along with electrical resistivity profile, pumping test wells, Vertical Electrical Soundings (VES) survey points and observation wells.

Figure 1 .
Figure 1.Location of study area along with electrical resistivity profile, pumping test wells, Vertical Electrical Soundings (VES) survey points and observation wells.

18 Figure 2 .
Figure 2. Schematic relation between aquifer parameter x with respective model output Y (Modified from Lenhart et al. [49]).

Figure 2 .
Figure 2. Schematic relation between aquifer parameter x with respective model output Y (Modified from Lenhart et al. [49]).

Figure 3 .
Figure 3. RMSD for the various aquifer parameters over their respective valid range for (a) h outputs and (b) C outputs.

Figure 3 .
Figure 3. RMSD for the various aquifer parameters over their respective valid range for (a) h outputs and (b) C outputs.

Figure 6 .
Figure 6.Variogram for generating layer-wise log n K xx and log n K yy fields.

Figure 7 .
Figure 7. Layer-wise heterogeneous Kxx and Kyy field.From Figure 7, it can be observed that Kxx and Kyy values range from 5.4-84.5 m/d and 5.6-74.8m/d, respectively.The K values in both x and y-axis decrease with depth indicating the top layer are relatively permeable.The Kxx field in layer 1 is relatively heterogeneous compared to other layers,

Figure 7 .
Figure 7. Layer-wise heterogeneous K xx and K yy field.From Figure 7, it can be observed that K xx and K yy values range from 5.4-84.5 m/d and 5.6-74.8m/d, respectively.The K values in both x and y-axis decrease with depth indicating the top layer are relatively permeable.The K xx field in layer 1 is relatively heterogeneous compared to other layers, where the absolute difference in K xx is 69.2 m/d.The K yy fields in layer 2 are relatively less

Figure 8 .
Figure 8. Spatial bh over 14 observation wells for the upscaled and transversely isotropic model.

Figure 8 .
Figure 8. Spatial b h over 14 observation wells for the upscaled and transversely isotropic model.The performance of the model is evaluated at the well Nos. 2 and 10 only (C values < 0.5 kg/m 3 are not considered).The temporal b C and temporal RMSE C at well Nos. 2 and 10 for simulation period of 1095 days are tabulated in Table7.The temporal b C in upscaled model output at well Nos. 2 and 10 are 2.33% and 31.165%less than transversely isotropic model output, respectively.The temporal RMSE C at well Nos. 2 and 10 are 2.44% and 28.9% lesser, respectively.Figure9illustrates the spatial b c at well No. 10 between upscaled and transversely isotropic model output.The spatial b h value varies from 0-0.75 kg/m 3 , which is lesser than the transversely isotropic model output.Therefore, it can be concluded that the C output from the upscaled conceptual model is better than the transversely Figure 9 illustrates the spatial b c at well No. 10 between upscaled and transversely isotropic model output.The spatial b h value varies from 0-0.75 kg/m 3 , which is lesser than the transversely isotropic model output.Therefore, it can be concluded that the C output from the upscaled conceptual model is better than the transversely isotropic model.The temporal observation C values at well No. 2 vary between 0.1-0.6 kg/m 3 but the values less than 0.5 kg/m 3 are ignored for RMSE and b C .Water 2018, 10, x FOR PEER REVIEW 14 of 18 is better than the transversely isotropic model.The temporal observation C values at well No. 2 vary between 0.1-0.6 kg/m 3 but the values less than 0.5 kg/m 3 are ignored for RMSE and bC.

Figure 9 .
Figure 9. Spatial bC at observation well No. 10 for the upscaled and transversely isotropic model.

Figure 9 .
Figure 9. Spatial b C at observation well No. 10 for the upscaled and transversely isotropic model.

Table 3 .
Sensitivity Index (SI) and relative ranking of various input aquifer parameters based on h and C output.

Table 5 .
Field measurements and estimated K from the regression model.

Table 6 .
Input data for generating the logn K field.

Table 6 .
Input data for generating the log n K field.

Table 7 .
Comparison of temporal performance between the upscaled model output and the transversely isotropic model output for h and C.

Output Performance Measure Hydraulic Head (m) Solute Concentration (kg/m 3 ) Mean b h Mean RMSE h b C at Well No. 2 RMSE C at Well No. 2 b C at Well No. 10 RMSE C at Well No. 10
Kxx is 69.2 m/d.The Kyy fields in layer 2 are relatively less heterogeneous compared to other layers because the Kyy values range from 10.5-12.1 m/d.The absolute difference in K values for layer 3 in both x and y-axis is less than 4.8 m/d, this indicates the presence of similar geological formation in both axes.

Table 7 .
Comparison of temporal performance between the upscaled model output and the transversely isotropic model output for h and C.