Assessment of Groundwater Flow Dynamics Using MODFLOW in Shallow Aquifer System of Mahanadi Delta (East Coast), India

: Despite being a biodiversity hotspot, the Mahanadi delta is facing groundwater salinization as one of the main environmental threats in the recent past. Hence, this study attempts to understand the dynamics of groundwater and its sustainable management options through numerical simulation in the Jagatsinghpur deltaic region. The result shows that groundwater in the study area is extensively abstracted for agricultural activities, which also causes the depletion of groundwater levels. The hydraulic head value varies from 0.7 to 15 m above mean sea level (MSL) with an average head of 6 m in this low-lying coastal region. The horizontal hydraulic conductivity and the speciﬁc yield values in the area are found to vary from 40 to 45 m/day and 0.05 to 0.07, respectively. The study area has been calibrated for two years (2004–2005) by using these parameters, followed by the validation of four years (2006–2009). The calibrated numerical model is used to evaluate the net recharge and groundwater balance in this study area. The interaction between the river and coastal unconﬁned aquifer system responds differently in different seasons. The net groundwater recharge to the coastal aquifer has been estimated and varies from 247.89 to 262.63 million cubic meters (MCM) in the year 2006–2007. The model further indicates a net outﬂow of 8.92–9.64 MCM of groundwater into the Bay of Bengal. Further, the outﬂow to the sea is preventing the seawater ingress into the shallow coastal aquifer system.


Introduction
The coastal aquifer is one of the most important water resources in coastal regions that supplies water to more than a billion people worldwide [1][2][3] and connects the world's oceanic and hydrologic systems [4]. The general hydrogeological characteristics of the coastal aquifers are influenced by geologic environments, mixing zones, long and shortterm sea fluctuations, and the density gradients due to differences in salinity [1]. In general, groundwater is an attractive source of water (15,300 × 10 3 km 3 ) as it is fresh and readily available [5,6]. However, the groundwater of coastal regions is more susceptible to deterioration due to several factors, such as rapid urbanization, intensified agricultural development, economic development, climate change, sea-level rise, and lack of sufficient surface water resources [7][8][9]. The effect of groundwater withdrawal on coastal aquifer systems with the smallest topographic gradient is more pronounced than the impact of sea-level rise and variation in groundwater recharge due to the dynamic nature of

Hydrogeology
The deltaic region of Jagatsinghpur belongs to the quaternary formation that covers recent sediments of flood plain deposits of the Mahanadi River and the Devi River. These are mainly comprised of gravel, sand (fine to coarse grain), silt, and clay (black, red, and yellow) materials [42]. These unconsolidated to semi-consolidated materials act as a good repository of groundwater resources. The lower part of this study area lying close to the coast is characterized by low lying wet plains, fine-grained sediments, tidal infected rivers, tidal creeks, swamps, ill drainage of land, and non-development of levees [37]. This coastal tract acts as a favorable zone of groundwater availability due to the large thickness of sediments of varying sizes deposited in this part of the study area. The aquifer system of the Jagatsinghpur area is divided mainly into two zones, i.e., a shallow aquifer zone (<50 m thickness) and deeper aquifer zone (50-300 m thickness) below ground level [35].

Methodology
The groundwater flow simulation can be performed through Visual MODFLOW, which integrates the modular 3D finite-difference groundwater flow code [20]. Several numerical codes are used to simulate groundwater flow both in local and regional groundwater systems [43]. The groundwater flow equation (Equation (1)) is used in MODFLOW to know the groundwater flow in three different directions for this study [12].
where Kx, Ky, and Kz refer to hydraulic conductivity in three different directions; Ss, h, and R represent specific storage, hydraulic head, and sink or source, respectively. Visual MODFLOW is based on the finite-difference mathematical equation (Equation (2)) with assumptions of constant density and viscosity of groundwater flow under transient state conditions [44].
Hydraulic head value does not change with time under steady-state conditions. This condition expresses as (Equation (3)).

