Aquifer-wide estimation of longitudinal dispersivity by the combination of empirical equations, inverse solution, and aquifer zoning methods

Longitudinal dispersivity is a crucial parameter for the numerical simulation of groundwater quality, which is highly variable. The use of empirical equations and the inverse solution are two main methods of estimating longitudinal dispersivity. This study investigates the estimation of value and aquifer-wide spatial distribution of longitudinal dispersivity parameters using a combination of the empirical equation, the inverse solution method, and the aquifer zoning technique. The combined approach is applied to Bandar-e-Gaz aquifer in northern Iran, and Willmott’s index of agreement was used to assess the simulation precision of total dissolved solids in this aquifer. The values of this criterion were 0.9985–0.9999 and 0.9756–0.9992 in calibration and validation periods showing the developed combined approach obtained high precision for both calibration and validation periods, and the simulation shows remarkable consistency. Also, the one-way sensitivity analysis indicates that the longitudinal dispersivity is more sensitive than the effective porosity in this simulation. The investigation of the spatial distribution of the estimated longitudinal dispersivity by the combined approach indicates that the value of the parameter has a decreasing trend from the south to the north (50–8 m) in the aquifer environment, which is consistent with the changes in the characteristics of porous media in this study area. Therefore, it concludes that the combined approach provides a reliable and appropriate estimation of the spatial distribution of longitudinal dispersivity.


