Modelling Groundwater Hydraulics to Design a Groundwater Level Monitoring Network for Sustainable Management of Fresh Groundwater Lens in Lower Indus Basin, Pakistan

The over-extraction of groundwater from thin fresh groundwater lenses is a threat to the livelihood of farmers in the Lower Indus Basin (LIB). It is essential to monitor and regulate this pumping to sustain fresh groundwater lenses. In this study, we applied a modelling approach in combination with geostatistical analysis to identify the critical locations to monitor the groundwater levels for sustaining fresh groundwater in the LIB. Our approach included four steps: (i) simulating temporal heads using a calibrated hydrogeological model; (ii) sampling monitoring locations using a hexagonal pattern of sampling; (iii) applying principal component analysis (PCA) of the temporal head observations, and selecting high scoring locations from the PCA; and (iv) minimizing the observation points to represent the water level contours. The calibrated model was able to replicate the hydro-dynamic behavior of the study area, with a root mean square of 0.95 and an absolute residual mean of 0.74 m. The hexagonal pattern of spatial sampling resulted in a 195 point network, but PCA reduced this network to 135 points and contour classification reduced it even further to 59 points. The 195, 135, and 59 point networks represented the water levels with average standard errors of 0.098, 0.318, and 0.610 m, respectively. Long-term simulations with increased pumping showed that the water levels would best be assessed by 195 monitoring points, although 135 and 59 points would represent the depleting area but would not capture the water logging area.


Introduction
Groundwater in the Lower Indus Basin (LIB) originates from the river system, which has flowed through the valley since late Tertiary times. The fresh water occurs in lenses near the Indus River or near the major canals. These freshwater lenses are underlain by dense saline water. This is due to the hydrogeological conditions in the LIB, which cause rapid mineralization of the groundwater. The largest freshwater lens occurs along the Indus River, and its quality deteriorates away from the river [1]. Inequity in surface water distribution and increasing irrigation intensities have led to an increase in groundwater dependence on fresh groundwater lenses. Government subsidies for also be coupled with optimization algorithms to design monitoring networks. Janardhanan et al. [12] used flow modelling coupled with the probabilistic framework to locate optimal locations for the piezometers. Similarly, Maqsood et al. [13] used numerical modelling and statistical analysis for designing a groundwater monitoring network. In the Indus Basin, and particularly the Central Indus Basin, groundwater modelling studies have been carried out to study the groundwater dynamics at regional and sub-regional levels. Ashraf and Ahmad [22] developed a regional flow model using the FEFLOW model to study the impact of adverse climatic conditions and plausible groundwater extractions on regional groundwater levels in the Upper Chaj Doab of the Indus Basin. Khan et al. [23] calibrated the MODFLOW [24] and MT3D [25] models for the Rechna Doab of the Indus Basin to assess the impact of dry conditions and increased extraction on regional groundwater systems. Similarly, Punthakey et al. [26] used modelling to explore options for equitable surface water and groundwater distribution to maximize the crop production and manage salinization in the Rechna Doab. In the Lower Indus Basin, groundwater modelling studies have been conducted at small spatial scales to study the groundwater hydraulics and hydro salinity behavior. Kori et al. [27] calibrated the MODFLOW and MT3D models to optimize the operation of selected tube wells in Sindh. Qureshi et al. [28] studied hydro salinity behavior using the calibrated MT3D model of the Kunner II irrigation distributary, Sindh. Chandio et al. [29] calibrated finite element groundwater models to investigate the impact of pumping and canal water levels on the waterlogging in the Khairpur district, Sindh. All of these studies were site-specific and lacked representation of the regional response of the aquifer. In this study, we calibrated a sub-regional MODFLOW model that provided a base tool for regional-level investigation in the freshwater zone in the LIB. This regional-level effort to numerically model freshwater lenses will be fundamental for informing evidence-based decisions for monitoring and managing groundwater systems. The approach presented in this paper is replicable in an irrigated area, where groundwater is used conjunctively with surface water for irrigation.
In this study, we used an approach that integrates hydrogeological flow simulation outputs with statistical analysis for finding critical monitoring locations in freshwater lenses in the LIB. The first aspect was hydrogeological characterization, so we characterized the aquifer properties using lithological logs, and calibrated a flow model for the study area. The calibrated model was simulated to generate monthly temporal heads. Then we applied a hexagonal pattern of sampling for continuous spatial functions on the model grid to select initial monitoring locations. Olea (1984) [30] stated that the semi-variance and drift are inherent in continuous spatial functions and are therefore unmanageable factors. The manageable factors in sampling continuous spatial functions are sampling pattern and sampling density. It is not economically feasible to have dense sampling. The only remaining factor which can be managed is the pattern of sampling. A hexagonal pattern of sampling has the least average standard error in the case of continuous spatial functions [30]. Average standard error is a measure of the sampling efficiency of continuous spatial functions. Since groundwater is also a continuous spatial function, a hexagonal pattern of sampling was selected for the sampling of observation grid cells in the flow simulation model. The second aspect, accounting for temporal variations, was considered using PCA. PCA is a dimensionality reduction technique, which uses an orthogonal transformation to eliminate the correlation in the observations of variables by converting them into linearly uncorrelated variables called principal components [3,4]. The third aspect, the representation of potentiometric surface with minimum wells, was addressed through groundwater potentiometric surface classification. The potentiometric head surface used for the head classification was obtained through an inverse distance weighting (IDW) technique. IDW is a technique of spatial interpolation that estimates the unknown values based on the known values by giving more weight to nearer points and less weight to farther points [31]. This interpolation technique does not abruptly alter the interpolated surfaces, so, in our study, this technique was applied to obtain the potentiometric head surface. Finally, we rationalized the identified locations by comparing standard errors between identified location and base case locations, which were based on 1 × 1 km model grids. We simulated the future levels on the historical climatic cycles to determine whether the monitoring network could replicate the future hotspots for management.