Hydrogeology
The deltaic region of Jagatsinghpur belongs to the quaternary formation that covers recent sediments of flood plain deposits of the Mahanadi River and the Devi River. These are mainly comprised of gravel, sand (fine to coarse grain), silt, and clay (black, red, and yellow) materials [42]. These unconsolidated to semi-consolidated materials act as a good repository of groundwater resources. The lower part of this study area lying close to the coast is characterized by low lying wet plains, fine-grained sediments, tidal infected rivers, tidal creeks, swamps, ill drainage of land, and non-development of levees [37]. This coastal tract acts as a favorable zone of groundwater availability due to the large thickness of sediments of varying sizes deposited in this part of the study area. The aquifer system of the Jagatsinghpur area is divided mainly into two zones, i.e., a shallow aquifer zone (<50 m thickness) and deeper aquifer zone (50-300 m thickness) below ground level [35].

Methodology
The groundwater flow simulation can be performed through Visual MODFLOW, which integrates the modular 3D finite-difference groundwater flow code [20]. Several numerical codes are used to simulate groundwater flow both in local and regional groundwater systems [43]. The groundwater flow equation (Equation (1)) is used in MODFLOW to know the groundwater flow in three different directions for this study [12].
where K x , K y , and K z refer to hydraulic conductivity in three different directions; S s , h, and R represent specific storage, hydraulic head, and sink or source, respectively. Visual MODFLOW is based on the finite-difference mathematical equation (Equation (2)) with assumptions of constant density and viscosity of groundwater flow under transient state conditions [44].
Hydraulic head value does not change with time under steady-state conditions. This condition expresses as (Equation (3)). Visual MODFLOW has several solvers, such as Preconditioned Conjugate Gradient (PCG), Strongly Implicit Procedure (SIP) package, WHS solver for Visual MODFLOW package (WHS), Slice Successive Over Relaxation (SOR) package, and Geometric Multigrid solver (GMG) package, are used to solve the numerical equation for groundwater flow simulation purposes. For this study, the WHS solver package with Bi-Conjugate Gradient Stabilized (Bi-CGSTAB) accelerator us used to resolve the partial differential equations through iterative procedures.

Development of the Model
The groundwater model starts with the development of a groundwater flow model for a particular study area, which represents its physical condition. Similarly, for this study, a model consisting of a single layer has been conceptualized based on geology, lithologs, river boundary conditions, and groundwater level data sets [45]. The shallow unconfined aquifer composed of unconsolidated formation up to 50 m depth has been considered as the modeled single layer. After preparing a conceptual model, it is translated into a numerical model with the help of the Visual MODFLOW software package. The numerical model is developed through several steps. The input parameters in the conceptual model are presented in Table 1.

Hydraulic Head Data
Hydraulic head data of eleven (11) different observation wells were collected from the Central Ground Water Board (CGWB) in the Jagatsinghpur coastal aquifer system for this study. The hydraulic head varies from 0.7 m near the sea coast to 15 m away from the shoreline. For groundwater flow simulation, 1st January 2004 has been taken as the initial time. Annually, four different periods of head data (from the year 2004 to 2009) have been used for both calibration and validation of the model. These head data have been categorized as post-monsoon (Rabi), pre-monsoon, monsoon, and post-monsoon (Kharif).

Boundary Conditions
The Visual MODFLOW simulates the groundwater flow followed by different types of boundary conditions. In the Jagatsinghpur coastal aquifer system, two types of boundaries have been used ( Figure 2). The Bay of Bengal is considered as a constant head boundary or Dirichlet boundary [46]. Two large perennial rivers, i.e., the Mahanadi River and its distributary the Devi River, flowing along the two flanks of the study area are considered as the Cauchy boundary or head-dependent flux boundary. The river conductance value can be determined from Equation (4) [29].
where Kr = hydraulic conductivity of the river bed (m/day), L = length of the reach/grid size (m) Wr = width of the river (m), and B = thickness of the river bed (m). The river bed conductance of two rivers is approximately the same, i.e., 30,000-35,000 m 2 /day [35].

Hydrological Parameters
The model domain is classified into three hydraulic conductivities and specific yield zones for this single unconfined aquifer system (Figure 2), which also belongs to the

Hydraulic Head Data
Hydraulic head data of eleven (11) different observation wells were collected from the Central Ground Water Board (CGWB) in the Jagatsinghpur coastal aquifer system for this study. The hydraulic head varies from 0.7 m near the sea coast to 15 m away from the shoreline. For groundwater flow simulation, 1 January 2004 has been taken as the initial time. Annually, four different periods of head data (from the year 2004 to 2009) have been used for both calibration and validation of the model. These head data have been categorized as post-monsoon (Rabi), pre-monsoon, monsoon, and post-monsoon (Kharif).

Boundary Conditions
The Visual MODFLOW simulates the groundwater flow followed by different types of boundary conditions. In the Jagatsinghpur coastal aquifer system, two types of boundaries have been used ( Figure 2). The Bay of Bengal is considered as a constant head boundary or Dirichlet boundary [46]. Two large perennial rivers, i.e., the Mahanadi River and its distributary the Devi River, flowing along the two flanks of the study area are considered as the Cauchy boundary or head-dependent flux boundary. The river conductance value can be determined from Equation (4) [29].
where K r = hydraulic conductivity of the river bed (m/day), L = length of the reach/grid size (m) W r = width of the river (m), and B = thickness of the river bed (m). The river bed conductance of two rivers is approximately the same, i.e., 30,000-35,000 m 2 /day [35].

Hydrological Parameters
The model domain is classified into three hydraulic conductivities and specific yield zones for this single unconfined aquifer system (Figure 2), which also belongs to the alluvial formation of the Mahanadi delta. The horizontal hydraulic conductivities (K h ) are 40 m/day, 42 m/day, and 45 m/day for ZONE I, ZONE II, and ZONE III, respectively, whereas the corresponding vertical conductivity (K v ) of the three zones is 4, 4.2, and 4.5 m/day. The different specific yield values of 0.05, 0.06, and 0.07 for three respective zones I-III were taken during the calibration of the groundwater model ( Table 2). As the water table is very close to the ground surface, some groundwater is extracted through the evapotranspiration process. Hence the evapotranspiration data have been taken into consideration for the groundwater simulation model. Further, the study area has been divided into eleven different recharge zones, in which monthly rainfall recharge values have been assigned.

Calibration and Validation of Model
Calibration is the process through which the calculated head value is subjected to match with the observed head value. The head values from the year 2004 to the year 2005 are taken for calibrating the model for this study. In the present study, the calibration of the model is done through trial and error in which the unknown hydrogeologic parameters are set to be fixed to minimize the head difference between the calculated head and observed head at steady-state conditions. Then, the model is run for 2 years from January 2004 to December 2005 under transient state conditions. The groundwater simulation is said to be a good fit when the computed head value is very close to the observed head value and this good match can be analyzed through calibration criteria, e.g., the mean error (ME) (Equation (5)), the mean absolute error (MAE) (Equation (6)), and the root mean squared error (RMSE) (Equation (7)) [12,47,48]. After calibration, the validation of the model is performed by taking the hydraulic head values from January 2006 to December 2009.

Mean Error
where h o refers to the observed head value, h c the calculated head value, and n the total number of observed data. A statistical analysis of calibrated model under steady-state conditions has been given ( Figure 3). Under transient state conditions, the calibrated and validated model indicates a good correlation between the observed head and calculated head in the study area (Figure 4a

Parameter Estimation (PEST) Model
The groundwater model parameters are also estimated by using PEST [49]. Knowling and Adrian (2016) used PEST to minimalize the weighted least squares objective function based on Tikhonov regularization [50]. In this study, the groundwater model has been auto-calibrated with the help of the PEST module of MODFLOW to optimize the aquifer

Parameter Estimation (PEST) Model
The groundwater model parameters are also estimated by using PEST [49]. Knowling and Adrian (2016) used PEST to minimalize the weighted least squares objective function based on Tikhonov regularization [50]. In this study, the groundwater model has been auto-calibrated with the help of the PEST module of MODFLOW to optimize the aquifer

Parameter Estimation (PEST) Model
The groundwater model parameters are also estimated by using PEST [49]. Knowling and Adrian (2016) used PEST to minimalize the weighted least squares objective function based on Tikhonov regularization [50]. In this study, the groundwater model has been auto-calibrated with the help of the PEST module of MODFLOW to optimize the aquifer parameters (hydraulic conductivity and specific yield). The calibrated model shows a good correlation between the observed head value and calculated head value with a correlation coefficient value of 0.993 ( Figure 5). The automated calibrated (PEST) aquifer parameters (conductivity and specific yield) have been compared to the manually calibrated aquifer parameters, as shown in Table 3. The uncertainty is the process through which the uncertainty on the estimated parameters is quantified to understand the risk associated with different groundwater management models [51].
Water 2022, 14, x FOR PEER REVIEW 9 parameters (hydraulic conductivity and specific yield). The calibrated model sho good correlation between the observed head value and calculated head value w correlation coefficient value of 0.993 ( Figure 5). The automated calibrated (PEST) aq parameters (conductivity and specific yield) have been compared to the man calibrated aquifer parameters, as shown in Table 3. The uncertainty is the process thro which the uncertainty on the estimated parameters is quantified to understand the associated with different groundwater management models [51].  The estimated parameters viz. specific yield and hydraulic conductivity by the P tool have been subjected to uncertainty analysis, which shows a 95% confidence int between 40.19 and 44.96 m/day (Table 4).  The estimated parameters viz. specific yield and hydraulic conductivity by the PEST tool have been subjected to uncertainty analysis, which shows a 95% confidence interval between 40.19 and 44.96 m/day (Table 4).

Interaction between Aquifer and River
Both inflow from the rivers into the aquifer system and outflow from the aquifer system to the rivers have been observed in a different time period. This unconfined coastal aquifer receives water from rivers in the pre-monsoon and post-monsoon period, whereas the excess amount of groundwater in the form of base flow discharges to the rivers during monsoon season, as shown in Figure 6a,b. The inflow from the river boundary has been estimated as 34 MCM in the post-monsoon and pre-monsoon period in the year 2006-2007. In the monsoon period, the unconfined coastal aquifer system supplies around 23 to 27 MCM of groundwater to the river system after irrigation (Figure 6a,b). This deltaic aquifer system is mostly recharged by rainfall during wet days. Figure 6c shows that the extraction of groundwater is different in different time periods. In the post-monsoon time period, withdrawal of groundwater is more than that of the pre-monsoon period to provide water for post-monsoon crop, though there is available of adequate amount of water in coastal areas to meet the monsoon period Kharif crop [52]. The extracted groundwater for the sustainability of agricultural productivity and livelihoods in the post-monsoon season has been estimated at around 180 MCM in the year 2006-2007 (Figure 6c). This implies that the heavy abstraction of groundwater for agriculture activity and other domestic uses declines the groundwater level of the aquifer system, which is also a respondent of river inflow into the aquifer system. The resultant river inflow of 33.92 MCM of water entering into this coastal aquifer is due to the pumping of groundwater. Similarly, during the pre-monsoon time, the groundwater is extracted to fulfill the water demand for Rabi crops, but the amount of water required is less than that of the post-monsoon season. It is calculated that about 100 MCM of groundwater is pumped out for agricultural activity and other needs, which also causes the inflow of water through the river boundary. When there is groundwater stress, whether less or more, it also affects the river system.
The river-aquifer interaction indicates a good relationship between the river stage and the outflow/inflow from the river boundary (Figure 7a,b).
As the outflow from the river boundary increases, the river stage also increases. The estimated base flow (7.81 MCM) to the river could be one of the factors contributing to the highest river stage value of 8 m in August of the monsoon period. In contrast, the inflow from river boundary to aquifer system during pre-and post-monsoon time causes declination of river stage from 8 to 3 m. According to the estimation, about 10 MCM of water from the river system enter into the aquifer systems during this season. The interpreted result shows the influence of groundwater flux on the river stage and also suggests a good interaction between the river and the coastal aquifer system [52].

Fluctuation in Groundwater Level
The spatiotemporal variation in groundwater level has been identified in this coastal aquifer system (Figure 8).
during the pre-monsoon time, the groundwater is extracted to fulfill the water demand for Rabi crops, but the amount of water required is less than that of the post-monsoon season. It is calculated that about 100 MCM of groundwater is pumped out for agricultural activity and other needs, which also causes the inflow of water through the river boundary. When there is groundwater stress, whether less or more, it also affects the river system.  The river-aquifer interaction indicates a good relationship between the river stage and the outflow/inflow from the river boundary (Figure 7a,b). As the outflow from the river boundary increases, the river stage also increases. The estimated base flow (7.81 MCM) to the river could be one of the factors contributing to the highest river stage value of 8 m in August of the monsoon period. In contrast, the inflow from river boundary to aquifer system during pre-and post-monsoon time causes declination of river stage from 8 to 3 m. According to the estimation, about 10 MCM of water from the river system enter into the aquifer systems during this season. The interpreted result shows the influence of groundwater flux on the river stage and also suggests a good interaction between the river and the coastal aquifer system [52].

Groundwater Recharge Estimation
In this area, rainfall, seepage from the riverbed, and irrigation return flow are the major sources of groundwater recharge. Groundwater recharge estimation plays a vital role in the optimal development and efficient management of fresh groundwater resources in coastal areas. The study area of the Jagatsinghpur district is divided into eleven recharge zones (Figure 9) in the form of Thiesen polygons to estimate the groundwater recharge [43]. The observation well (OW 6) is showing the highest hydraulic head rise of 3.07 m, situated away from the coastal tract of the study area. The lowest hydraulic head rise of 0.75 m has been observed in the observation well (OW 10), which is close to the Bay of Bengal. This spatial variation in hydraulic head rise depends on the topography, rainfall intensity, type of soil, and land-use patterns. Groundwater stress has been observed during the pre-monsoon period. This also results in a decline in groundwater levels due to the absence of rainfall events in between days 1 and 89, as shown in Figure 8. In the case of monsoon and post-monsoon periods (243 and 334 days), the hydraulic head increases from place to place, suggesting that the groundwater is being mostly recharged by rainfall and regained the groundwater potentiality in the study area. The contour lines in the upper part of the study area are very close, which reveals the presence of a high hydraulic gradient and encounters high groundwater movement [53]. Gradually, large spacing between two equipotential lines has been observed close to the sea shoreline, indicating the sluggish movement of groundwater as the presence of a low hydraulic gradient (Figure 8). However, there is an average rise of 1.84 m of hydraulic head in the monsoon period as compared to the pre-monsoon period, which implies that the shallow aquifer of this coastal region is recharged quickly due to rainfall events. This suggests that sustainable management of coastal aquifers is required during dry periods as there is an absence of a primary recharge source (i.e., rainfall).

Groundwater Recharge Estimation
In this area, rainfall, seepage from the riverbed, and irrigation return flow are the major sources of groundwater recharge. Groundwater recharge estimation plays a vital role in the optimal development and efficient management of fresh groundwater resources in coastal areas. The study area of the Jagatsinghpur district is divided into eleven recharge zones (Figure 9) in the form of Thiesen polygons to estimate the groundwater recharge [43].

PEER REVIEW
13 of 17 Figure 9. Recharge zones in the study area.
In the groundwater simulation model, the recharge was manually optimized to minimize the head difference between the observed head and calculated head [46]. Modeled recharge rates vary across the recharge zones and also from year to year due to large variations of rainfall over the simulation period. The temporal and spatial distribution of annual rainfall from the year 2004 to 2009 varies from 875 to 1229 mm. The average net recharge (rainfall recharge-groundwater draft) in the area varies from 247.89 to 262.63 million cubic meters (MCM) in the year 2006-2007. The net recharge in post-monsoon is estimated to be less than that of pre-monsoon and monsoon periods ( Figure 10). Around 6-15 MCM of water recharge the aquifer system during the post-monsoon period, whereas 33-54 MCM of water percolates into the aquifer system in the pre-monsoon period. Huge extraction of groundwater leads to less net recharge in the post-monsoon season, though natural rainfall happens to have occurred in the study area. As compared to the post and pre-monsoon season, the coastal aquifer system gets highly recharged by rainfall in the monsoon period, estimated in between 180.26 and 223.08 MCM. In the groundwater simulation model, the recharge was manually optimized to minimize the head difference between the observed head and calculated head [46]. Modeled recharge rates vary across the recharge zones and also from year to year due to large variations of rainfall over the simulation period. The temporal and spatial distribution of annual rainfall from the year 2004 to 2009 varies from 875 to 1229 mm. The average net recharge (rainfall recharge-groundwater draft) in the area varies from 247.89 to 262.63 million cubic meters (MCM) in the year 2006-2007. The net recharge in post-monsoon is estimated to be less than that of pre-monsoon and monsoon periods ( Figure 10). Around 6-15 MCM of water recharge the aquifer system during the post-monsoon period, whereas 33-54 MCM of water percolates into the aquifer system in the pre-monsoon period. Huge extraction of groundwater leads to less net recharge in the post-monsoon season, though natural rainfall happens to have occurred in the study area. As compared to the post and pre-monsoon season, the coastal aquifer system gets highly recharged by rainfall in the monsoon period, estimated in between 180. 26

Groundwater Outflow to the Bay of Bengal
The groundwater flow direction in the study area is towards the Bay of Bengal as shown in Figure 8. The model simulation results show that the outflow to the Bay of Bengal varies from 8.92 to 9.64 MCM on an annual basis (2006-2007) ( Figure 11). The unconfined coastal aquifer of the studied area discharges nearly 50% of its groundwater during the monsoon period, while the rest discharges during dry periods [29]. This is due to the groundwater recharge by rainfall and reduction of groundwater abstraction during the monsoon season. The resultant outflow from the Bay of Bengal also prevents seawater ingress into the aquifer system.

Conclusions
The groundwater dynamics in the coastal area of the Mahanadi delta was studied using a modeling technique with the help of Visual MODFLOW. The simulated heads

Groundwater Outflow to the Bay of Bengal
The groundwater flow direction in the study area is towards the Bay of Bengal as shown in Figure 8. The model simulation results show that the outflow to the Bay of Bengal varies from 8.92 to 9.64 MCM on an annual basis (2006-2007) ( Figure 11). The unconfined coastal aquifer of the studied area discharges nearly 50% of its groundwater during the monsoon period, while the rest discharges during dry periods [29]. This is due to the groundwater recharge by rainfall and reduction of groundwater abstraction during the monsoon season. The resultant outflow from the Bay of Bengal also prevents seawater ingress into the aquifer system.

Groundwater Outflow to the Bay of Bengal
The groundwater flow direction in the study area is towards the Bay of Bengal as shown in Figure 8. The model simulation results show that the outflow to the Bay of Bengal varies from 8.92 to 9.64 MCM on an annual basis (2006-2007) ( Figure 11). The unconfined coastal aquifer of the studied area discharges nearly 50% of its groundwater during the monsoon period, while the rest discharges during dry periods [29]. This is due to the groundwater recharge by rainfall and reduction of groundwater abstraction during the monsoon season. The resultant outflow from the Bay of Bengal also prevents seawater ingress into the aquifer system.

Conclusions
The groundwater dynamics in the coastal area of the Mahanadi delta was studied using a modeling technique with the help of Visual MODFLOW. The simulated heads

Conclusions
The groundwater dynamics in the coastal area of the Mahanadi delta was studied using a modeling technique with the help of Visual MODFLOW. The simulated heads matched significantly with the observed heads characterized by different statistical parameters, showing the development of a robust model for this study area. The aquifer parameters (hydraulic conductivity and specific yield) estimated by the PEST tool and trial-and-error method through modeling are approximately the same. Further, the parameters have been subjected to uncertainty analysis, yielding a 95% confidence interval between 39.70 and 48.43 m/day for a particular zone.
The shallow coastal aquifer was influenced by the river system. The estimated groundwater abstraction (180 MCM) led to inflow from the river during the non-monsoon time. Further, the excess amount of groundwater as base flow during the monsoon period recharges the river. This is one of the controlling factors that also affects the river stage. The good connectivity between river and shallow aquifer indicates the regular water circulation, exchange from groundwater to surface water or vice versa. This can improve the aquatic environmental condition. On the other hand, the decline of groundwater was due to the heavy extraction of groundwater to grow non-monsoon agricultural products. Furthermore, the temporal variation of the hydraulic head reveals that the shallow coastal aquifer is very sensitive to rainfall as it quickly responds to rainfall events.
Based on the estimated aquifer parameters, i.e., hydraulic conductivity and specific yield, the net groundwater recharge to this coastal aquifer was estimated to vary between 247.89 and 262.63 MCM. The results derived from the groundwater modeling indicate that there is a net outflow of groundwater into the Bay of Bengal. The outflow varies between 8.92 and 9.64 MCM, which may prevent the seawater ingress into the coastal aquifer of Jagatsinghpur, Odisha. Further, this scientific evidence may prove to be significant for the development of resilient water development plans for dynamic coastal aquifers, so that future generations can utilize the precious resource against climate-related hazards.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.