Introduction
The quality of the groundwater reduces due to the various contaminants entering the aquifer from different sources. The simulation and prediction of the temporal and spatial distribution of pollution dispersion in a groundwater aquifer is a good tool for managing the consumption of these resources (Rossetto et al. 2018). Simulation of groundwater quality could be considered an appropriate tool for aquifer management and planning (Ansarifar et al. 2020a, b). Typically, the governing equation (the related partial differential equations) is used to predict contaminant dispersion's spatial and temporal distribution in an aquifer (Singh et al. 2021).
These governing equations are often solved using a numerical approach by finite element or finite difference methods. They can estimate the variations in concentration in aquifers over time (Ding et al. 2021). Different models have been developed to estimate the variation of the spatial and temporal distribution of the contaminant in the aquifer environment. The Modular 3-Dimensional Transport model with the Multi-Species structure (MT3DMS) is one of the developed models in this field, which researchers have used considerably (Adeoye et al. 2018;Ameur et al. 2021). The MT3DMS model is used in conjunction with the MODFLOW model.
In other words, the MODFLOW model simulates flow behavior in an aquifer environment. Then, the MT3DMS model can be used to predict the temporal and spatial pattern of the contamination in an aquifer. Different studies have been done on applying the MT3DMS model to predict contamination changes in an aquifer in various regions. The MT3DMS model was used for the simulation tracer movement in the LaSalle Beauvais aquifer located in northern Paris (Zghibi et al. 2016), prediction of salt concentrations 1 3 14 Page 2 of 12 to current and future development conditions in Dehloran in the western region of Iran (Ehtiat et al. 2018), regional groundwater quality management study conducted in West Faisalabad of Pakistan (Shakoor et al. 2017). Also, it was used by El Khattabi et al. (2018) to study the impact of fertilizer on the groundwater quality in northern France, Golzar et al. (2018) to analyze of FE3O4 transport for twodimensional porous media, Sidiropoulos et al. (2018) for nitrate transport modeling in Lake Karla watershed as an over-exploited aquifer located in Greece and Khayyun et al. (2018) for modeling of the radioactive Cobalt-60 migration from LAMA nuclear facility in Iraq. There are different developed models to simulate the behavior of the contaminant in groundwater. The models require longitudinal dispersivity estimation since they include the application of an advection-dispersion equation plus a numerical solution. Different researchers have investigated methods for the assessment of longitudinal dispersivity. Pickens and Grisak (1981) investigated the longitudinal dispersivity magnitude in a stratified granular using laboratory column and field tracer examinations. They reported that dispersivity is a linear function of the mean travel distance. Arya (1986) investigated the variation of macroscopic and megascopic dispersions versus heterogeneity, diffusion coefficient, autocorrelation, and aspect ratio. He concluded that the dispersivity increases with travel distance, and the diffusion is insignificant in the field scale. Neuman (1990) suggested a universal scaling law for dispersivities in geologic media by treating the medium as uniform. The weighted least-squares technique was applied to estimate dispersivity as a function of the field scale by Xu and Eckstein (1995). The empirical equations provided by these researchers have been used in many studies. The Arya equation was applied by Salmo et al. (2022) to model unstable water-oil displacement by Adegbite and Al-Shalabi (2022) to study the effect of heterogeneity on engineered water injection in carbonates. Also, the Pickens and Grisak equation was used by Di Gianfilippo (2018) to evaluate the recycling potential of alkaline waste as a filler for road sub-bases and by McLean et al. (2019) for statistical modeling of groundwater contamination monitoring data and by Madie et al. (2022) for modeling pollutants in a porous medium with finite thickness. The Neuman empirical equation was applied to simulate sewage leakage and transport in a carbonate aquifer (Dvory et al. 2018) and to estimate macrodispersivity for groundwater transport (Zech et al. 2022). The Xu and Eckstein empirical equation was considered for chloride transportation in urban riparian aquifers (Ledford et al. 2016), estimating anisotropic heterogeneous hydraulic conductivity and dispersivity (Priyanka et al. 2018) and three-dimensional modeling of the heterogeneous coastal aquifers (Priyanka and Mohan Kumar. 2019. The expressed empirical equations for estimations of longitudinal dispersivity in the aquifer are widely used, but the following significant constraints influence these empirical equations.
1. These equations use only one scale variable to estimate the longitudinal dispersivity. At the same time, factors such as particle size mean and range and uniformity coefficient can also significantly affect it. 2. Estimated values of the longitudinal dispersivity by different equations in most cases have a significant difference. 3. Different reports show that longitudinal dispersivity values estimated in aquifers are considerably greater than those estimated by empirical equations.
Appropriate spatial estimation of longitudinal dispersivity is essential for conducting accurate and reliable water quality simulations of groundwater resources. The inverse solution method is also suitable for parameter estimation in different models, which can also be used to estimate longitudinal dispersivity in groundwater quality studies. This method uses an optimization algorithm to minimize the simulation error and estimate the parameter, and the results may not be appropriate and reasonable. An aquifer zoning technique can be used in which the application of this technique can reasonably help achieve proper and reliable parameter estimation. This study aims to evaluate the efficiency of a combined approach by combining the inverse solution and empirical equations methods and aquifer zoning technique to estimate the spatial distribution of the longitudinal dispersivity parameter. The Bandar-e-Gaz aquifer in northern Iran is studied to investigate the results of this combined approach.

The study aquifer
The study area is Bandar-e-Gaz unconfined aquifer in Iran, consider ing geog raphic coordinates of 53 °50′13"-54 °05′10" E longitude and 36 °41′48"-36 °48′29" N latitude. The Gorgan gulf is located in the north of this aquifer, which causes the freshwater resources of the Bandar-e-Gaz aquifer to be exposed to saltwater intrusion (Ansarifar et al. 2020a, b). There are also many towns and villages in the Bandar-e-Gaz aquifer area, which generally use absorption wells for sewage disposal, so the freshwater resources of this aquifer are seriously at risk of water quality degradation. The size of the Bandar-e-Gaz aquifer deposits is reduced from south to north so that in the southern regions, especially in the alluvial fans regions, the coarse-grained sedimentary deposits include rubble, gravel, and sand. However, the fine-grained sedimentary deposits in northern regions include fine sand, silt, and clay.
The following figure presents the geographic location of the Bandar-e-Gaz aquifer.

Inverse solution method
A combination of parameters and variables should be used for quantitative and qualitative groundwater modeling. Model parameters are usually unknown and are estimated using the inverse solution method. The inverse solution method tries to estimate the proper values of the model parameters by considering the model's variables and the initial values for the parameters by reducing the error between observed and computed values. The process of reducing errors between observed and computed values is performed using optimization algorithms. The MODFLOW model uses the PEST algorithm for optimization, which can estimate the values of parameters such as specific yield and hydraulic conductivity. In the MT3DMS model, the inverse solution method is trial and error, or in other words, non-automated. The objective function in the inverse solution method in MODFLOW and MT3DMS models is to minimize the rootmean-square error (RMSE) function considering the logical values and behavior of different variables.

The numerical models:
The MODFLOW groundwater model was applied for groundwater level simulation, and the MT3DMS model, in the form of an integrated model with MODFLOW, was used for aquifer-wide simulation of spatial distribution and concentration variation of total dissolved solids (TDS). A proper cell grid must be created following the conditions of the aquifer to simulate using the MODFLOW and MT3DMS models (Zong et al. 2021). Various packages, including boundary conditions, quantitative and qualitative observation wells, hydrodynamic parameters (hydraulic conductivity, specific yield, longitudinal dispersivity, and effective porosity), recharge, and discharge, must be prepared to simulate groundwater's quantitative and qualitative behavior. The three-dimensional groundwater flow follows the governing equation solved by the MODFLOW model. The unknown velocity and head in each cell are computed at a point or node at the center of the cell (Lalehzari et al. 2013).
In the mentioned equation, K zz , K yy , and K xx are values of principal hydraulic conductivity along the z, y, and x coordinate axes (LT −1 ), respectively; W is a volumetric flux per unit volume introducing sinks and/or sources of water (L 3 T −1 ); h is the potentiometric head (L); S s is the specific storage of the porous material (T −1 ); and t is time (T). The MT3DMS model solves numerically the following partial differential equation governing the 3-D solute transport of a single chemical constituent in groundwater, considering chemical reactions, dispersion, and advection: where C is the dissolved contaminant concentration (ML −3 ); Dij is the hydrodynamic dispersion parameter tensor (L 2 T −1 ); t is time (T); q s is the volumetric flow rate per unit volume of the aquifer and represents fluid sources and sinks (T −1 ); Vi is the pore water velocity (LT −1 ); Cs denotes the concentration of the fluid source or sinks flux (ML −3 ); θ is the effective porosity; and R k is relative to chemical reactions. The MODFLOW and MT3DMS models implement a finite difference method to numerically solve the partial differential equations (Abraham and Manikannan, 2021).

The model precision criterion and sensitivity analysis
The aquifer-wide simulation of groundwater level and total dissolved solids is conducted for two validation and calibration periods. In this study, Willmott's index of agreement criterion (Willmott, 1981) was applied to evaluate the precision of the results. The equation for Willmott's index of agreement criterion is presented below, and this criterion has a value between 0 and 1, the value of one representing the best fit (Ghadim et al. 2020;Guo et al. 2019).
The used parameters in different models and various simulation processes do not have the same impact magnitude, and it is necessary to estimate the significance of each of them. Sensitivity analysis is a suitable method and test for measuring the effectiveness of the parameters used in the model on model accuracy. The one-way sensitivity analysis was used to measure the model's accuracy by considering each parameter's increasing and decreasing changes relative to the optimal state (Romeo et al. 2022).

The field data, aquifer zoning technique, and empirical estimations of the longitudinal dispersivity
The field data and aquifer zoning technique There are eleven and four groundwater levels and quality observation wells located in the aquifer-wide environment with recorded monthly water level TDS, as presented in Fig. 1. For this research, 24-month recorded data of these observation wells in the time domain of (October 2011-September 2013) are used, and the statistical characteristics of these data are presented in Table 1. The aquifer environment is different from laboratory conditions and has variations in the various characteristics of the porous media. Therefore, the aquifer study as a zone can lead to illogical estimates. By dividing the aquifer into different zones, one can find a more suitable and reliable assessment of the behavior of the aquifer and its characteristics. According to available data, the division of the aquifer into different zones is subject to a practical limitation. The aquifer was divided into four zones considering the locations of water quality observation wells, the aquifer-wide pattern of sediment size variations, and the direction of groundwater flow. In Fig. 2, the location of these zones is presented.

Estimation of the longitudinal dispersivity using empirical estimation
Since molecular diffusion and mechanical dispersion processes cannot be separated in aquifer-wide groundwater flow, a single dispersion parameter named the hydrodynamic dispersion parameter, or, in short, the dispersion parameter, is used in groundwater studies. This dispersion parameter has the effects of both molecular diffusion and mechanical dispersion. It is the combined output of two porous media's mass transport processes. Both molecular diffusion and mechanical dispersion will force the contaminant to move longitudinally, taking into account the mean flow direction. In addition, the longitudinal dispersion parameter in the porous medium can be written as follows:   For m = 1, the al parameter has a length dimension. In this case, it is called the longitudinal dispersivity parameter, and it expresses the ability of the porous media to mix and transport the contaminant. Researchers have presented various equations for estimating longitudinal dispersivity based on laboratory and field data analysis. In these studies, field scale is defined as "the distance traveled from the source for ambient conditions or the distance between injection and observation wells for the case of an induced flow configuration." Gelhar et al. (1992) carefully classified 106 data points into three classes of varying reliability levels: high, intermediate, or low. The data points and their associated reliability are presented in Fig. 3.
The equations presented below are the most important empirical equations for estimating longitudinal dispersivity, widely used in various studies. Pickens and Grisak (1981): Arya (1986): Neuman (1990): Xu and Eckstein (1995): 10) a L = 0.32L 0.83 100 m < L where L is traveled distance and al is longitudinal dispersivity. Experimental results show that longitudinal dispersivity mainly depends on the range and mean particle size, the uniformity coefficient, and some other characteristics of porous media (Sun and Sun 2013). The longitudinal dispersivity becomes greater with the uniformity coefficient and varies from laboratory to field experiments. The explanation for much larger longitudinal dispersivity in a field test is that the uniformity coefficients of the media in the field are larger than that in the laboratory (Sun and Sun 2013). Other reports, such as Sudicky et al. (1983), Anderson (1979), andFried (1975), indicated that longitudinal dispersivity values estimated using field examinations are several orders of magnitude higher than the laboratory ones. In empirical equations, the magnitude of longitudinal dispersivity only depends on the scale of experiments, while the impacts of other factors have not been considered. Also, the largest value for higher reliability of studied data for empirical equations extraction was only 250 (m). These limitations result in the estimates made by the empirical equations being merely an initial value for longitudinal dispersivity and the acceptable and logical values estimated by the researchers by additional investigations.

Combination of the inverse solution, empirical equations methods, and aquifer zoning technique
The models for the simulation of groundwater behavior often include parameters that cannot easily be measured in (11) a L = 0.83(log L) 2.214

Fig. 3
Longitudinal dispersivity against scale with data classified by reliability (Gelhar et al. 1992) nature and are not available in most practical studies (Hamraz et al. 2015). Typically, the values of used parameters in groundwater models are adjusted non-automatic (trial-anderror method) or automatically by optimization algorithms. Non-automatic calibration methods include estimating the parameter values by the modeler, and the success of this method is entirely dependent on the modeling experience and knowledge (Sadeghi-Tabas et al. 2017). Implementing automated calibration methods of the models is also easy and widely used, so automatic optimization algorithms have been developed in recent decades. The automatic optimization methods to find the global optimal value of the objective function in the range of parameter variations face severe constraints, and the estimation of the parameters by different optimization algorithms is accompanied by mistakes (Miao et al. 2019). Also, suppose a series of different parameters have the same performance for modeling during the calibration period. In that case, there will be no criterion for choosing the best parameter series, and predictions outside the calibration period from this series of parameters will not give similar results. The distributed numerical groundwater models require estimating the spatial distribution of the aquifer's hydrodynamic parameters. These hydrodynamic parameters are subject to aquifer physical conditions and vary in aquifer space.
On the one hand, using the aquifer zoning technique can lead to an improved logical estimation of the spatial distribution of hydrodynamic parameters and more reliable results. On the other hand, applying field data and empirical equations developed in estimating hydrodynamic parameters can lead to an appropriate initial estimation of these parameters. In the approach developed in this study, the aquifer environment is zoned based on available data. Then, initial estimates of the hydrodynamic parameters are obtained using the experimental equations for different zones. In the following, the value of the parameters in each zone is estimated using the inverse solution method. After implementing this combined approach, the parameter estimation results should be evaluated, and the fitting of the values and the spatial distribution of the estimated parameters to the field and physical realities in the region must be reasonable and appropriate. The flowchart of this combined approach is presented in Fig. 4.

Results
The results of this study are presented in the following three sections. Quantitative modeling is the base of qualitative modeling; therefore, high-precision quantitative modeling can have a perfect effect on the results of qualitative modeling. The results of quantitative modeling are presented in the first section. The results of aquifer environment zoning and empirical equations were introduced in the second section. The results of the combination of the inverse solution, empirical equation methods, and aquifer zoning technique in estimating the spatial distribution of longitudinal dispersivity were introduced in the third section.

Results of quantitative modeling
In the process of quantitative modeling, the groundwater level simulations in both calibration and validation periods were conducted, and Willmott's index of agreement criterion was applied to evaluate the precision of the simulation results. Their values in calibration and validation periods are provided in Fig. 5. An assessment of this plot shows that simulations had high precision in both calibration and validation periods. In the calibration period, the hydraulic conductivity and specific yield parameters were estimated using an inverse solution method, where the spatial distribution of specific yield is presented in Fig. 6. The investigation of this map shows that the aquifer-wide range of specific yield in Bandar-e-Gaz aquifer is (0.022-0.038) and has a decreasing trend from the south to the north.

Results of the longitudinal dispersivity estimation based on aquifer zoning technique and experimental equations
The longitudinal dispersion is a function of molecular diffusion, longitudinal dispersivity, and flow velocity in the porous media. The flow velocity can be estimated using the hydraulic gradient, hydraulic conductivity, and effective porosity. In this study, molecular diffusion has no significant weight compared to mechanical dispersion and is neglected. The most critical parameter to be estimated is the longitudinal dispersivity, and the aquifer environment was zoned for the initial estimation of this parameter. Regarding the number and location of water quality observation wells, the aquifer-wide changes in sediment size, and the direction of groundwater flow in the aquifer environment, the Bandar-e-Gaz aquifer was divided into four zones. The values of the longitudinal dispersivity in each zone are estimated using the different empirical equations and presented in Table 2.
The estimated values for the parameter of longitudinal dispersivity by using empirical equations led to different estimations. It is necessary to note that the longitudinal dispersivity parameter in nature has highly variable behavior. The lack of access to an appropriate initial approximation of this parameter can lead to significant mistakes in simulating the groundwater quality components.

Results of the spatial estimation of longitudinal dispersivity by the combination of aquifer zoning technique, experimental equations, and inverse solution
The empirical equations based on estimated values of the longitudinal dispersivity parameter are used to approximate the initial value and the logical range of this parameter for calibrating the MT3DMS model. Using the trial-and-error (non-automated) calibration of the MT3DMS model, the longitudinal dispersivity parameter in the aquifer environment was estimated. The values of Willmott's agreement index were computed in both validation and calibration periods, as presented in Fig. 7, to check the precision of the simulation results. The values of this precision criterion imply that the results of both validation and calibration periods have perfect precision. In other words, the appropriate estimation of the longitudinal dispersivity parameter of the MT3DMS led to a consistent behavior for this model (Fig. 8).
The empirical equation method led to the various estimation of the longitudinal dispersivity parameter in different zones and the difference between hydraulic gradients, effective porosity, and hydraulic conductivity in the studied zones. Using a combined approach of the aquifer zoning technique, the application of empirical equations and inverse solution methods resulted in a reasonable and logical spatial distribution for this parameter following the hydrogeological characteristics of the Bandar-e-Gaz aquifer. The aquifer-wide values and spatial distribution of the longitudinal dispersivity parameter in the studied aquifer are shown in Fig. 9. The presented values in this map show that the longitudinal dispersivity parameter decreases from the south to the north, consistent with the sediment size variation pattern, the changes in the hydraulic gradient, and the effective porosity in this area. It can be deduced that the estimations of the longitudinal dispersivity by Arya and Xu & Eckstein equations are in accordance with the hydrogeological characteristics of the studied aquifer, considering the results of the combined approach. In contrast, the results of Neuman and Pickens & Grisak's equations are not acceptable.

Results of the sensitivity analysis
Sensitivity analysis of estimated parameters can be a suitable tool to examine the study's results. On the one hand, the  analysis results identify which parameters have more effects on the simulation. On the other hand, it can determine the magnitude of the impact of variation in optimized parameters on model simulation precision. The one-way sensitivity analysis is applied in this study for two effective porosity and longitudinal dispersivity parameters, and the results of this analysis are presented in Fig. 10. The sensitivity analysis plot indicated that longitudinal dispersivity is more sensitive than effective porosity in both decreasing and increasing directions. Moreover, according to this plot, the simulation is more responsive for longitudinal dispersivity in a declining direction than in a rising direction, but vice versa for effective porosity.

Conclusion
Hydrodynamic dispersion is the most crucial parameter in modeling groundwater resource quality. It depends on several factors, including molecular diffusion, flow gradient, hydraulic conductivity, and longitudinal dispersivity. Estimating the longitudinal dispersivity parameter is more important and less reliable than other factors. Two methods are widely used to estimate longitudinal dispersivity: the first method uses experimental equations, and the second method is the inverse solution. The most important advantage of the empirical equation method is that it is based on laboratory and field experiment data, but it also has some limitations. One limitation of the empirical equation method is that different equations lead to varying estimations for longitudinal dispersivity. Another one is that these equations only calculate the longitudinal dispersivity based on the travel distance of the contaminant. At the same time, this parameter is dependent on different characteristics of the porous media. By minimizing the simulation error, the inverse solution method results in the estimation of longitudinal dispersivity, but this method's results may also be problems. The estimated longitudinal dispersivity by the inverse solution method is not necessarily accurate. Another common problem is the spatial distribution of this parameter, which may not follow the characteristics of the aquifer. In other words, the results may not be reliable or acceptable.
In some cases, these problems cause the simulation results to be inappropriate in the validation period, and the model does not have the proper consistency. In this study, a combined approach, including a combination of inverse solution method, empirical equations methods, and aquifer zoning techniques, was developed to improve the precision of the spatial distribution of longitudinal dispersivity. The recorded data and information of the Bandar-e-Gaz aquifer in northern Iran are used to investigate the combined approach's efficiency, and the total dissolved solids' behavior in this aquifer's groundwater is simulated. In this developed approach, considering the location of observation wells, the direction of flow, and the pattern of changes in aquiferwide sediment size, the aquifer area was divided into four zones. The zoning of the aquifer was carried out so that one observation well was located in each zone, which increased the reliability of the combined approach. The values of longitudinal dispersion are estimated using the empirical equation for these four zones. They are used as an initial approximation and a reasonable logic range for applying the inverse solution method. The value of the parameter of longitudinal dispersivity in different zones was estimated using the inverse solution method, and its spatial distribution was determined. The simulation was conducted in calibration (October 2011-September 2012) and validation (October 2012-September 2013) periods, and Willmott's index of agreement was used for model precision analysis. The results showed that this criterion was in the ranges of (0.9985-0.9999) and (0.9756-0.9992) in two calibration and validation periods, respectively, and these results indicated a significant consistency of the model. The considerable precision and consistency of the model seem to be due to the appropriate estimation of the spatial distribution of longitudinal dispersivity. Moreover, the application of the combined approach showed that the result of the Arya and Xu and Eckstein empirical equations compared to other equations led to a better estimation of the longitudinal dispersivity in the aquifer. The sensitivity analysis of the qualitative model confirmed that longitudinal dispersivity has more sensitivity than effective porosity. Investigation of the spatial distribution of longitudinal dispersivity showed that the value of this parameter from the south side to the north (upstream to downstream) has a decreasing pattern, which conformed to the pattern of variations in the sediment size of the porous media. The proper accordance between longitudinal dispersivity's spatial distribution and aquifer characteristics confirmed that the developed approach has led to the reliable aquifer-wide estimation of values and spatial distribution.
Funding The author(s) received no specific funding for this work.

Data availability
The data supporting this study's findings are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.