Assessment of contamination dispersion in a porous aquifer environment using explicit and implicit methods in the MODFLOW–MODPATH model

Abstract

This study presents a meticulously edited article that assesses the uncertainty associated with contamination dispersion in a porous aquifer environment. The article employs a conceptual model incorporating finite difference flow modeling and particle tracking methods grounded in particle movement theory within groundwater flow. The groundwater flow model utilizes an automatic calibration technique, PES, and independent verification, to minimize statistical discrepancies in the optimized parameters. The statistical summary of the model output reveals the average pollution concentration in the study area, specifically the landfill aquifer in northern Iran. The results indicate that the initial contamination value at the boundaries is zero, whereas 91,802 units are discharged into the sea. Moreover, approximately 3,317,290,003 pollution units have been extracted from pumping wells along the contamination plume over 10 years. Furthermore, the drainage network releases 1,778,680,001 pollution units, resulting in a net outflow of 5,082,740,007 units from the same drainage network. Additionally, 9,912,180,003 units were discharged from the leaky boundaries. It should be noted that these statistics can be reversed depending on the specified period. The difference between the input and output pollution values in the study area indicates a net accumulation of 8192 units. To analyze the sensitivity of contamination dispersion in the porous aquifer environment and the variations in coastal drainage patterns obtained from the finite difference simulation, the explicit and implicit methods, as well as the random walk particle tracking method, were investigated in a portion of the basin area using the South Side Implicit method. One of the main objectives of this sensitivity analysis is to evaluate the results of the uncalibrated qualitative model and examine the particle dispersion in the adjacent area of the contaminated landfill region to assess the error of the proposed model by developing a limited test. The calculated error or computational difference between the explicit and implicit methods and the applet model is considered negligible. For example, the calculated error between the explicit and backward methods and the applet model ranges from 0.086 to 0.372 units for a concentration interval of 0.512 units and 0.680 units for a time interval 0.537. Therefore, the computational differences are not significant. Furthermore, based on theoretical sensitivity analysis of contamination dispersion, the laboratory model demonstrates that the flow direction in this modeling is from west to east, as expected.

Introduction

The decline in groundwater levels in plains has significant effects, leading to increased costs of water extraction and increased consumption of nitrate fertilizers to enhance productivity, consequently reducing water quality. However, uncertainty in assessing pollution release in porous environments remains a significant factor of distrust in particle dispersion models used for their exploitation. Nonetheless, scientific sources use these models to develop groundwater simulations as a decision-making tool. Each factor contributing to groundwater pollution also has secondary negative consequences.

For example, the reduction in the extractable water volume of groundwater sources and increased costs of water extraction result in a decrease in cultivated area and an increase in product prices. Another factor contributing to the critical state of groundwater resources is the alteration of natural recharge conditions and groundwater flow regime. In some plains, improper hydraulic constructions on rivers have reduced the appropriate recharge of downstream aquifers. Some dams, irrigation and drainage projects, and river regulation projects have undesirable effects on the natural recharge of groundwater reservoirs (Ghodrati 2012).

The concentration of population in a specific area, despite having reliable water, will increase due to the decline in water resources and the difficulty in accessing it, leading to the destruction of a significant portion of the capacity for groundwater flow and replenishment. Changes like groundwater flow will fundamentally alter the water situation (Ayedi et al. 2018). The ultimate decline in water quality will undoubtedly lead to environmental crises, making it one of the most significant impacts of this phenomenon.

Understanding the actual behavior of a natural system necessitates long-term research for each specific region. Considering the current situation and the existing constraints, using numerical models as valuable and efficient tools have been suggested and have yielded satisfactory results. In recent years, numerous projects centered around groundwater utilization, protection, and purification have primarily focused on developing groundwater models. The role of these models has also grown with the continuous advancements in computer hardware and software availability (Chien et al. 2018).

Managing groundwater in a particular location demands the maximum possible supply of groundwater to meet the needs of all users within the study area while considering predetermined conditions such as aquifer conditions, restricting development and exploitation, and addressing specific legal, social, political, and economic limitations. The management approach should be designed to prevent excessive pumping’s actual and potential adverse effects on water quality deterioration and overexploitation (Gaji et al. 2020).