Study Area
The study area includes the northern Rohri canal command area of Sindh province, Pakistan, which lies in the Lower Indus Basin (LIB). Two irrigation divisions, Dad and Moro, were studied within the administrative boundaries of two districts, Shaheed Benazirabad and Nausheroferoze ( Figure 1). The climatic regime for the study area is classed as subtropical continental low land type, which is typically characterized by hot summers and mild winters. The hottest months are May and June, when temperatures can exceed 45 • C; and the coldest month is January, when temperatures can fall below 10 • C. The mean annual rainfall varies from 100 to 200 mm, with most falling during the monsoon season, i.e., July, August, and September. The occurrence of fresh water lenses in this area of the LIB is unique, as most of the groundwater in the basin is not suitable for drinking and irrigation. In this area, at the tail reaches the community is entirely dependent on groundwater and groundwater is extensively used for irrigation. The conservation of these fresh water lenses is essential for sustaining the livelihood of the communities and ecosystems in this area.

Hydrogeology
The unconsolidated deposits in the area are the alluvial fill of Quaternary age. The Indus River forms the western boundary, and the Rohri canal flows at the east. Floodplain, streambed, and meander belt deposits are dominant in the area. Lithologically, the area consists mainly of sands of various grades with silt and clay down to several hundred meters [32]. Sand is predominant, and is highly transmissive and constitutes the potential aquifer in the area [32]. The groundwater flow in the region is divided into two directions, with the hydrological divide seen near the Rohri canal. At the west side of the Rohri canal, water flows to or from the Indus River, and to the east of the Rohri canal, the flow is towards the southeast. The aquifer is fairly transmissive, thick, and suitable for the installation of high capacity wells. Groundwater occurs under water-table conditions and depth to water varies from 1.5 to 6 m [33]. Freshwater lenses are found along the river, and total dissolved solids in the freshwater zones are up to 1000 mg/L, and deteriorate gradually with depth and distance away from the river [34].

Aquifer Properties
Aquifer properties were calculated from 41 lithological logs [32]. These logs were digitized in the geo-modelling software Rockworks 16, and a 3D geological model was created ( Figure 2). Initial values of permeability (i.e., kh, Kv = 0.33 kh) and specific yield (i.e., Sy) were assigned based on the log material. There was a continuous thick aquifer in the study area, so the layer thicknesses were defined based on the tube well depth. The first layer was considered to be 35 m thick. The bottom of the aquifer (i.e., the bottom of layer 2) was assigned based on the thick layer of clay or sequence of alternating clay, and sand occurred at the bottom of bore logs, which was in the range of 180 to 200 m. An initial estimate of hydraulic conductivities (kh) showed that the kh values for the first 35 m (layer 1) ranged from 25 to 32 m/d, with a spatial average of 28 m/d. The distribution of higher kh values was greater on the north side of the model. The kh values for the second layer (layer 2), which had variable depth, ranged from 24 to 37 m/d, with a spatial average of 33 m/d. The distribution of high kh values for layer 2 was in most of the study area. The spatial distribution of specific yield (Sy) for layer 1 was reasonably high with an average value of 0.15. The gridded range of Sy was narrow, ranging from 0.14 to 0.16. The range of this layer's specific storage (Ss) varied from 1.2 × 10 −3 to 6.66 × 10 −3 . For layer 2, the Ss values were multiplied by the layer thickness to calculate the storage coefficient. Initial estimates are consistent with the literature; generally, it is estimated that the hydraulic conductivity in LIB ranges between 1 and 20 m/d and specific yields range from 5 to 15% [34][35][36][37].

2.4.Groundwater Flow Simulation
Flow simulation was performed via calibration of the groundwater flow model (i.e., MODFLOW 2005 [24]) from October 2010 until April 2014. A monthly water balance was quantified. MODFLOW 2005 uses a continuity equation for water balance as shown in Equation (1) and a finite difference scheme to solve it numerically [24].

