Application of Monte Carlo Simulation to Find Travel Time of Groundwater in the Iraqi Western Desert

In this study, a computerized mathematical method represented by Monte Carlo simulation was used to predict the travel time of the groundwater flow in the Iraqi western desert. During the run of the simulations, all the hydraulic parameters of Darcy’s Law were fixed but the hydraulic conductivity. The input data of the hydraulic conductivity is compared to the triangular distribution function to find the best number of iteration to run the simulations. The results showed that an iteration number of 5000 was enough to achieve best match between the input data of the hydraulic conductivity and the fitted distribution function. In addition, the estimated travel time of the groundwater flow is broadly varied through the entire area and ranges from 1983 years to 113 741 years based on 10 000 m of travel distance. Furthermore, hydraulic conductivity of the aquifer has high impact on the estimated travel time of the groundwater flow. However, head difference of groundwater elevation among the selected wells considerably influences the expected travel time of the groundwater flow.


Introduction
Increasingly, reliable predictions of groundwater flow and other hydraulic properties associated with uncertainties have been widely sought.Hence, various stochastic models to evaluate the uncertainty of groundwater problems have been developed.Uncertainty is defined as a situation of having lack of confidence in future outcomes as a result of unknown or inadequate input variables (Singh, Jain, & Tyagi, 2007).In general, there are two types of method to analyze the uncertainty; analytical which uses models to propagate the uncertainty to estimate the probability distributions and empirical which uses computer simulation experiments to generate the distributions of the outputs (Doctor, Jacobson, & Buchanan, 1988).Uncertainty of a system can be caused by natural activities such as those of environmental, hydrological, metrological or other natural processes which cause significant fluctuations in the space and time.In addition, inadequate data for a number of parameters which are required for risk analysis is a major source of the uncertainty in the system (Singh et al., 2007).To find or predict groundwater flow of an aquifer, it is necessary to identify the hydraulic properties of this aquifer.The parameters such as hydraulic conductivity or transmissivity, dispersivity, porosity, and hydraulic gradient, which influence the groundwater flow and transport, are taken into account as the main variables in constructing a stochastic model.However, lack of data of the hydraulic properties as a result of the heterogeneity of the natural porous media leads to the uncertainty of the groundwater flow problems (Neshat, Pradhan, & Javadi, 2015).Additionally, it is well known that such parameters are characterized by spatial distribution throughout the aquifer and their values are randomly varied according to the point measurements (Crestani, Camporese, & Salandin, 2015;Hassan, Bekhit, & Chapman, 2009;Zhang & San, 2000;Hoeksema & Clapp, 1990).To produce a stochastic travel time of groundwater in a natural aquifer, a differential form of Darcy's law which is assumed to be valid for most groundwater flow is used (Zhang, Shi, Chang, & Yang, 2010;Fitts, 2002;Zhang & San, 2000).According to Darcy's law, travel time of groundwater flow is a function of travel distance, porosity, retardation factor, gradient, and hydraulic conductivity.Application of Darcy's law in groundwater flow problems and its advantages and limitations is well elaborated in literature.random samples mostly lies within the range of the input distribution.Monte Carlo analysis involves different types of distribution functions to generate a large number of random samples (Doctor et al., 1988).Some of the probability distributions are Triangular, Log Logistic, Beta General, Pearson5, Inve Gauss, and Log normal.In each distribution function, data can be simulated in different iterations up to 10 000 iterations.Several researchers reported different numbers of iterations to represent the sample population.For example, Driels & Shin (2004) indicated that running the simulation for 4684 iterations is accurate to 5% and gives a confidence interval of 95% for weapon effectiveness.Moreover, the number of the simulations also can be changed which results in a higher probability.The simulation is stopped when the desired confidence interval with a specific error has been achieved to produce a clear statistical description of dependent variables.However, using a larger number of simulations to achieve higher accuracy by generating more random numbers and increasing the sample size will generally increase the time consumed by the computer to finish the simulated run, which leads to high costs in complicated transient systems (Pasetto, Guadagnini, & Putti, 2011;Singh et al., 2007).The results of the simulated data can be compared in many ways, one of which is by using fit comparison to the input or output data.A fit comparison graph (density or cumulative curves) combines both the input data and the fitted distribution in one graph in order to determine the location(s) of the best match area of the fitted distribution to the input data.A higher number of iterations gives a closer fit to the data (Palisade Corporation, 2008).The number of simulations can be fixed whenever the difference between the results of two or more simulations becomes insignificant.
Over the years, application of Monte Carlo simulation methods to estimate parameter uncertainty in groundwater hydrology has been widely used.Hoeksema and Clapp (1990), for instance, applied a Monte Carlo simulation to perform model calibration of a groundwater flow.Fu and Gómez-Hernández (2009), in another study, used Monte Carlo Markov chains to assess the uncertainty of groundwater flow and mass transport of a synthetic aquifer with three variables; conductivity, piezometric head, and travel time.They found that all three parameters impacted the uncertainty, and the piezometric head had the highest effect on uncertainty reduction.Similarly, a Monte Carlo Markov Chain method was also used by Hassan et al. (2009) to reduce the range of the input parameter distributions of a two dimensional groundwater flow model.For more studies see Zhang et al. (2010), Hassan et al. (2009), andPasetto et al. (2011).
In this study, a computerized mathematical method represented by the Monte Carlo technique is used to estimate the likelihood of the travel time of the groundwater flow in several regions of the Iraqi western desert.The simulation included a sampling of hydraulic conductivity of the aquifers in the considered area.