Consequently, understanding the groundwater status can assist policymakers and decision-makers in protecting and managing limited water resources in the long term (Main et al. 2019). Planning for the sustainable extraction and utilization of groundwater requires a realistic and accurate assessment of groundwater depletion’s impact on water tables and associated natural ecosystems (Gaji et al. 2020). Evaluations of changes aim to preserve water for ecosystems as defined by environmental requirements; however, sustainable groundwater security also necessitates maintaining sustainable groundwater quality for consumers and recipients, such as streams and groundwater-dependent ecosystems (Gaji et al. 2020).

The speed and magnitude of changes in the conditions governing groundwater resources highlight the need for computer models, especially in assessing water quality affected by various factors, such as recharge events that replenish aquifers or human-induced recharge from wells (Roy et al. 2020). The degradation of groundwater quality is closely tied to the extent and type of groundwater level decline resulting from human activities. Diminished recharge capacity and the obstruction of saturated flow extension are fundamental factors in reducing aquifer vulnerability (Arzami et al. 2019).

In cases where numerous uncertainties exist in studying the abrupt changes of saturated aquifer structure, a numerical method can calculate the system’s capacity reliability by analyzing sensitive and anomalous regions. This is because capacity reliability is related to inherent hydraulic inadequacy (Azimi et al. 2019a, b).

Phenomena or events that cannot be accurately predicted, such as the decline in the water table, are considered risky or uncertain. The term “risk” is often used to describe situations where having available probabilities is necessary to define the likelihood of various phenomena or consequences accurately. Uncertainty in information is inherent in future planning activities (Kashtgar et al. 2020). This uncertainty arises from insufficient information, incorrect assumptions, and the variability of natural processes over time and space. Water managers often need to determine the uncertainty and variations in system performance indicators due to changes in input data and predicted parameter values. They must assess this level of uncertainty for the practical development of engineering projects (Zhang et al. 2020).

In the plain’s region, due to localized declines caused by human activities, it is essential to implement prohibition measures to prevent the destruction of its recharge capacity. However, providing a decision-making model with minimal uncertainty requires using accurate numerical tools to reduce uncertainty to an acceptable level. Flexible screening tools capable of simultaneously analyzing large and diverse regional data will identify spatial or temporal changes in water quality, making this approach increasingly important for groundwater scientists today (Chien et al. 2018).

Understanding the seasonal and island-like aspects of natural systems, along with modeling and significant parameter deviations, can lead to destructive managerial decisions that consider the conflicting interests of stakeholders and decision-makers. That is why the precise determination and documentation of vulnerabilities significantly reduce the cost of achieving a sustainable ecosystem (Mine et al. 2020).

In a local decision-making system, determining final conditions requires well-documented studies. Documenting a specific region to define a vulnerability map of a threatened plain due to potential quality degradation events can only be accomplished using spatial maps. A practical solution is to use maps with quantitative figures to determine the likelihood of reducing secondary risks for each management decision. Therefore, using groundwater models as a cost-effective and efficient method in studying the movement, analyzing the budget, and managing the operation of groundwater has gained considerable importance.

The speed and magnitude of changes in conditions governing groundwater resources highlight the need for computer models, especially in assessing water quality impacted by various factors, such as recharge events that replenish aquifers or human-induced recharge from wells (Rai et al. 2020). The degradation of groundwater quality is closely tied to the extent and type of groundwater level decline resulting from human activities. Diminished recharge capacity and the obstruction of saturated flow extension are fundamental factors in reducing aquifer vulnerability (Arzami et al. 2019).

Depending on the assumptions regarding the desired phenomenon, various models, such as black box, conceptual, and mathematical models (analytical and numerical), exist as system representations. Among these models, conceptual models can be considered closer to reality than black box models. These models are prepared based on the existing conditions in the aquifer, such as geological conditions, flow types (sheet-like or turbulent), input and output flow components, boundary conditions, and more (Water and Wastewater Engineering Organization 2011).