Model Conceptualization and Boundary Conditions
The model was conceptualized as a two-layer version as shown in Figure 3. The first layer included the river, canal, and private tube wells. The second layer covered the deep tube wells that drain the saline water from the area. The model simulated a monthly stress period with a minimum of three time steps for model convergence. The model was discretized on a 1000 × 1000 m grid. The model boundaries were set up by investigating the depth to water level contours ( Figure 1). On the north side, we followed a 42 m contour line to provide a general head boundary. On the east side, a water divide was observed at which groundwater flows southeast and southwest, respectively. We demarcated the boundary at the water divide and assigned no flow to the east side of the model domain. On the west side is the Indus River, and we designated the river as the western boundary of the model domain. The south of the model was assigned as a flux boundary. Figure 4 shows the model active area and model boundary conditions that cover the Dad division and a portion of the Moro division.

Groundwater Flow Simulation
Flow simulation was performed via calibration of the groundwater flow model (i.e., MODFLOW 2005 [24]) from October 2010 until April 2014. A monthly water balance was quantified. MODFLOW 2005 uses a continuity equation for water balance as shown in Equation (1) and a finite difference scheme to solve it numerically [24].

Model Conceptualization and Boundary Conditions
The model was conceptualized as a two-layer version as shown in Figure 3. The first layer included the river, canal, and private tube wells. The second layer covered the deep tube wells that drain the saline water from the area. The model simulated a monthly stress period with a minimum of three time steps for model convergence. The model was discretized on a 1000 × 1000 m grid. The model boundaries were set up by investigating the depth to water level contours ( Figure 1). On the north side, we followed a 42 m contour line to provide a general head boundary. On the east side, a water divide was observed at which groundwater flows southeast and southwest, respectively. We demarcated the boundary at the water divide and assigned no flow to the east side of the model domain. On the west side is the Indus River, and we designated the river as the western boundary of the model domain.
The south of the model was assigned as a flux boundary. Figure 4 shows the model active area and model boundary conditions that cover the Dad division and a portion of the Moro division.

Evapotranspiration
Major sink/source terms in our study area were evapotranspiration (ET) from shallow water tables, recharge from rainfall and irrigation return flows, river leakages, and pumping. We used a Simplified Surface Energy Balance product [38] for an initial estimate of actual ET, and extinction depth was calculated via soil type in each cell. Data was downloaded for each stress period in the study area, and a temporal composite was used as a model input. In June and July, the maximum actual evapotranspiration (Eta) exceeded 300 mm per month, and varied spatially, with highest values Appl. Sci. 2020, 10, 5200 7 of 19 occurring at farm fields. The extinction depths assigned to each log were adopted from Shah et al. [39]. The interpreted extinction depths for each log were then gridded and a value for each model cell was obtained. The spatial variation of the extinction depths ranged from a low of 1.46 m to a high of 4.63 m, with a mean of 2.81 m and a standard deviation of 0.45 m.

Evapotranspiration
Major sink/source terms in our study area were evapotranspiration (ET) from shallow water tables, recharge from rainfall and irrigation return flows, river leakages, and pumping. We used a

Evapotranspiration
Major sink/source terms in our study area were evapotranspiration (ET) from shallow water tables, recharge from rainfall and irrigation return flows, river leakages, and pumping. We used a

Recharge
Recharge was calculated with Equation (2). Two components (recharge from rainfall and recharge from irrigations flows) were estimated, the latter of which was further divided into two subcomponents: recharge from canal irrigation and recharge from pumping irrigation. Rainfall recharge was estimated from Gridded Climate Products (GCPs) obtained from the Climate Forecast System Reanalysis (CFSR) dataset [40]. The GCP dataset was re-gridded over the groundwater model domain to represent the maximum potential rainfall recharge amount in each cell. Basin and furrow irrigation methods are widely used for supplying water to wheat, cotton, vegetable, fruit, and fodder crop fields during Rabi and Kharif seasons in the study area. The irrigation water allowance of the canal command area was determined in the modeled area on a monthly basis ( Figure 5). Since surface irrigation methods, which are prevalent in the area, were 40 to 60 percent efficient, it was estimated that 50 to 60 percent of irrigation water in the modeled area contributed to the aquifer as recharge. The amount was adjusted spatially and temporally during the calibration process. Based on the water allowance, each grid cell was assigned returns from irrigation on the monthly stress period. Initial recharge values were then assigned based on Equation (2), and initial factors f 1 , f 2 , and f 3 of 0.6, 0.5, and 0.5 were used respectively, and were later tuned in the calibration process to estimate net recharge to groundwater.
where f 1 , f 2 , f 3 = fitting parameters; R i,j = rainfall in cell; Ic i,j = return from canal irrigation; Ip i,j = return from pumping irrigation.
depth was calculated via soil type in each cell. Data was downloaded for each stress period in the study area, and a temporal composite was used as a model input. In June and July, the maximum actual evapotranspiration (Eta) exceeded 300 mm per month, and varied spatially, with highest values occurring at farm fields. The extinction depths assigned to each log were adopted from Shah et al. [39]. The interpreted extinction depths for each log were then gridded and a value for each model cell was obtained. The spatial variation of the extinction depths ranged from a low of 1.46 m to a high of 4.63 m, with a mean of 2.81 m and a standard deviation of 0.45 m.

Recharge
Recharge was calculated with Equation (2). Two components (recharge from rainfall and recharge from irrigations flows) were estimated, the latter of which was further divided into two subcomponents: recharge from canal irrigation and recharge from pumping irrigation. Rainfall recharge was estimated from Gridded Climate Products (GCPs) obtained from the Climate Forecast System Reanalysis (CFSR) dataset [40]. The GCP dataset was re-gridded over the groundwater model domain to represent the maximum potential rainfall recharge amount in each cell. Basin and furrow irrigation methods are widely used for supplying water to wheat, cotton, vegetable, fruit, and fodder crop fields during Rabi and Kharif seasons in the study area. The irrigation water allowance of the canal command area was determined in the modeled area on a monthly basis ( Figure 5). Since surface irrigation methods, which are prevalent in the area, were 40 to 60 percent efficient, it was estimated that 50 to 60 percent of irrigation water in the modeled area contributed to the aquifer as recharge. The amount was adjusted spatially and temporally during the calibration process. Based on the water allowance, each grid cell was assigned returns from irrigation on the monthly stress period. Initial recharge values were then assigned based on Equation (2), and initial factors f1, f2, and f3 of 0.6, 0.5, and 0.5 were used respectively, and were later tuned in the calibration process to estimate net recharge to groundwater.

River Leakage
River leakage was calculated with Equation (3). The canal network, consisting of main and branch canals, distributaries, minors, sub-minors, and water courses, were too extensive to model at a regional scale, where grid scales were relatively coarse. To rationalize the modeling, a river package was used to simulate the main canals, branch canals, and distributaries. Initially, the hydraulic conductivity of the river cells was considered to be in the range of 0.002 to 0.004 m/d, with an initial bed thickness of 1 m, which was further adjusted during the calibration process.

River Leakage
River leakage was calculated with Equation (3). The canal network, consisting of main and branch canals, distributaries, minors, sub-minors, and water courses, were too extensive to model at a regional scale, where grid scales were relatively coarse. To rationalize the modeling, a river package was used to simulate the main canals, branch canals, and distributaries. Initially, the hydraulic conductivity of the river cells was considered to be in the range of 0.002 to 0.004 m/d, with an initial bed thickness of 1 m, which was further adjusted during the calibration process.
where k i,j = river bed conductivity; l l,j = length of river in cell; w i,j = width of river in the cell; m i,j = river bed thickness; hriv i,j = head in river; haq i,j = head in aquifer.

Pumping
Pumping density was estimated via field surveys and government tube well records, and pumping density was used to calculate pumping in each cell using Equation (4), which was then given as the input to the model. An initial survey estimated that pumping was concentrated at the tail reaches, with a density of five (5) wells per square kilometer. Most of the private tube wells have a capacity of 1-2 cusecs and pump for an average of 8 h per day. At the tail reaches, groundwater was the sole source of irrigation, and at the mid reaches, farmers used groundwater during periods of surface water shortages.
where Cf = unit conversion for m 3 d .

Model Calibration
The groundwater flow model was calibrated by adjusting aquifer parameters k, S (Equation (1)), stress factors f 1 , f 2 , f 3 (Equation (2)) and river parameters Kriv ij , m ij (Equation (3)) using a trial and error approach. We reduced the root mean square value and absolute residual mean between simulated and observed head values of 35 SCARP's Monitoring Organization (SMO) observation points from post-monsoon 2010 until pre-monsoon 2014. Multiple simulations were performed to achieve the best suited combination of aquifer and stress factors that represented the temporal and spatial water levels with minimum absolute residual mean and RMS of greater than 0.9.

Hexagonal Pattern of Sampling
The grid cells of the groundwater flow simulation model were considered as observation wells in the area. The active model area was 3415 km 2 , where the groundwater was available in the shallow layer of 35 m. A total of 195 hexagons of 17 km 2 each were made. Since the model area was discretized by forming squares of 1000 × 1000 m, the hexagons were made accordingly (Figure 6a). The central grid cell of each hexagon was considered as an observation well unless the central grid cell was not in either the river, a canal, or a distributary. In such cases, the adjacent grid cell was considered to be an observation well.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 20 most representative of the overall groundwater head variance, thereby minimizing the required number of monitoring wells. Varimax rotation was then applied to further reduce the representing wells. Varimax rotation additionally cleaned the principal components. This maximized the sum of the variances of the square loadings and preserved the invariance of the variables. The application of varimax rotation further optimized the principle components and this became evident in the results. Indeed, prior to undertaking varimax rotation, a total of 148 wells were optimal in the first two components, and these two components explained 98.489% of the variance. After varimax rotation, the first two components explained the same amount of variance but with 135 optimal wells. In this study, only wells that had a factor loading matrix score equal to or greater than 0.90 were prioritized in the rotated loading factor matrix. The threshold of 0.90 was selected because we had to remove the correlation in the data. This threshold only retained those wells where the groundwater head fluctuations were high [3].