Description of Study Area
The boundary of the considered area in this study is the Euphrates River from the east and north east and Iraq's border with Syria, Jordan, Saudi Arabia, and Kuwait from west and south as shown in Figure (1).The Iraqi western desert has several unconfined aquifers which are hydraulically connected to each other as well as some confined aquifers (Al-Fatlawi & Jawad, 2011).The acceptable quality of the groundwater which is easy to access due to its closeness to the surface of the earth makes these aquifers a potential good source of water in this area (Al-Mussawi, 2014).Topographically, the north part of the area which is mostly covered with sandstone, clastic, and marl has a slope of about 0.002 from west to east.However, the southern part of the area which is generally covered with residual soil with calcium carbonate and chert has a lower slope toward the east.Geologically, limestone, dolomitic limestone, dolostone, marly limestone, sandstone, siltstone, claystone, marl, and evaporate are the most common rocks of the aquifers in the region of the Iraqi western desert (Al-Jiburi & Al-Basrawi, 2009;Al-Jiburi & Al-Basrawi, 2007).The rocks' diverse geological formations in the aquifers resulted in a wide range in the hydrogeological properties and hence uncertain hydraulic conductivity throughout this region.Therefore, the hydraulic conductivity of the aquifers in the Iraqi western desert has a broad range which is between 0.03 m/day and 100 m/day.The static groundwater level ranges from several meters to more than 300 m below the ground surface.The main source of the recharge to these aquifers is rainfall.The flow of the groundwater in the designated area is in the direction of the east and northeast (from west and southwest towards the Euphrates River).Additional information on the studied area is described elsewhere (Al-Jiburi & Al-Basrawi, 2009;Al-Jiburi & Al-Basrawi, 2007).

Methodology
The available data of several wells located in the Iraqi western desert which are designated by the yellow color in Fig. 1 and listed in Table ( 1) as well as many existing wells in this region, which are described elsewhere (Al-Jiburi & Al-Basrawi, 2009;Al-Jiburi & Al-Basrawi, 2007), were used to run a computerized mathematical method represented by Monte Carlo simulation to estimate the travel time of the groundwater in this area.The locations of the selected wells represent nearly all the area of the Iraqi western desert.As previously mentioned, according to Darcy's law, travel time (T) is the function of the distance (L) for groundwater to travel, hydraulic gradient (dh/dl), porosity, and hydraulic conductivity (K), and is expressed in the following formula (Todd & Mays, 2005): In addition, it is well-known that the most affected parameter by the wide diversity of the natural lithologic material formation is the hydraulic conductivity the values of which cover an extensive range, even with very small hydrologic units.Therefore, in this study, all the parameters of Darcy's Law apart from the hydraulic conductivity were fixed.Again, Monte Carlo simulation can evaluate models with iterations up to 10 000 and it is used to analyze the uncertain propagation of a model.The simulation involved the use of triangular distribution function to form the statistical analysis of the hydraulic conductivity.Triangular distribution is an easy probability function that is used in hydrologic and environmental engineering applications (Singh et al., 2007).To use triangular distribution function, three values, minimum, most likely, and maximum, of the hydraulic conductivity are needed.The hydraulic conductivity of all wells close to the selected wells in Figure ( (2) At the beginning, two different wells (RW-5) and (B7-13) were used to run several simulations by using different numbers of iterations to achieve smooth hydraulic conductivity and determine the best number for further simulations.A fit comparison graph for each run was used to determine the best match between the input data and the fitted distribution.If the best match is around the mean or along with tails, good results are most likely achieved (Palisade Corporation, 2008).The best number of iterations was used to assess the travel time of the groundwater for all selected pairs of wells.

Results and Discussion
As previously mentioned, wells (RW-5) and (B7-13) were used to run several simulations with different numbers of iterations to determine the best number that can be used for further simulations.Histograms of Figures 2, 3, 4, 5, 6, and 7 show the simulations that have been run with iteration number of 100, 500, 1000, 2000, 5000, and 10 000, respectively, to estimate the hydraulic conductivity between well (RW-50) and well (B7-13).The simulations, as previously mentioned, involved the use of the triangular function as a fitted distribution to compare it to the input data of the hydraulic conductivity.However, in Monte Carlo simulation, travel time within given hydrologic parameters most likely follows a lognormal distribution (Doctor et al., 1988).Therefore, lognormal density was used as a fitted distribution to compare it to the input data of the travel time.Histograms of Figures 2 through 6 explained that when the number of iteration was gradually increased, an obvious increase in the match between the input data of the hydraulic conductivity and the triangular distribution function was observed and the best fit was achieved at an iteration number of 5000.However, increasing the number of iteration to 10 000 showed no significant difference to the 5000 iterations as shown in Figure 7 and compared to Figure 6.As a result, a number of 5000 iterations was used to run all the simulations to estimate the likelihood of the travel time of the groundwater flow of each selected pair of wells.Figures A1 through A8   (3).The units of the estimated travel time shown in Table (3) which are expressed in days were converted to years to make it easier in the comparison.In addition, a fit comparison graph to determine the best match between the input data of travel time (between well (RW-50) and well (B7-13)) and lognormal density as a fitted distribution is illustrated in Figure 9. Furthermore, Figures B1 through B8 (Appendix B) display fit comparisons for the selected pairs of the wells which are listed in Table (2) and were used to run the software of the Monte Carlo simulation.Table (3) explains that the predicted travel time of the groundwater flow is varied through the entire Iraqi western desert which is the result of the wide range in hydrogeological properties of the geological formation in the aquifers.For example, the estimated mean travel time of the groundwater flow of the north part between well (RW-50) and well (B7-13) (which is 33 682 years with travel distance of about 134 000 m) is about half of the estimated mean travel time of the southern part between well (5616) and well (4910) (which is 67 033 years with travel distance of about 127 080 m).This is attributed to the higher values of the hydraulic conductivity of the geological formation of the southern part of the selected area.Moreover, although the travel distance between well (5383) and well (5361) which is 139 000 m is higher than that between well (5383) and well (5351) which is 109 670 m, the estimated travel time of the groundwater flow between the first two wells (103 396 years) is much lower than that of the second two wells (1247 406 years).This is more likely because of the higher head difference of the elevation (97.8 m) between well (5383) and well (5361) compared to that of 4.5 m between well (5383) and well (5351) and as displayed in Table (2).Therefore, in addition to the hydraulic conductivity of the aquifers, head difference of the groundwater elevation highly influences the outputs of the simulations to predict travel time of groundwater flow.Over all, the histograms of the computerized mathematical method which is represented by Monte Carlo simulation revealed that the estimated travel time of the groundwater flow for the selected wells in the Iraqi western desert ranges from 1983 years to 113 741 years based on 10 000 m of travel distance.Moreover, fit comparison graphs to determine the best match between the input data of travel time and fitted distribution function show the same behavior for all selected wells.One of the advantages of assessing groundwater flow time to travel is to estimate travel time for a contaminant to reach discharge boundary.Therefore, it is recommended to apply Monte Carlo simulation on studying contaminant plumes in groundwater transport in this area.

Conclusions
Following conclusions are drawn from this study: 1) An iteration number of 5000 was enough to achieve best match between the input data of the hydraulic conductivity and the triangular distribution function.
2) As the hydrogeological properties of the Iraqi western desert are significantly varied, the expected travel time of the groundwater flow is widely varied as indicated by the histograms which resulted from Monte Carlo simulations.
3) As indicated by several studies, hydraulic conductivity of the aquifer has a high impact on the estimated travel time of the groundwater flow.However, head difference of groundwater elevation among the selected wells considerably influences the expected travel time of the groundwater flow.Fit Comparison for Travel Time, T (day) between well (5579) and well ( 5510) RiskLognorm(3050748.3,4449629

Figure 1 .
Figure 1.Location of the Considered Area (Adopted from Al-Jiburi & Al-Basrawi (2009) and Al-Jiburi & Al-Basrawi (2007)) Figures 2, 3, 4, 5, 6, and 7  show the simulations that have been run with iteration number of 100, 500, 1000, 2000, 5000, and 10 000, respectively, to estimate the hydraulic conductivity between well (RW-50) and well (B7-13).The simulations, as previously mentioned, involved the use of the triangular function as a fitted distribution to compare it to the input data of the hydraulic conductivity.However, in Monte Carlo simulation, travel time within given hydrologic parameters most likely follows a lognormal distribution(Doctor et al., 1988).Therefore, lognormal density was used as a fitted distribution to compare it to the input data of the travel time.Histograms of Figures 2 through 6 explained that when the number of iteration was gradually increased, an obvious increase in the match between the input data of the hydraulic conductivity and the triangular distribution function was observed and the best fit was achieved at an iteration number of 5000.However, increasing the number of iteration to 10 000 showed no significant difference to the 5000 iterations as shown in Figure7and compared to Figure6.As a result, a number of 5000 iterations was used to run all the simulations to estimate the likelihood of the travel time of the groundwater flow of each selected pair of wells.FiguresA1 through A8(Appendix A) display fit comparisons between the input data of the hydraulic conductivity and the triangular distribution function for the selected pairs of the wells which are listed in Table(2).

Figure 8 .Fit
Figure 8. Probability Distribution of Travel Time of Groundwater Flow between Well (RW-50) and Well (B7-13) with an Iteration Number Of 5000

Figure B7 .
Figure B7.Fit Comparison for Travel Time between Well (5579) and Well (5510)

Table 1 .
Selected Wells to Estimate Time of Travel for Groundwater in the Iraqi Western Desert

Well No. Ground surface above sea level (m) Static water level below ground surfaces (m) Hydraulic conductivity, K (m/day)
Basrawi (2007)were used to find the minimum, most likely, and maximum of the hydraulic conductivity in the area around and between wells (RW-5) and (B7-13).The other values of the minimum, most likely, and maximum of the hydraulic conductivity as well as the head difference, length (L), and hydraulic gradient (dh/dl) of each selected pair of wells are listed in Table(2) and were used to run the software of the Monte Carlo simulation to predict the travel time of the groundwater.Equation (1) which is used in Monte Carlo simulation to predict travel time of groundwater flow is expressed as follows:

Table 2 .
Hydraulic Parameters of the Selected Pairs of Wells Used in the Monte Carlo Simulation

Table 3 .
Summary of the Predicted Travel Time of Groundwater Flow for the Selected Pairs of Wells Used in the Monte Carlo Simulation Fit Comparison for Travel Time, T (day) between well (5006) and well (5621) .1,RiskShift(923750))Fit Comparison for Travel Time, T (day) between well (5616) and well(4910)