From the perspective of water quality parameters, the spatial domain’s density dependence emphasizes the need for integrated and continuous analysis that considers this factor. Although these models significantly contribute to the initial awareness and understanding of the prevailing conditions in the aquifer, they cannot provide numerical results due to their descriptive nature. Therefore, they are not utilized in aquifer management. Consequently, mathematical models are employed to understand the aquifer better.

However, the complex structure of aquifers makes simple multivariable models less capable of generating predictive outputs for future intervals. Hence, numerical codes executed in a modified network of finite difference cells can more effectively represent various decision-making conditions on aquifer capacity in prediction scenarios (Datagupta et al. 2020). The results of this research can contribute to a more accurate prediction of future aquifer conditions under probable scenarios and facilitate planning and management with minimal adverse effects on water quality. Regional managers can utilize it as a decision-making mechanism for sustainable performance with lower uncertainty.

Materials and methods

Study area

The selected area encompasses the second-order basins between the Sefidroud and Haraz Rivers, the Haraz River, and the basins between Haraz and Gorganrood Rivers (Fig. 1). It extends across the Central Alborz Mountains’ northern slopes and the Mazandaran Sea’s southern coast, from the Sefidroud Delta to Bandar Gaz. The geographical coordinates of the study area range from 49°45′40″ E to 54°44′22″ E longitude and from 35°44′04″ N to 37°23′38″ N latitude. The total area of these regions is approximately 29,668 square kilometers, divided into 7 study zones and 78 sub-basins. This basin is considered part of the Mazandaran Sea watershed in the overall hydrological division of Iran. The Mazandaran Sea borders it to the north, the Sefidroud basin to the west, the Central Desert and the Salt Lake watershed to the south, and the Gorganrood basin to the east (Fig. 2).

Fig. 1

Geographic location of the study area—Babol landfill

Fig. 2

Schematic structure of the MODFLOW model for Babol plain

This study presents a carefully developed conceptual model of a limited difference network for groundwater flow. The computational framework employed for executing the groundwater flow model is MODFLOW (2005), which ensures stability and progressive analysis. The research acknowledges the prevailing uncertainty of specific input parameters exhibiting high errors. The PCGN computational engine was chosen as the preferred option for both stable and unstable conditions to address this. The computational process involved 100 iterations of Outer and Inner loops, with a convergence change of the critical limit set at 1 m and a convergence error critical limit of 1 m³/day. The subsequent execution of the model included calibration and validation over 10 years in a stable form. The entire process is visually illustrated in Fig. 3, showcasing the homogeneous form of sensitivity analysis conducted for the examined analysis parameters using the PEST code.

Fig. 3

Comparative sensitivity analysis of all parameters in the varying calibration stage

This study utilized a combined method of sensitivity analysis and calibration stages. The sensitivity analysis results emphasize the significant influence of hydraulic conductivity parameters and hydraulic conductivity heterogeneity, especially in specific pilot areas. Moreover, a particular group of rivers within the basin has been identified as highly sensitive. Based on this, adjustments were implemented to the boundaries that encounter high hydraulic load within the conceptual model. These modifications significantly reduced the overall calibration error during the final stage of the four-step calibration process, achieving a desirable minimum. The corresponding error is indicated in Fig. 4.

Fig. 4

Reduction of error in the final stage of unstable calibration

The RMS (Root-Mean-Square) or RMSE (Root-Mean-Square Error), as shown in the figure, equals 2.59, indicating high simulation accuracy. It needs to be standardized to understand this relative error’s significance better. In order to determine whether the flow is being transferred into the aquifer from the landfill site or not, the MODPATH model was utilized. However, due to the relatively short computational period in the current model, predicting the advancement of contamination into the aquifer requires more detailed investigations, especially considering that a decrease in water level practically increases the likelihood of infiltration.

We employed a particle transport model to investigate the dispersion of pollutants into the aquifer. Our primary objective was to analyze the contamination of the adjacent river caused by wastewater disposal. We specifically examined the re-entry of this contamination into isochronous regions where the flow converges. With the objectives of conducting a sensitivity analysis of pollution dispersion in the porous aquifer environment and analyzing changes in the drainage system of the border region, as indicated by the diagram obtained from the finite difference simulation model, we investigated the Explicit and Implicit methods, along with the random walk particle transfer method, for a selected portion of the plain area.