Potentiometric Head Classification
An inverse distance weighting (IDW) technique was used to obtain the potentiometric surface of the active model area. The point data used for interpolation was the average head data of prioritized wells over three and half years, based on the results of PCA. The logic for taking the average point data for obtaining flow patterns in the area was based on the cyclic trend of groundwater recovery. The potentiometric surface was then used to classify the groundwater head. The groundwater head variations in the study area were in the range of 22.89 to 41.37 m. A total of 19 classes were made with groundwater head variation of 0.97 m within each class. Three observation wells in each class were selected as prioritized wells. The aim of the study was to rationally reduce the number of observation wells, so we had two choices: either select the three wells with the highest scores in the PCA analysis or select three wells based on the spatial coverage of the groundwater head class. Each choice was tried, and the average standard error was determined for each network. The average standard error when selection was based on spatial coverage was lower than when it was based on the highest PCA analysis scores. Thus, the selection of wells in each class was based on the spatial locations of the wells. Figure 6b shows an example of how the wells can be selected using a head classification method. The contour classes (different colors) are of 0.97 m head variations. The wells shown in yellow are wells prioritized during the PCA stage. Of these yellow wells, red circled wells were selected as prioritized wells based on their spatial location in the groundwater head class (one well at the start of class, one well in the middle of the class, and one well at the end of the class). In addition, a few additional wells were selected based on the interpolation results, such as a few pockets of water which were part of a given class but were not continuous within that class. If there were already three or fewer wells in any class, they were left unaltered.

Principal Component Analysis
The calibrated model was run for 42 months and provided groundwater head data in the sampled wells on a monthly basis. Hence, a total of 42 temporal observations from 195 observation wells were available to inform the further prioritization of observation wells in the study area. PCA was run on all 195 observation wells and the 42 monthly temporal observations. The wells were considered as variable (columns) and temporal head values were considered as observations (rows). Since the PCA identifies axes which record the highest variation in the data, these axes are called principal components and are orthogonal to each other. The first component records the maximum variance of the data and the second component records the maximum variance in the remaining data and so forth. The total number of principal components can be equal to the number of variables. The intention of applying PCA in this study was to identify the locations of monitoring wells that were most representative of the overall groundwater head variance, thereby minimizing the required number of monitoring wells. Varimax rotation was then applied to further reduce the representing wells. Varimax rotation additionally cleaned the principal components. This maximized the sum of the variances of the square loadings and preserved the invariance of the variables. The application of varimax rotation further optimized the principle components and this became evident in the results. Indeed, prior to undertaking varimax rotation, a total of 148 wells were optimal in the first two components, and these two components explained 98.489% of the variance. After varimax rotation, the first two components explained the same amount of variance but with 135 optimal wells. In this study, only wells that had a factor loading matrix score equal to or greater than 0.90 were prioritized in the rotated loading factor matrix. The threshold of 0.90 was selected because we had to remove the correlation in the data. This threshold only retained those wells where the groundwater head fluctuations were high [3].

Potentiometric Head Classification
An inverse distance weighting (IDW) technique was used to obtain the potentiometric surface of the active model area. The point data used for interpolation was the average head data of prioritized wells over three and half years, based on the results of PCA. The logic for taking the average point data for obtaining flow patterns in the area was based on the cyclic trend of groundwater recovery. The potentiometric surface was then used to classify the groundwater head. The groundwater head variations in the study area were in the range of 22.89 to 41.37 m. A total of 19 classes were made with groundwater head variation of 0.97 m within each class. Three observation wells in each class were selected as prioritized wells. The aim of the study was to rationally reduce the number of observation wells, so we had two choices: either select the three wells with the highest scores in the PCA analysis or select three wells based on the spatial coverage of the groundwater head class. Each choice was tried, and the average standard error was determined for each network. The average standard error when selection was based on spatial coverage was lower than when it was based on the highest PCA analysis scores. Thus, the selection of wells in each class was based on the spatial locations of the wells. Figure 6b shows an example of how the wells can be selected using a head classification method. The contour classes (different colors) are of 0.97 m head variations. The wells shown in yellow are wells prioritized during the PCA stage. Of these yellow wells, red circled wells were selected as prioritized wells based on their spatial location in the groundwater head class (one well at the start of class, one well in the middle of the class, and one well at the end of the class). In addition, a few additional wells were selected based on the interpolation results, such as a few pockets of water which were part of a given class but were not continuous within that class. If there were already three or fewer wells in any class, they were left unaltered.

Groundwater Flow Simulation Model
Model calibration was performed for 42 stress periods from October 2010 through April 2014. Head measurements in October 2010 (post-monsoon) were taken as the initial head condition. Observed heads were arranged for post-and pre-monsoon season for each year from 2010 until 2014. A total of eight temporal points were used in the calibration. Figure 7 shows the calibration results of observed heads measured at over 35 wells across the model domain. A root mean square value of 0.95 and an absolute mean error of 0.74 m were obtained, which shows acceptable calibration of the simulated heads over the model domain. Figure 8 shows the residuals (i.e., observed minus simulated heads) at the end of the model calibration. Most piezometers in the study area (i.e., 23 piezometers) showed residuals in the range of 0-1 m, while 12 piezometers showed residuals in the range of 1-2 m. The model's performance was lower in the regions near the Indus River than in other regions of the model domain. Hydrographs of the simulated head are presented in Figure 9, showing that the model was able to replicate the response to external stress, and that simulated trends followed the temporal observed trend.

Hexagonal Pattern of Sampling
The hexagonal pattern of spatial sampling resulted in a total of 195 hexagons of 17 km 2 each. Thus, a total of 195 observation locations were selected from 3415 (base case scenario) through the hexagonal pattern of sampling. The potentiometric head surfaces were generated using the groundwater head values in these 195 wells for pre-monsoon (October) and post-monsoon (April) seasons using an IDW interpolation. These potentiometric head surfaces were compared with the base case. The standard error raster was generated by subtracting the potentiometric head surface obtained by the 195 observation wells from the potentiometric head surfaces obtained by the base case (i.e., 3415 observations). The mean error was found to be 0.098 m with a standard deviation of 0.58 m ( Table 1). These results can be considered satisfactory when the reduction in observation wells is more than 94% (195 out of 3415 wells).

Groundwater Flow Simulation Model
Model calibration was performed for 42 stress periods from October 2010 through April 2014. Head measurements in October 2010 (post-monsoon) were taken as the initial head condition. Observed heads were arranged for post-and pre-monsoon season for each year from 2010 until 2014. A total of eight temporal points were used in the calibration. Figure 7 shows the calibration results of observed heads measured at over 35 wells across the model domain. A root mean square value of 0.95 and an absolute mean error of 0.74 m were obtained, which shows acceptable calibration of the simulated heads over the model domain. Figure 8 shows the residuals (i.e., observed minus simulated heads) at the end of the model calibration. Most piezometers in the study area (i.e., 23 piezometers) showed residuals in the range of 0-1 m, while 12 piezometers showed residuals in the range of 1-2 m. The model's performance was lower in the regions near the Indus River than in other regions of the model domain. Hydrographs of the simulated head are presented in Figure 9, showing that the model was able to replicate the response to external stress, and that simulated trends followed the temporal observed trend.

Hexagonal Pattern of Sampling
The hexagonal pattern of spatial sampling resulted in a total of 195 hexagons of 17 km 2 each. Thus, a total of 195 observation locations were selected from 3415 (base case scenario) through the hexagonal pattern of sampling. The potentiometric head surfaces were generated using the groundwater head values in these 195 wells for pre-monsoon (October) and post-monsoon (April) seasons using an IDW interpolation. These potentiometric head surfaces were compared with the base case. The standard error raster was generated by subtracting the potentiometric head surface obtained by the 195 observation wells from the potentiometric head surfaces obtained by the base case (i.e., 3415 observations). The mean error was found to be 0.098 m with a standard deviation of 0.58 m ( Table 1). These results can be considered satisfactory when the reduction in observation wells is more than 94% (195 out of 3415 wells).

Principal Component Analysis
The potentiometric head data of 195 observation wells (selected from sampling) were analyzed using PCA in IBM SPSS (version 22.0) software. The major output of the PCA results were two matrices i.e., "total variance explained" and "rotated loading factor" matrix. The total variance explained matrix provided statistics about the variance of the data explained by each component. In our case, the first two components of PCA explained more than 98% of the variation of the groundwater head in the area. Hence, only the first two components were selected for further analysis. The rotated loading factor matrix gave a score between 0 and 1 based on the importance of a given well in explaining the groundwater head variation in the study area. In our case, 135 of the 195 wells scored greater than or equal to 0.90 in the rotated loading factor matrix. It was found that these 135 wells could represent the study area with an average error of 0.318 m and standard deviation of 0.95 m ( Table 1). The number of observation wells was reduced from 195 to 135 (i.e., about a 31% reduction) with an insignificant loss of information. This reduction is less compared to that shown in previous literature, in which the reduction in randomly installed monitoring wells was more than 70% [3,23]. The smaller reduction in observation wells occurred because the network was spatially optimized using a hexagonal pattern before the application of PCA.

Principal Component Analysis
The potentiometric head data of 195 observation wells (selected from sampling) were analyzed using PCA in IBM SPSS (version 22.0) software. The major output of the PCA results were two matrices i.e., "total variance explained" and "rotated loading factor" matrix. The total variance explained matrix provided statistics about the variance of the data explained by each component. In our case, the first two components of PCA explained more than 98% of the variation of the groundwater head in the area. Hence, only the first two components were selected for further analysis. The rotated loading factor matrix gave a score between 0 and 1 based on the importance of a given well in explaining the groundwater head variation in the study area. In our case, 135 of the 195 wells scored greater than or equal to 0.90 in the rotated loading factor matrix. It was found that these 135 wells could represent the study area with an average error of 0.318 m and standard deviation of 0.95 m ( Table 1). The number of observation wells was reduced from 195 to 135 (i.e., about a 31% reduction) with an insignificant loss of information. This reduction is less compared to that shown in previous literature, in which the reduction in randomly installed monitoring wells was more than 70% [3,23]. The smaller reduction in observation wells occurred because the network was spatially optimized using a hexagonal pattern before the application of PCA.

Potentiometric Head Classification
The wells that were prioritized in PCA were used to obtain a potentiometric surface of the study area. The average head values of these 135 prioritized wells over three and half years were considered as points for IDW interpolation. IDW interpolation results consisted of a minimum GW head value of 22.89 m and a maximum head value of 41.37 m. The potentiometric head surface was classified into 19 classes with a variation of 0.97 m in each class, and three wells in each class (first at the head, second at the middle, and third at the tail) were selected as proposed locations for the installation of

Potentiometric Head Classification
The wells that were prioritized in PCA were used to obtain a potentiometric surface of the study area. The average head values of these 135 prioritized wells over three and half years were considered as points for IDW interpolation. IDW interpolation results consisted of a minimum GW head value of 22.89 m and a maximum head value of 41.37 m. The potentiometric head surface was classified into 19 classes with a variation of 0.97 m in each class, and three wells in each class (first at the head, second at the middle, and third at the tail) were selected as proposed locations for the installation of observation wells. In addition, a few additional wells were also selected based on the interpolation results, such as, a few pockets of water which were part of some classes but were not continuous with that class. If there were already three or fewer wells in any class, they were left unaltered. A total of 59 wells were selected as the critical locations of the study area.
Potentiometric head surfaces of these 59 wells were generated using an IDW interpolation for different time series and compared against the base case scenario. The results showed that these 59 observation wells could represent the study area with an average error of 0.61 m and a standard deviation of 1.15 m ( Table 1). The groundwater head classification helped in reducing the number of observation wells from 135 to 59, which was more than a 56% reduction, with an insignificant loss of information. The overall performance of the network was satisfactory as it could represent the area with an average error of 0.61 m.

Discussion
The spatial distribution of error in the study area shows the effectiveness of the designed monitoring network. Figure 10a (195) shows the spatial distribution of error in the study area for post-monsoon 2010, when the area was monitored with 195 observation wells. It can be seen from the figure that 195 observation wells can represent the study area with a negligible degree of error, except at one location in the south of the model area where the groundwater head changes very rapidly within a small spatial distance, thus creating a nugget effect. The analysis proves the robustness of a hexagonal pattern of sampling in the case of continuous spatial functions. When the wells were reduced with PCA, a degree of error appeared in the east of the study area (Figure 10a (135)). The error in the east appeared because there was no variation in the groundwater head in that area. This area was waterlogged or had shallow water tables. Such conditions occur in most of the head reaches of the canal system, and this area is immediately adjacent to the main Rohri canal. PCA therefore does not consider this area important, so even with this stipulation, wells in the entire area did not score greater than or equal to 0.90 in the rotated factor loading matrix. The error also increased in the south of the model due to the removal of the wells in the area of very low groundwater head. Finally, when the wells were reduced to 59 through head classification, the error remained almost the same in the south of the model (Figure 10a (59)) but increased to some extent in the area with no flux boundary due to the exclusion of wells in that area. Furthermore, groundwater head contours were drawn in the study area by taking 3415, 195, 135, and 59 wells into account, respectively (Figure 10b). Results showed that the change in contours with the reduction of wells was insignificant. In general, contours were almost the same for all three designs and flow direction was conserved.
The other important factor to rationalize the proposed network was that it can identify the depleting water tables for different pumping regimes, which may exist in the future. This is done by simulating the system for longer periods and assessing whether the proposed network can represent the depletion of freshwater lenses. Figure 11 shows the depth to the water table (i.e., groundwater table from the top of the natural surface) for three pumping regimes in 2035. Figure 11a shows the regime with a pumping density of 5 wells/100 ha, Figure 11b shows the regime with a pumping density of 10 wells/100 ha, and Figure 11c shows the regime with a pumping density of 20 wells/100 ha. Distinct water table zones were visible for each of the pumping regimes. Waterlogged and shallow water table zones could be seen at head of the system, near the Rohri canal. Then there was a transition zone of moderate water tables at the middle of the system, and deep water tables near the Indus River. As the pumping increases in the region, the zone of the deep water table expands. An objective of the monitoring network was to capture the changes in the deep water tables, and ensure that propagation of this zone was represented completely by the monitoring network for different pumping regimes. The 195 well monitoring network described the distinct zones in all pumping regimes, while the 135 well monitoring network underestimated the waterlogging and shallow depth zone but represented the deep water zone. The 59 well network overestimated the deep water zone as the pumping increased in the south-east part of the study area. This was because all of the screening processes in the monitoring design were based on the removal of those wells that were under shallow water depths, and does not reflect temporal fluctuations. As per the focus of the monitoring network design for this study, the 59 well network will represent the depleting areas. If decision makers also need to monitor water logging areas, then the network should be expanded to 195 wells.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 20 boundary due to the exclusion of wells in that area. Furthermore, groundwater head contours were drawn in the study area by taking 3415, 195, 135, and 59 wells into account, respectively (Figure 10b). Results showed that the change in contours with the reduction of wells was insignificant. In general, contours were almost the same for all three designs and flow direction was conserved. The other important factor to rationalize the proposed network was that it can identify the depleting water tables for different pumping regimes, which may exist in the future. This is done by simulating the system for longer periods and assessing whether the proposed network can represent the depletion of freshwater lenses. Figure 11 shows the depth to the water table (i.e., groundwater table from the top of the natural surface) for three pumping regimes in 2035. Figure 11a shows the regime with a pumping density of 5 wells/100 ha, Figure 11b shows the regime with a pumping density of 10 wells/100 ha, and Figure 11c shows the regime with a pumping density of 20 wells/100 ha. Distinct water table zones were visible for each of the pumping regimes. Waterlogged and shallow water table zones could be seen at head of the system, near the Rohri canal. Then there was a transition zone of moderate water tables at the middle of the system, and deep water tables near the Indus River. As the pumping increases in the region, the zone of the deep water table expands. An objective of the monitoring network was to capture the changes in the deep water tables, and ensure that propagation of this zone was represented completely by the monitoring network for different pumping regimes. The 195 well monitoring network described the distinct zones in all pumping as the pumping increased in the south-east part of the study area. This was because all of the screening processes in the monitoring design were based on the removal of those wells that were under shallow water depths, and does not reflect temporal fluctuations. As per the focus of the monitoring network design for this study, the 59 well network will represent the depleting areas. If decision makers also need to monitor water logging areas, then the network should be expanded to 195 wells.

Conclusions
In Pakistan, there have been several shifts in groundwater management. Initially the focus was on managing the high water tables in the canal command areas, followed by a diversion in focus to encouraging private tube well development for increasing agricultural income [41]. Presently, the focus is on protecting aquifers from pollution and unsustainable pumping. National Water Policy (NWP) 2018 [42] emphasizes the establishment of monitoring networks to set sustainable groundwater yields and avoid the vertical movement of saline water. As per the NWP, provincial governments have to regulate and develop policy frameworks for sustainable groundwater management. The first step towards the formulation of effective policy is to generate knowledge about groundwater use and demarcate critical areas for groundwater management.
Of the freshwater lenses in the northern Rohri canal command, Sindh is a critical area for groundwater management as the livelihood of its community is highly dependent on the use of fresh groundwater. Its overexploitation will have huge impacts on the livelihood of the community. It is essential that this lens be monitored and regulated for sustaining the livelihood of the community. In this study, we designed monitoring networks for freshwater lenses using a four-step approach. In the first step, groundwater levels were simulated. In the second step, sampling was performed in a hexagonal pattern. In the third step, PCA was applied at temporal levels, and wells with low variance were removed; and in the last step, minimization of wells was performed for representation of potentiometric contours. The key outcomes, limitations, and future prospects from the study are: • A hexagonal sampling pattern resulted in a 195 well monitoring network, which represented the area with a mean error of 0.098 m compared to the base case. • PCA resulted in a 135 well monitoring network, which represented the area with a mean error 0.318 m.

•
Well reduction based on the contour classification resulted in a 59 well monitoring network, which represented the area with 0.61 m of error.

•
Three networks are presented for monitoring and management of the freshwater lenses, which can be established based on the available budget and monitoring targets. One with 195 wells will fully represent the water levels in the freshwater lenses. Designs with 135 and 59 wells will represent the depleting area, but they will not capture the water logging area. As the objective is to monitor the over-exploitation of the freshwater lenses, monitoring networks with 59 wells will help to achieve this objective.

•
The objective of this study was to model groundwater hydraulics and design a groundwater level monitoring network that can represent depleting and increasing water levels. As per the physical setting of the aquifer in the LIB, saline water up-coning will take place if the thickness of shallow freshwater lenses becomes thin following over-extraction [27,28,43,44]. It is essential to set a water level threshold that should be kept in freshwater lenses to avoid saline water intrusion.
In future studies, we will extend the flow model to perform freshwater-saline water interaction to set a groundwater level threshold. Water level measurements in combination with thresholds will provide regulators with a decision tool to regulate pumping in freshwater lenses by only monitoring levels in the piezometers.