Next Article in Journal
Comparison of Flow-Dependent Controllers for Remote Real-Time Pressure Control in a Water Distribution System with Stochastic Consumption
Previous Article in Journal
Reflections: Contested Epistemologies on Large Dams and Mega-Hydraulic Development
Previous Article in Special Issue
Spectral Decomposition and a Waveform Cluster to Characterize Strongly Heterogeneous Paleokarst Reservoirs in the Tarim Basin, China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Civil Engineering, Indian Institute of Science, Bengaluru 560 012, India
2
Associate faculty, ICWaR, RBCCPS and IFCWS, Indian Institute of Science, Bengaluru 560 012, India
*
Author to whom correspondence should be addressed.
Water 2019, 11(3), 421; https://doi.org/10.3390/w11030421
Submission received: 10 December 2018 / Revised: 5 January 2019 / Accepted: 11 January 2019 / Published: 27 February 2019
(This article belongs to the Special Issue Heterogeneous Aquifer Modeling: Closing the Gap)

Abstract

:
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.

1. 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 (Kxx and Kyy), 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
  • To determine the relative sensitive aquifer parameters in modelling a layered coastal phreatic aquifer.
  • To estimate the anisotropic heterogeneous aquifer parameters by using 2D resistivity data.
  • To upscale the anisotropic hydraulic parameter from local measurements to aquifer scale.
  • 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.

2. 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 km2. 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 (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.
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.

3. 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.

3.1. Development of Conceptual Model

3.1.1. 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].

3.1.2. 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.

3.1.3. 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).

3.2. Sensitivity Analysis (SA)

3.2.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 (Kzz) is assumed to be 10% of Kxx. 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].

3.2.2. 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 (Kxx, Kyy, ε and αL). The mathematical expression to calculate the Sensitivity Index (SI) is given by Lenhart et al. [49] as:
SI = | Y 2 - Y l | Y 0 2 Δ x x 0
where Y is the vector of output state variables (hydraulic head and solute concentration) and x is an aquifer parameter (Kxx, Kyy, ε and αL). Y0 is the model output vector calculated with respect to the initial aquifer parameters value x0. 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 (Yn) from the initial output (Y0) is estimated using Root Mean Square Deviation (RMSD) which is given as:
RMSD = t = 1 T ( Y 0 , t Y n , t ) 2 T
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.

3.3. 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.
α L = 0.83 ( log 10 L ) 2.414
where L = field length from source (m) and αL = longitudinal dispersivity (m).
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.

3.4. 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.

4. Results

4.1. 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 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 Kxx (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 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.

4.2. 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 (R2) = 0.9244.
K = 1487.2 × 1 ρ + 5.255

4.3. 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.9 m. 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 layer-wise 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 Kyyy) 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 layer-wise Kzz field is 10% of Kxx 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.8 m/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, where the absolute difference in 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.

4.4. 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 = KyyKzz) 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.
The performance of the model is evaluated at the well Nos. 2 and 10 only (C values < 0.5 kg/m3 are not considered). The temporal bC and temporal RMSEC at well Nos. 2 and 10 for simulation period of 1095 days are tabulated in Table 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/m3, 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/m3 but the values less than 0.5 kg/m3 are ignored for RMSE and bC.

5. 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; Figure 8 and Figure 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.713 kg/m3 and −0.381 kg/m3, 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/m3 and 0.433 kg/m3, 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.

6. 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 R2 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.

Author Contributions

P.B.N. has contributed in conceptualization, methodology, validation, formal analysis, investigation, resources, data curation, visualization, writing—original draft preparation, writing—review and editing. M.S.M.K. is a supervisor of the research work and has contributed in conceptualization, validation, resources, writing-review and editing.

Funding

This research received no external funding.

Acknowledgments