One of the main goals was to examine the sensitivity of the uncalibrated qualitative model results and assess the neighboring environment of the contaminated area (Landfill) to particle dispersion for evaluating the accuracy of the proposed model through developing a limited test. Therefore, the dispersion model diagrams were initially obtained at time intervals of 0 and 10. The time step of 10 was achieved with 100 iterations and a threshold of one-tenth. The remaining parameters were consistent with the problem statement and default file.

In this research, the post-processing flow model MODFLOW, known as MODPATH, will be used as the basis of calculations, and the sensitivity analysis of pollution concentration will be conducted using one-dimensional (up to multi-dimensional) theoretical dispersion methods. Essentially, in the execution of numerical models for pollutant transport analysis, the selected parameters in the Advection and Dispersion packages can significantly deviate the results from the ground reality. In this study, by employing a theoretical approach in extensive page environments, an attempt was made to measure the impact of the selected advection approach on the structure of the finite difference network model. The particle dispersion behavior theoretical investigation was also conducted using the Random Walk and probabilistic approaches.

Results

Investigation of pollutant dispersion coefficients using the explicit and implicit methods

A portion of the plain in the southern boundary was examined to analyze the sensitivity of pollutant dispersion in a porous aquifer environment and the variations in the coastal drainage system through finite difference simulation. This examination utilized the Explicit and Implicit methods and the random walk particle transport method. The primary objective was to investigate the sensitivity of the unvalidated qualitative model results and assess the error of the proposed model by examining the environment adjacent to the contaminated area, specifically the Landfill, through a limited test.

The scatter plots of the dispersion model at time steps 0 and 10 are shown in the following figures. The time step was set at 10, with 100 iterations and a threshold of 0.1. The remaining parameters were consistent with the problem statement and default settings.

Based on this method, the first step is to determine the coordinate position of the Mass Moment center. The concept of the mass column center, which refers to the center of the particle column, can be obtained by calculating the first moment of mass. The equation for this is as follows:

$$\frac{1}{{M_{{0}} }}\sum\limits_{i}^{m} {n \cdot C_{i} } \cdot x_{i} \Delta x\Delta y\Delta z$$

(1)

where M₀ is the total mass (M), n is the density (practical) (−), C_i is the concentration at node (ML-3), and x_i is the coordinate of x at the i^th node (L).

The RW method or random walk derivative can determine the particles’ coordinates (x and y). The averages of $\frac{1}{nr}\sum_{i}^{nr}{x}_{i}$ and $\frac{1}{nr}\sum_{i}^{nr}{y}_{i}$ It can be employed for the x and y coordinates, respectively. In the graph of this method, the mass center can be represented by a schematic square. The coordinates are presented as follows, based on the problem statement (Table 1):

Table 1 Position of the mass center of the plume; first mass moment

When developing the approximation graph for particle dispersion, linear regression is employed to determine the slope of the variations. The spread of plumes, in the general one-dimensional case, is dependent on the standard deviation of the individual particle coordinates. The standard deviation in the y direction is similar to that in the x direction. This value is calculated as follows:

$$\sigma_{x} = \sqrt {\frac{1}{{{\text{nr}}}}\mathop \sum \limits_{i}^{{{\text{nr}}}} \left( {x_{i} - x_{{i - {\text{AVG}}}} } \right)^{2} }$$

(2)

A circle can be swiftly drawn around the center of mass to represent the standard deviation of the plumes, which is equivalent to the mass moment. The circle’s radius depicts the spread of the plumes, taking into account the average spread in the x and y directions (Table 2).

Table 2 Spreading of the plume; second Mass Moment

In theoretical reports, the displacement (∆x) resulting from the diffusion process (including propagation) is explained to be equal to $\sqrt{\frac{2D\Delta t}{R}}$. It is assumed that only the diffusion process is of interest, and based on this assumption, the fundamental equation of particle transport will serve as the basis for calculation. At the tenth step of the calculation, the value is 0.354. By dividing the σavg (average standard deviation) obtained at the end of computational step 10 by the number of steps, a displacement delta of 0.317 indicates a difference of < 0.04 (value). The reason for dividing the calculated σavg by 10 is to convert it into computational units by ∆t. Figures 4, 5, and 6 show that the calculations’ output considers the particle propagation condition using the Finite Differences method.

Fig. 5

Display the position and state of dispersion at the start of the calculations

Fig. 6

Display the position and state of dispersion at the end of the calculations

Figures 4, 5, 6, 7, 8, and 9 illustrate the execution of computations at different time steps, with a threshold change of 5 units. Figure 10 represents the variation graph for two intermediate computation cases, specifically at 12.5 and 37.5 units, for all the methods from the release point. This graph can examine the computational differences between the selected methods.

Fig. 7

Execution for explicit and backward methods

Fig. 8

Execution for explicit and central methods

Fig. 9

Execution for implicit and central method

Fig. 10

Compares propagation in the three specified methods within two intervals

In order to obtain reliable results, these experiments were conducted within the applet’s environment. It is assumed that the soil can absorb salts. The water concentration, denoted as “c,” represents the mass of the solute per unit volume of water. On the other hand, the absorbed concentration refers to the solute’s mass per unit of soil. The governing equation for this scenario, assuming a homogeneous porous medium, is as follows:

$$\frac{\partial \theta c}{{\partial t}} + \frac{\partial \rho s}{{\partial t}} = D\theta \frac{{\partial^{2} c}}{\partial x} - v\theta \frac{\partial c}{{\partial x}}$$

(3)

θ represents the porosity, ρ represents the bulk density, and v is the flow velocity. The term “D” stands for the longitudinal hydrodynamic dispersion coefficient. We assume that absorption can be modeled as a linear equilibrium process, which implies that s = Kdc, where “Kd” is referred to as the distribution or partition coefficient. The above-governing equation assumes the inlet boundary is at x = 0, and the system extends infinitely in the x-direction.

Initially, a uniform concentration of “C_i” is present throughout the system, which, in this case, is equal to 1. At the time “t = 0,” salinity with a concentration of “C_f” is introduced into the column and maintained for a duration of “t_p,” also known as the pulse duration (Fig. 11).

Fig. 11

Schematic representation of the applets model execution

The model settings in the applets are shown in Fig. 12, displaying the concentration values constrained along all longitudinal coordinates for two-time intervals: 12.5 and 37.5 (Fig. 13 and Table 3).

Fig. 12

Applets model settings for extracting data in the 12.5-time interval

Fig. 13

Output data of the applets model for the time intervals 12.5 and 37.5

Table 3 Validation error values of the applets model compared to Excel environment data

Therefore, the error or computational difference between the Excel method and the applets model can be considered minimal, as shown in Table 4. For example, the calculated error between the Explicit and Backward methods data and the applets model for the time interval 12.5 is approximately 0.086 concentration units, and for the time interval 37.5, it is approximately 0.068 concentration units. Hence, the computational difference is not significant.

Table 4 Standard input values for hydraulic conductivity parameter in GMS and applets

To develop a sensitivity analysis model based on recent execution results, the value of the hydraulic conductivity coefficient in a test section was input into the flow model using the applet’s exponential distribution. However, modifications are required to accommodate the natural hydraulic conductivity GMS requires, as the applets utilize normal conductivity distributions for input. The system input modifications include the following:

$${\text{Mean}} = \ln \left( \mu \right) - \frac{{{\text{sd}}^{2} }}{2}$$

(4)

$${\text{sd}} = \sqrt {\ln \left( {1 + \frac{{\sigma^{2} }}{{\mu^{2} }}} \right)}$$

(5)

These averages and standard deviations (SD) can be utilized in the program. The numbers 15 and 3 were converted to 2.688 and 0.198, respectively, for the mean and standard deviation in the model, using μ as the mean of k and σ as the standard deviation of the k distribution (actual). Table 4 illustrates this conversion.

Figure 14 depicts the applet’s environment for determining the spatial distribution of hydraulic conductivity values (prior to conversion) in an area measuring 1500 m by 500 m. The output values of network cells were transformed into actual hydraulic conductivity values. This conversion and calculation process is described below.

Fig. 14

Applets settings for a more realistic hydraulic conductivity coefficient shape determination

Upon applying realistic coefficient values to the flow model, variations in hydraulic conductivity were obtained, as shown in the figure below. The water level at the western boundary remained constant at 10 m, while the average groundwater level in the region was 9.25 m. The specifications for the water level and other parameters at the eastern boundary were consistent with the default river model, with a water surface elevation of 6.75 m, a bed elevation of 4.75 m, a cell width of 10 m, and a riverbed resistance of 5 days.

The mass loading value, which indicates pollution concentration, was set at a constant value of 3500 g/l per 100 m² throughout the entire landfill area on the first day. The flow direction in this model indicated an expected movement from west to east. Thus, a sample well was placed in a cell in the eastern part of the landfill area to define the first scenario for pollution mitigation through extraction wells.

Based on these variations, it can be determined that if a significant amount of pollution remains at the end of the 2-year model execution, additional wells, such as the second, third, and so on, can help eliminate the remaining pollution.

The presented model effectively demonstrates the impact of a well with the desired flow rate on pollution dispersion. The pollution calculations generated by the MT3DMS model indicate that this well has contributed to removing 459,826 g of infiltrated pollution from the landfill area into the aquifer during the 2-year simulation. However, a significant amount of pollution remains stored in the aquifer. By examining the pollution contour plots of neighboring cells, it can be concluded that although pollution reduction in the aquifer is more significant when moving toward the southwest (presumably due to closer proximity to the pollution source), in one case, the desired well should still be located within the landfill area. As a result, cell 01-82-20 is considered preferable over other cells since it is positioned in an average southwestern location.

Figure 15 illustrates the positions of two additional wells incorporated into the original sound system to mitigate pollution through pumping. The placement of these two wells (specified by cell coordinates) was determined by referencing the previous model implementation, which identified a single well as the primary source of pollution even after 2 years. Consequently, despite a single well in the initial cell, the model indicated the continued existence of two pollution plumes. These two pollution plumes were subsequently chosen as the locations for the two additional wells, following the same methodology for constructing the original well. The figures provided indicate the cell coordinates for these newly established wells. Moreover, the impact of groundwater flow from these two wells can be observed after 2 years, equivalent to 730 days.

Fig. 15

illustrates the impact of three drilled wells on the groundwater table

Development of a realistic model through optimization of pollution dispersion coefficients

Figure 16 illustrates the flow balance sheet for all specified elements within the defined range, assuming the presence of three wells at the end of the modeling period. According to this assumption, 450 m³ of water is extracted from the aquifer daily. The water is discharged into the Landfill through wells with optimized coordinates. It can be observed that at the beginning of the period, with the entry of heavy pollution of 3500 g (on the first day) from each 100 m² cell into the aquifer, the entire landfill area becomes saturated with particles. Therefore, in the sensitivity analysis of particle transport from the landfill area toward the permeable boundaries of the Babol aquifer and the route based on the MODPATH model, we can rely on the findings presented in the finite difference model. Based on the provided sensitivity analysis, the method used in the final and realistic model is the isotropic linear equilibrium, considering the effect of radioactivity decay or degradation factors. The input parameter values in this package were assigned for a bulk density of 53,500 mg/m³, an equivalent of 0.00000585 for the first sorption constant, a numerical value of 0.0001 for the soluble rate constant, and the same value for the adsorption rate constant. These values remain unchanged based on the software default.

Fig. 16

Groundwater flow balance over 2 years considering pollution mitigation factors

Figure 17 shows the location of the Babol landfill outside the boundaries of the finite difference model cells. Therefore, it was necessary to define the distance between the Landfill and the model using the boundary condition in the conceptual model development of particle transport. In the particle transport simulation, the pollution source is the drainage network and permeable boundaries that pollution will directly affect during rainy periods, particularly during flood events. The concentration limit was defined as variable numbers based on recharge and precipitation within the range and set equal to the initial sustainable threshold of the Landfill itself in the model.

Fig. 17

Implementation of the simulation of polluted particle movement in the limited difference network

The pollution cloud concentration diagram demonstrates a consistent increase in concentration throughout the predicted period in the qualitative modeling of the plain. Initially, the rate of pollution increase was relatively low as the wastewater disposal site was outside the plain. However, the concentration accumulation can be attributed to permeable boundaries and a substantial decline in groundwater levels within the area. In other words, as the groundwater level significantly drops in the saturated zone of the plain, it transports a higher volume of polluted water toward the aquifer at an accelerated rate. The initial months of the prediction period play a crucial role in introducing pollution into the area. The crisis will inevitably escalate without artificial feeding programs or significant groundwater resource restrictions. Additionally, the reduced access to water will lead to an increased reliance on surface-feeding resources, resulting in the depletion of exploitation sources and the densification of urban areas downstream of this region, ultimately leading to a modified state.

By implementing the quantitative and pollution transport models, the results indicate that the pollution cloud expands to the locations depicted in Fig. 18. Considering the soil structure, slope, and flow velocity within this plain, three sections can be observed where the pollution cloud accumulates and expands, with the most significant dispersion occurring in the southern region.

Fig. 18

The final state of excessive pollution in the landfill damage zone wide display

The pollution level recorded at the final phase of the borders with a constant load is zero. However, 91,802.98 units have entered the sea. Approximately 3,331,729,000 pollution units and the extracted water have been extracted from the wells, contributing to pollution expansion. The waterway network has released 1,177,868,000 pollution units over 2 years, with 7,508,274,000 units discharged from the same network. Furthermore, 3,991,218,000 units are discharged from the permeable borders. It should be noted that these statistics can be reversed depending on the specified period. In conclusion, the difference between the input and output of pollution indicates that 8192 units entered the Landfill Damage Zone (Fig. 19).

Fig. 19

Statistical pollution dispersion summary at the modeling period’s end

By confirming the transfer of particles to the sea border, we consider the possibility of their return to the aquifer, which is equivalent to the return of seawater or the penetration of saline water into the aquifer.

In the two runs, where the second case reduced the visibility of seawater, the possibility of water penetration was examined. A numerical investigation of salinity diffusion using the Budget engine can determine changes in mass and concentration in the second model. Figure 20 shows these changes, indicating that, based on the flow, 3,069,952 kg of salinity has been reduced from the overall aquifer range in the final step of the model period. This means that 835,584 kg of salinity has accumulated in the aquifer compared to the initial model in this step. It should be noted that although this figure is less than one percent of the total mass of salinity in the modeling, if the scenario of reducing the Caspian Sea level is not implemented, a much more significant accumulation of salinity would occur. During this period, only 16,131.25 kg of salinity from the sea entered the aquifer range (all layers). However, in the realistic representation of the Caspian Sea level decrease, the decrease in sea level had a significant impact. The figure shows a decreasing trend of 1,595.07 kg of salinity penetrating the aquifer only in the final step.

Fig. 20

Implementation of the saline water infiltration model

Discussion and conclusion

This study developed a conceptual model using a finite difference groundwater flow model to evaluate the uncertainty of pollution dispersion in a porous aquifer environment. The model was developed utilizing the best software approach and processed raw data. Considering an automated calibration and validation approach, the groundwater flow model produced only minor statistical deviations from the optimized parameters. The particle transport model was also accurately implemented based on the predicted groundwater flow model. In the Babol aquifer, the waste disposal area was determined to be saturated outside the aquifer. The concentration at the aquifer’s boundaries and the discharge network were estimated by adjusting the simulation parameters for the particle transport of pollutant plumes. This implies that in the particle transport simulation, the pollutant’s spreading source, the discharge network, and the hydraulic boundaries will directly affect pollution during precipitation periods, especially flood events.

The concentration limit was represented as variable numbers based on recharge and precipitation within the range equivalent to the initial stable threshold and provided to the model for the landfill area. The cumulative concentration output graph at the boundaries of the aquifer and the river in the southeastern region also indicates the trend of changes in pollutant concentration in the groundwater quality modeling, which increases linearly during the prediction period. Although the initiation of pollution entry into the area is noteworthy from the early months of the prediction period, the expansion of the crisis will be inevitable without artificial recharge programs or significant restrictions on groundwater resources. Furthermore, extraction through exploitation sources and urban density downstream of this area, in the event of reduced access to water, will increase the utilization of surface water resources, resulting in modification.

The spatial representation of the pollution zone in the form of a finite difference network demonstrates excessive entry of pollutant volumes from several cells in a specific area in the southern part, near the location of waste disposal.

Although the concept of blocking permeable boundaries is theoretically proposed, it does not guarantee the delay or prevention of contamination transfer to the aquifer reservoir. This representation also demonstrates that the direction of groundwater flow from southern regions toward the northern plains exposes areas with lower elevations to a higher risk. Urban and rural land uses are practically situated along the contamination transport path due to their reliance on water resources. As a result, the Babol Plain is on the verge of a significant water quality crisis.

Although the particle transport model implementation assumes a homogeneous introduction of contamination at the boundary and a uniform flow expansion during the simulation period, the output of contamination dispersion has been randomly distributed and unpredictable in certain sections. This means that events such as climate change or hydrological network destruction in future periods can be expected to involve more sectional areas in contamination propagation.

The difference between the total inflow and outflow from the area was calculated as the particle dispersion factor to compare the flow balance in two sections. Based on the results, it can be observed that the maximum difference in flow between the two models is up to 900 m³/day, which decreases at the end of the simulation period. Similar calculations for the western outlets indicate a completely reversed trend compared to the eastern areas, which is undoubtedly due to the higher elevation of the drainage area.

Finally, taking into consideration the objectives of this study, which aim to analyze the sensitivity of pollutant dispersion in a porous aquifer environment and examine the changes in the drainage pattern of the coastal region obtained from the finite difference simulation model, we investigated the Explicit and Implicit methods, as well as the Random Walk Particle Tracking method, for a specific area of the active plain region. One of the main objectives was to assess the sensitivity of the uncalibrated qualitative model results and investigate the environment adjacent to the contaminated area, the Landfill, to evaluate the proposed model’s error through a limited test. Furthermore, the sensitivity analysis of the laboratory model, based on theoretical computations of pollutant dispersion, indicates an expected eastward flow direction in this modeling. In Scenario 1, the pollution was eliminated through Extraction Wells, where a hypothetical well sample was positioned east of the Landfill region to assess the spread rate of pollution plumes. The positions (cell coordinates) of these two wells were chosen based on the locations where, in the previous model execution with only one well, the highest pollution was still present after 2 years. This implies that although the model, with one well in the initial cell, eventually released a certain amount of pollution and created two persistent pollution plumes, these plumes were identified by the method used to construct the first well for the other two wells.

The impact of groundwater flow can also be observed from the two wells after 2 years (730 days). After analyzing the flow budget for all defined elements within the area and assuming the presence of three wells at the end of the modeling period, it was determined that 450 m³ of water per day is being extracted from the aquifer. The water extracted from the wells, optimized for their coordinates, is discharged into the landfill leachate.

At the beginning of the period, it can be observed that with the introduction of a heavy pollution load of 3500 g (on the first day) from each 100 m² cell into the aquifer, the entire landfill area becomes saturated with particles. Therefore, in the sensitivity analysis of particle transfer from the landfill area to the boundaries of the Babol aquifer, and based on the routing using the MODPATH model, the results obtained from the finite difference model can be considered reliable.

The results of this research can contribute to a more accurate prediction of future aquifer conditions under probable scenarios and facilitate planning and management with minimal adverse effects on water quality. Regional managers can utilize it as a decision-making mechanism for sustainable performance with lower uncertainty. This study employed a theoretical model to address the current uncertainty in the base groundwater flow model, which serves as the foundation for particle transport modeling. The model assumes a scaled small-scale distribution of structural parameters of the plain, such as hydraulic conductivity. The implementation of the model was based on the Exp and Imp methods, the RWPT method, and the App model, resulting in similar patterns of pollutant dispersion as the original model. Thus, while confirming the model’s accuracy, scenario-based particle discharge was also examined by providing the possibility of defining hypothetical wells.