The authors thank Johan Hoareau, Jean-Michel Vouillamoz (University Grenoble Alpes, IRD, CNRS, Grenoble INP, IGE, CS407000, 38058 Grenoble CEDEX 9, France) and Mahesha Amai (Department of Applied Mechanics and Hydraulics, NITK, Surathkal, Mangalore 575 025, India) for sharing the field measured data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ferguson, G.; Gleeson, T. Vulnerability of coastal aquifers to groundwater use and climate change. Nat. Clim. Chang. 2012, 2, 342–345. [Google Scholar] [CrossRef]
  2. Lee, C.-H.; Cheng, R.T.-S. On seawater encroachment in coastal aquifers. Water Resour. Res. 1974, 10, 1039–1043. [Google Scholar] [CrossRef]
  3. Segol, G.; Pinder, G.F. Transient simulation of saltwater intrusion in southeastern Florida. Water Resour. Res. 1976, 12, 65–70. [Google Scholar] [CrossRef]
  4. Gupta, D.A.; Yapa, P.N.D.D. Saltwater encroachment in an aquifer: A case study. Water Resour. Res. 1982, 18, 546–556. [Google Scholar] [CrossRef]
  5. Essaid, H.I. A multilayered sharp interface model of coupled freshwater and saltwater flow in coastal systems: Model development and application. Water Resour. Res. 1990, 26, 1431–1454. [Google Scholar] [CrossRef]
  6. Xue, Y.; Xie, C.; Wu, J.; Liu, P.; Wang, J.; Jiang, Q. A three-dimensional miscible transport model for seawater intrusion in China. Water Resour. Res. 1995, 31, 903–912. [Google Scholar] [CrossRef]
  7. Paniconi, C.; Khlaifi, I.; Lecca, G.; Giacomelli, A.; Tarhouni, J. A modelling study of seawater intrusion in the Korba Coastal Plain, Tunisia. Phys. Chem. Earth Part B Hydrol. Oceans Atmos. 2001, 26, 345–351. [Google Scholar] [CrossRef]
  8. Werner, A.D.; Gallagher, M.R. Characterisation of sea-water intrusion in the Pioneer Valley, Australia using hydrochemistry and three-dimensional numerical modelling. Hydrogeol. J. 2006, 14, 1452–1469. [Google Scholar] [CrossRef]
  9. Cobaner, M.; Yurtal, R.; Dogan, A.; Motz, L.H. Three dimensional simulation of seawater intrusion in coastal aquifers: A case study in the Goksu Deltaic Plain. J. Hydrol. 2012, 464–465, 262–280. [Google Scholar] [CrossRef]
  10. Chang, Y.; Hu, B.X.; Xu, Z.; Li, X.; Tong, J.; Chen, L.; Zhang, H.; Miao, J.; Liu, H.; Ma, Z. Numerical simulation of seawater intrusion to coastal aquifers and brine water/freshwater interaction in south coast of Laizhou Bay, China. J. Contam. Hydrol. 2018, 215, 1–10. [Google Scholar] [CrossRef] [PubMed]
  11. Priyanka, B.N.; Kumar, M.M.; Amai, M. Estimating anisotropic heterogeneous hydraulic conductivity and dispersivity in a layered coastal aquifer of Dakshina Kannada District, Karnataka. J. Hydrol. 2018, 565, 302–317. [Google Scholar] [CrossRef]
  12. Zamrsky, D.; Oude Essink, G.H.P.; Bierkens, M.F.P. Estimating the thickness of unconsolidated coastal aquifers along the global coastline. Earth Syst. Sci. Data Discuss. 2018, 10, 1–19. [Google Scholar] [CrossRef]
  13. Hamby, D.M. A review of techniques for parameter sensitivity analysis of environmental models. Environ. Monit. Assess. 1994, 32, 135–154. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Saltelli, A.; Tarantola, S.; Campolongo, F. Sensitivity analysis as an ingredient of modeling. Stat. Sci. 2000, 15, 377–395. [Google Scholar]
  15. Frey, C.H.; Patil, S.R. Identification and review of sensitivity analysis methods. Risk Anal. 2002, 22, 553–578. [Google Scholar] [CrossRef] [PubMed]
  16. Song, X.; Zhang, J.; Zhan, C.; Xuan, Y.; Ye, M.; Xu, C. Global sensitivity analysis in hydrological modeling: Review of concepts, methods, theoretical framework, and applications. J. Hydrol. 2015, 523, 739–757. [Google Scholar] [CrossRef] [Green Version]
  17. Pianosi, F.; Beven, K.; Freer, J.; Hall, J.W.; Rougier, J.; Stephenson, D.B.; Wagener, T. Sensitivity analysis of environmental models: A systematic review with practical workflow. Environ. Model. Softw. 2016, 79, 214–232. [Google Scholar] [CrossRef] [Green Version]
  18. Borgonovo, E.; Plischke, E. Sensitivity analysis: A review of recent advances. Eur. J. Oper. Res. 2016, 248, 869–887. [Google Scholar] [CrossRef]
  19. Ferretti, F.; Saltelli, A.; Tarantola, S. Trends in sensitivity analysis practice in the last decade. Sci. Total Environ. 2016, 568, 666–670. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Zimmerman, D.A.; Hanson, R.T.; Davis, P.A. A Comparison of Parameter Estimation and Sensitivity Analysis Techniques and Their Impact on the Uncertainty in Ground Water Flow Model Predictions (No. NUREG/CR-5522; SAND–90-0128); Nuclear Regulatory Commission: Washington, DC, USA; Division of High-Level Waste Management, Sandia National Labs: Albuquerque, NM, USA, 1991.
  21. Lathashri, U.A.; Mahesha, A. Predictive simulation of seawater intrusion in a tropical coastal aquifer. J. Environ. Eng. ASCE 2015, 142, 2–13. [Google Scholar] [CrossRef]
  22. Kolbe, T.; Marçais, J.; Thomas, Z.; Abbott, B.W.; de Dreuzy, J.R.; Rousseau-Gueutin, P.; Aquilina, L.; Labasque, T.; Pinay, G. Coupling 3D groundwater modeling with CFC-based age dating to classify local groundwater circulation in an unconfined crystalline aquifer. J. Hydrol. 2016, 543, 31–46. [Google Scholar] [CrossRef] [Green Version]
  23. Saha, G.C.; Li, J.; Thring, R.W. Understanding the effects of parameter uncertainty on temporal dynamics of groundwater-surface water interaction. Hydrology 2017, 4, 28. [Google Scholar] [CrossRef]
  24. Archie, G.E. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 1942, 146, 54–62. [Google Scholar] [CrossRef]
  25. Kelly, W.E. Geoelectric sounding for estimating aquifer hydraulic conductivity. Ground Water 1977, 15, 420–425. [Google Scholar] [CrossRef]
  26. Frohlich, R.K.; Kelly, W.E. Estimates of specific yield with the geoelectric resistivity method in glacial aquifers. J. Hydrol. 1988, 97, 33–44. [Google Scholar] [CrossRef]
  27. Heigold, P.C.; Gilkeson, R.H.; Cartwright, K.; Reed, P.C. Aquifer transmissivity from surficial electrical methods. Groundwater 1979, 17, 338–345. [Google Scholar] [CrossRef]
  28. Niwas, S.; Singhal, D.C. Aquifer transmissivity of porous media from resistivity data. J. Hydrol. 1985, 82, 143–153. [Google Scholar] [CrossRef]
  29. Frohlich, R.K.; Fisher, J.J.; Summerly, E. Electric-hydraulic conductivity correlation in fractured crystalline bedrock: Central Landfill, Rhode Island, USA. J. Appl. Geophys. 1996, 35, 249–259. [Google Scholar] [CrossRef]
  30. Chandra, S.; Ahmed, S.; Ram, A.; Dewandel, B. Estimation of hard rock aquifers hydraulic conductivity from geoelectrical measurements: A theoretical development with field application. J. Hydrol. 2008, 357, 218–227. [Google Scholar] [CrossRef]
  31. Tizro, A.T.; Voudouris, K.S.; Salehzade, M.; Mashayekhi, H. Hydrogeological framework and estimation of aquifer hydraulic parameters using geoelectrical data: A case study from West Iran. Hydrogeol. J. 2010, 18, 917–929. [Google Scholar] [CrossRef]
  32. Perdomo, S.; Ainchil, J.E.; Kruse, E. Hydraulic parameters estimation from well logging resistivity and geoelectrical measurements. J. Appl. Geophys. 2014, 105, 50–58. [Google Scholar] [CrossRef]
  33. Kazakis, N.; Pavlou, A.; Vargemezis, G.; Voudouris, K.S.; Soulios, G.; Pliakas, F.; Tsokas, G. Seawater intrusion mapping using electrical resistivity tomography and hydrochemical data. An application in the coastal area of eastern Thermaikos Gulf, Greece. Sci. Total Environ. 2016, 543, 373–387. [Google Scholar] [CrossRef] [PubMed]
  34. Wen, X.H.; Gómez-Hernández, J.J. Upscaling hydraulic conductivities in heterogeneous media: An overview. J. Hydrol. 1996, 183, ix–xxxii. [Google Scholar] [CrossRef]
  35. Vidstrand, P. Comparison of upscaling methods to estimate hydraulic conductivity. Groundwater 2001, 39, 401–407. [Google Scholar] [CrossRef]
  36. Gutjahr, A.L.; Gelhar, L.W.; Bakr, A.A.; MacMillan, J.R. Stochastic analysis of spatial variability in subsurface flows: 2. Evaluation and application. Water Resour. Res. 1978, 14, 953–959. [Google Scholar] [CrossRef]
  37. Tompson, A.F.; Ababou, R.; Gelhar, L.W. Implementation of the three-dimensional turning bands random field generator. Water Resour. Res. 1989, 25, 2227–2243. [Google Scholar] [CrossRef]
  38. Bellin, A.; Rubin, Y. HYDRO_GEN: A spatially distributed random field generator for correlated properties. Stoch. Hydrol. Hydraul. 1996, 10, 253–278. [Google Scholar] [CrossRef]
  39. Soraganvi, V.S.; Ababou, R.; Mohan Kumar, M.S. Analysis and upscaling of unsaturated flow through randomly heterogeneous soil. J. Hydrol. Eng. 2016, 22, 04016063. [Google Scholar] [CrossRef]
  40. Saldanha, J.P. Characterization of coastal aquifer system—A case study. MTech. Thesis, Department of Applied Mechanics and Hydraulics, National Institute of Technology, Karnataka, Surathkal, India, 2008. [Google Scholar]
  41. Hoareau, J. Utilisation d’une approche couplée hydrogéophysique pour l’étude des aquifers—Applications aux contextes de socle et côtier sableux. Ph.D. Thesis, Université Joseph-Fourrier-Grenoble 1, Grenoble, France, 2009. [Google Scholar]
  42. Vyshali. Studies on Saltwater Intrusion in the Coastal D.K District. Ph.D. Thesis, Department of Applied Mechanics and Hydraulics, National Institute of Technology, Karnataka, Surathkal, India, 2008. [Google Scholar]
  43. Honnanagoudar, S.A. Studies on Aquifer Characterization and Seawater Intrusion Vulnerability Assessment of Coastal Dakshina Kannada District, Karnataka. Ph.D. Thesis, Department of Applied Mechanics and Hydraulics, National Institute of Technology, Karnataka, Surathkal, India, 2014. [Google Scholar]
  44. Usha, A. Characterization of Large Diameter Wells in Shallow Coastal Unconfined Aquifers. MTech. Thesis, Department of Applied Mechanics and Hydraulics, National Institute of Technology, Karnataka, Surathkal, India, 2013. [Google Scholar]
  45. Diersch, H.-J.G. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer Science & Business Media: Berlin, Germany, 2014. [Google Scholar]
  46. Freeze, R.A.; Cherry, J.A. Groundwater; Prentice-Hall Inc.: Englewood Cliffs, NJ, USA, 1979; Volume 7632, p. 604. [Google Scholar]
  47. Udayakumar, G. Subsurface Barrier for Water Conservation in Lateritic Formations. Ph.D. Thesis, Department of Applied Mechanics and Hydraulics, National Institute of Technology, Karnataka, Surathkal, India, 2008. [Google Scholar]
  48. Gelhar, L.W.; Welty, C.; Rehfeldt, K.R. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 1992, 28, 1955–1974. [Google Scholar] [CrossRef]
  49. Lenhart, T.; Eckhardt, K.; Fohrer, N.; Frede, H.-G. Comparison of two different approaches of sensitivity analysis. Phys. Chem. Earth 2002, 27, 645–654. [Google Scholar] [CrossRef]
  50. Xu, M.; Eckstein, Y. Use of weighted least-squares method in evaluation of the relationship between dispersivity and field scale. Groundwater 1995, 33, 905–908. [Google Scholar] [CrossRef]
  51. Sanz, E.; Voss, C.I. Inverse modeling for seawater intrusion in coastal aquifers: Insights about parameter sensitivities, variances, correlations and estimation procedures derived from the Henry problem. Adv. Water Resour. 2006, 29, 439–457. [Google Scholar] [CrossRef]
  52. Ataie-Ashtiani, B.; Rajabi, M.M.; Ketabchi, H. Inverse modelling for freshwater lens in small islands: Kish Island, Persian Gulf. Hydrol. Process. 2013, 27, 2759–2773. [Google Scholar] [CrossRef]
  53. Mazáč, O.; Kelly, W.E.; Landa, I. A hydrogeophysical model for relations between electrical and hydraulic properties of aquifers. J. Hydrol. 1985, 79, 1–19. [Google Scholar] [CrossRef]
  54. Kelly, W.E.; Reiter, P.F. Influence of anisotropy on relations between electrical and hydraulic properties of aquifers. J. Hydrol. 1984, 74, 311–321. [Google Scholar] [CrossRef]
  55. Ataie-Ashtiani, B.; Volker, R.E.; Lockington, D.A. Tidal effects on seawater intrusion in unconfined aquifers. J. Hydrol. 1999, 216, 17–31. [Google Scholar] [CrossRef]
  56. Benini, L.; Antonellini, M.; Laghi, M.; Mollema, P.N. Assessment of water resources availability and groundwater salinization in future climate and land use change scenarios: A case study from a coastal drainage basin in Italy. Water Resour. Manag. 2016, 30, 731–745. [Google Scholar] [CrossRef]
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. Location of study area along with electrical resistivity profile, pumping test wells, Vertical Electrical Soundings (VES) survey points and observation wells.
Water 11 00421 g001
Figure 2. Schematic relation between aquifer parameter x with respective model output Y (Modified from Lenhart et al. [49]).
Figure 2. Schematic relation between aquifer parameter x with respective model output Y (Modified from Lenhart et al. [49]).
Water 11 00421 g002
Figure 3. RMSD for the various aquifer parameters over their respective valid range for (a) h outputs and (b) C outputs.
Figure 3. RMSD for the various aquifer parameters over their respective valid range for (a) h outputs and (b) C outputs.
Water 11 00421 g003
Figure 4. Linear regression analysis of hydraulic conductivity (K) and electrical conductivity (EC).
Figure 4. Linear regression analysis of hydraulic conductivity (K) and electrical conductivity (EC).
Water 11 00421 g004
Figure 5. Spatial variation of longitudinal dispersivity.
Figure 5. Spatial variation of longitudinal dispersivity.
Water 11 00421 g005
Figure 6. Variogram for generating layer-wise logn Kxx and logn Kyy fields.
Figure 6. Variogram for generating layer-wise logn Kxx and logn Kyy fields.
Water 11 00421 g006
Figure 7. Layer-wise heterogeneous Kxx and Kyy field.
Figure 7. Layer-wise heterogeneous Kxx and Kyy field.
Water 11 00421 g007
Figure 8. Spatial bh over 14 observation wells for the upscaled and transversely isotropic model.
Figure 8. Spatial bh over 14 observation wells for the upscaled and transversely isotropic model.
Water 11 00421 g008
Figure 9. Spatial bC at observation well No. 10 for the upscaled and transversely isotropic model.
Figure 9. Spatial bC at observation well No. 10 for the upscaled and transversely isotropic model.
Water 11 00421 g009
Table 1. Model input parameters.
Table 1. Model input parameters.
ParametersValues
Freshwater density1000 kg/m3
Seawater density1025 kg/m3
Molecular diffusion0.0000864 m2/d
Dynamic viscosity280,985.76 kg/m/yr.
Table 2. Layer wise flow and solute transport aquifer parameter ranges to perform sensitivity analysis (SA).
Table 2. Layer wise flow and solute transport aquifer parameter ranges to perform sensitivity analysis (SA).
LayerKxx (m/d)Kyy (m/d)ε (%)αL (m)
10.35–10
(5.175)
0.35–10
(5.175)
30–50
(40)
10–65
(37.5)
2 and 310–70
(40)
10–70
(40)
5–30
(17.5)
10–65
(37.5)
The parameter values in the bracket ( ) indicates initial values.
Table 3. Sensitivity Index (SI) and relative ranking of various input aquifer parameters based on h and C output.
Table 3. Sensitivity Index (SI) and relative ranking of various input aquifer parameters based on h and C output.
ParameterSIhRankSICRank
Kyy0.007110.5212
Kxx (with Kzz = 0.1 × Kxx)0.005820.4293
ε5.4 × 10−833 × 10−54
αL (with αT =0.1 × αL)--2.9391
Table 4. Estimated αL value at well points based on Xu and Eckstein [50] regression formula.
Table 4. Estimated αL value at well points based on Xu and Eckstein [50] regression formula.
Well NameField Length from Source (L), mLongitudinal Dispersivity (αL), mWell NameField Length from Source (L), mLongitudinal Dispersivity (αL), m
VES12006.205PW14258.555
VES22707.088PW45359.361
VES33457.860PW6>100011.772
VES41855.987PW7>100011.772
VES577110.730PW94808.976
VES64008.349PW103007.414
VES7>100011.772PW116309.960
VES865010.077PW12>100011.772
VES967010.191PW13>100011.772
VES101254.959PW143007.414
VES114508.751PW1583511.044
VES1275010.623
Table 5. Field measurements and estimated K from the regression model.
Table 5. Field measurements and estimated K from the regression model.
VES No.LocationAquifer Resistivity, ρ from VES (ohm-m)EC (EC = 1/ρ (mho/m))Estimated K from Equation (4) (m/d)
213°1′14.3″ N210.50.00475112.3056
74°47′58.8″ E
1413°01′14″ N200.640.00498412.65213
74°47′55″ E
513°00′24.2″ N29.60.03378455.3967
74°47′44.8″ E
713°00′30.9″ N159.60.00626614.5543
74°48′16.2″ E
813°00′57.82″ N232.3750.00430311.6419
74°47′34.87″ E
313°1′6.2″ N957.60.0010446.8047
74°47′23.1″ E
1513°00′59″ N1361.1780.0007356.3452
74°47′45″ E
Table 6. Input data for generating the logn K field.
Table 6. Input data for generating the logn K field.
ParametersValues
Length of study area in x-axis (Lx)3088 m
Length of study area in y-axis (Ly)4166.6 m
logn Kxx
ParametersLayer 1Layer 2Layer 3
Arithmetic mean (μx)3.4732.612.024
Variance (σx)0.0550.0710.0089
Correlation length (λx)40.02 m26.776m23.483 m
logn Kyy
ParametersLayer 1Layer 2Layer 3
Arithmetic mean (μy)3.4832.6452.035
Variance (σy)0.03730.0760.0079
Correlation length (λy)49.288 m357.68 m12.775 m
Table 7. Comparison of temporal performance between the upscaled model output and the transversely isotropic model output for h and C.
Table 7. Comparison of temporal performance between the upscaled model output and the transversely isotropic model output for h and C.
Model OutputPerformance Measure
Hydraulic Head (m)Solute Concentration (kg/m3)
Mean bhMean RMSEhbC at Well No. 2 RMSEC at Well No. 2bC at Well No. 10RMSEC at Well No. 10
Upscaled−1.812.26−0.7130.716−0.3810.433
Transversely isotropic−2.362.71−0.730.734−0.55350.609

Share and Cite

MDPI and ACS Style

B.N., P.; Mohan Kumar, M.S. Three-Dimensional Modelling of Heterogeneous Coastal Aquifer: Upscaling from Local Scale. Water 2019, 11, 421. https://doi.org/10.3390/w11030421

AMA Style

B.N. P, Mohan Kumar MS. Three-Dimensional Modelling of Heterogeneous Coastal Aquifer: Upscaling from Local Scale. Water. 2019; 11(3):421. https://doi.org/10.3390/w11030421

Chicago/Turabian Style

B.N., Priyanka, and M.S. Mohan Kumar. 2019. "Three-Dimensional Modelling of Heterogeneous Coastal Aquifer: Upscaling from Local Scale" Water 11, no. 3: 421. https://doi.org/10.3390/w11030421

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop