Taiwan's mid-20th-century Blackfoot disease epidemic highlighted the public health implications of geogenic arsenic in groundwater and prompted extensive water-supply interventions. Despite these measures, arsenic persists in aquifers, continuing to affect drinking water, irrigation, and aquaculture. This study developed a scalable machine learning model to identify areas across Taiwan at risk of groundwater arsenic concentrations exceeding the WHO guideline of 10 µg/L. The model demonstrated robust performance (mean AUC 0.93, balanced accuracy 0.87). A detailed evaluation and interpretation of the principal predictor variables governing arsenic dissolution provides insights into geochemical and hydrogeological controls on mobilization. Furthermore, the resulting hazard map of Taiwan was used to estimate the human population at risk as well as the locations and types of exposed agriculture and aquaculture. Hotspots are clustered in the Chianan, Pingtung and Lanyang plains, where many rural communities still rely on untreated groundwater. Population exposure is substantial, with an estimated 148,000 people at risk in 2023, which is about half the 303,000 reported in 1998. Agricultural impacts are also considerable: >80 % of aquaculture areas, and 33 % of paddy rice fields occur within high-risk zones. Several high-risk zones lie outside Taiwan's Groundwater Control Areas, revealing regulatory blind spots. By connecting Taiwan's historical arsenic crisis to novel predictive tools, this work supports evidence-based strategies to reduce current exposures and prevent future public arsenic-related health impacts.

ABSTRACT Multiple successive tracer tests are often conducted to obtain reliable breakthrough curve results under regional groundwater flow, especially when the accuracy is crucial. In such cases, the period of rest between the end of the first divergent tracer test and the initiation of the second divergent tracer test allows the tracer from the first test to travel along with the background regional flow, thereby influencing the distribution of residual tracer concentration. This residual tracer could potentially interfere with breakthrough curve results from the tracer injection in the second tracer test. Additionally, the conventional analytical solution used for the divergent tracer test considers only radial flow; regional flow and consecutive tracer tests are ignored. Consequently, interpreting the behaviour of the tracer in consecutive divergent tracer tests under regional flow conditions is challenging using conventional measures because of background regional concentration. This study proposes a new semi‐analytical solution, considering the effects of divergent radial and regional flows in consecutive tracer tests, addressing a critical gap in the conventional analytical solutions that, despite their practical necessity, have not been previously developed. The proposed semi‐analytical solution was subjected to parameter studies under various scenarios. In our case studies, the conventional analytical solution based on a single tracer test can be safely used for parameter estimation only in cases where the injected mass for the subsequent tracer test is approximately six‐fold that of the first tracer test or if the drift time is longer than 10 days.

Chlorinated solvents, common groundwater contaminants, can cause coexistence of the original contaminant and its degradation products during the transport process. Practically applicable analytical models for reactive transport are essential for simulating the plume migration of chlorinated solvent contaminants and their degradation products within a complex chemical mixture. Although several analytical models have been developed to solve advection–dispersion equations coupled with a series of decay reactions for simulating transport of the coexisting chlorinated solvent contaminants, the majority assume static, time-invariant inlet boundary conditions. Such time-invariant inlet boundary conditions may fail to adequately represent the temporal evolution of dissolved source discharge concentration concerning mass reduction, especially in the context of diverse DNAPL source remediation strategies. This study seeks to derive analytical models for three-dimensional reactive transport of multiple contaminants, specifically addressing the challenges posed by dynamical, time-varying inlet boundary conditions. The model development incorporates two distinct inlet functions: exponentially decaying and piecewise constant. Analytical solutions are obtained using three integral transform techniques. The accuracy of the newly developed analytical models is verified by comparing them with solutions derived from existing literature using multiple illustrative examples. By incorporating two distinct time-varying inlet boundary conditions, the models exhibit strong capabilities in capturing the complex transport dynamics and fate of contaminants within groundwater systems. These features make the models valuable tools for improving the understanding of subsurface contaminant behavior and for quantitatively evaluating and optimizing a range of remediation strategies.

A permeable reactive barrier (PRB) containing zero-valent iron (ZVI) is an in situ groundwater remediation technology that passively intercepts and treats contaminated groundwater plumes. Over time, secondary mineral precipitation within the PRB diminishes porosity and hydraulic conductivity, altering flow paths, residence times, and sometimes causing bypass of the reactive zone. This study utilizes the THMC software to simulate porosity reduction in a PRB, capturing the coupled effects of fluid flow and geochemical interactions. The simulation results indicate that porosity loss is most significant at the PRB entrance and stabilizes beyond 0.2 m. Porosity reduction is primarily caused by aragonite, siderite, and ferrous hydroxide precipitating in pore spaces. The model further elucidates the influence of groundwater chemistry, demonstrating that variations in bicarbonate concentrations significantly impact mineral precipitation processes, thereby leading to porosity reduction. Furthermore, the study highlights reaction kinetics, with anaerobic iron corrosion rates being critical in controlling porosity reduction via mineral precipitation. THMC software effectively simulates porosity reduction in PRBs, identifies key factors driving clogging, and informs design optimization for long-term remediation.

Most existing multispecies transport analytical models primarily focus on inlet boundary sources, limiting their applicability in real-world contaminated sites where contaminants often arise from multiple internal sources. This study presents a novel semi-analytical model for simulating multispecies contaminant transport driven by multiple time-dependent internal sources. The model incorporates key transport mechanisms, including advection, dispersion, rate-limited sorption, and first-order degradation. In particular, the inclusion of rate-limited sorption addresses limitations in traditional equilibrium-based models, which often underestimate pollutant concentrations for degradable species. The derivation of this semi-analytical model utilizes the Laplace transform, finite cosine Fourier transform, generalized integral transform, and a sequence of inverse transformations. Results indicate that the concentrations of contaminants and their degradation products are highly sensitive to the variations in time-dependent sources. The model’s most significant contribution lies in its capability to simulate the contaminant transport from multiple internal pollution sources at a contaminated site under the influence of rate-limited sorption. By enabling the representation of multiple time-varying sources, this model fills a critical gap in analytical approaches and provides a necessary tool for accurately assessing contaminant transport in complex, realistic pollution scenarios.

Chlorinated solvents can degrade to generate transformation products sequentially. The presence of such transformation products must be considered for the health risk assessment. Recently, we have developed a software package MUSt (MUltiSpecies transport Analytical Models) in which FORTRAN executable files based on our newly developed multispecies transport analytical solutions are equipped with an interactive graphical user interface (GUI). The multispecies transport analytical solutions embedded in MUSt have been further combined with health risk module for more reasonable health risk assessment. This study assesses the health risk for chlorinated solvent contaminated groundwater in northern Taiwan. The geographical distribution of non-carcinogenic and carcinogenic health risk is depicted for appropriate action for reducing groundwater concentration level of chlorinated solvent contaminants and protecting human health.

Traditional numerical models have been widely employed to simulate the transport of multispecies reactive contaminants in groundwater systems; however, their high computational cost limits their applicability in real-time or large-scale scenarios. Recent advances in artificial intelligence (AI) offer promising alternatives, particularly data-driven machine learning techniques, for accelerating such simulations. This study presents the development of a surrogate model based on artificial neural networks (ANNs) to simulate the transport and decay of interacting multispecies contaminants in groundwater. High-fidelity training datasets are generated through finite difference-based reactive transport simulations across a wide range of environmental and geochemical conditions. The ANN model is trained to learn the complex nonlinear relationships governing the multispecies transport and transformation processes. Model validation reveals that the ANN surrogate accurately reproduces the spatial–temporal concentration profiles of both original and degradation species, capturing key dynamic behaviors with high precision. Notably, the ANN model achieves up to a 100-fold reduction in computational time compared to traditional analytical or semi-analytical solutions. These results highlight the ANN’s potential as an efficient and accurate surrogate modeling tool for groundwater contamination assessment, offering a valuable advancement for decision-making in environmental risk analysis and remediation planning.

Most of the semi-analytical and analytical models employed to depict multidimensional, multispecies transport of sequentially degrading reaction products are built upon solving a set of coupled advection -dispersion equations (ADEs). Within these equations, sorption is considered as being equilibrium-controlled. However, it has been demonstrated that more realistic predictions of the transport of contaminants in the groundwater could be obtained by the use of a rate-limited sorption process instead of an equilibrium sorption assumption. This study is thus designed to devise semi-analytical models for the two-dimensional multispecies transport of a chemical mixture comprised of a parent compound and its degradation-daughter products which is influenced by ratelimited sorption subject to arbitrary time-dependent inlet boundary conditions. Three integral transforms are applied to generate a set of linear algebraic equations (AEs) resulting from the reductions of the ADEs. The contaminant concentration of each species is calculated by solving these AEs and then retransforming the solutions back to the original time-space domain. The simulation results obtained with our newly developed semianalytical model are nearly identical to those generated using a numerical model based on the Laplace transform finite difference (LTFD) method. This confirms the validity of the new semi-analytical model. An investigation of the effect of rate-limited sorption on the migration of the contaminant plume is carried out using various timedependent inlet boundary conditions. The sorption rate values vary from low to high, specifically, 0.05, 0.5, 5, and 50 years  1. The results show that the predicted concentrations of all contaminants within in the decay chain decrease as the sorption rate constant increases, for both constant and exponentially time-dependent boundary sources. However, under pulse loading boundary conditions, the concentrations of the later degradation products tend to increase with the sorption rate constant. The semi-analytical models created in this study, allow for different inlet boundary conditions, so can be utilized to simulate the transport of sequentially degrading reaction products. This greater accuracy makes them useful for many applications.

Multispecies transport analytical models that solve advection-dispersion equations (ADEs) are efficient tools for evaluating the transport of decaying contaminants and their sequential products. This study develops a novel semi-analytical model to simulate the multispecies transport of decaying contaminants, considering nonequilibrium sorption and decay in both dissolved and sorbed phases. First-order reversible kinetic sorption equations with decay processes are coupled to ADEs. Recursive analytical solutions, using the Laplace transform and generalized integral transform, are developed to address the mathematical complexity of the governing equations. The model's simulation results show excellent agreement with both numerical models and existing analytical solutions. Applied to a four-member radionuclide decay chain, the model reveals that including decay in the sorbed phase results in a lower concentration of the first member and avoids underestimating the radioactivity concentrations of daughter elements. These differences in dissolved radioactivity concentrations between models with and without sorbed phase decay may impact health risk assessments for radioactive waste disposal. Finally, this study provides a more sophisticated mathematical tool for analyzing multispecies transport in real field conditions where nonequilibrium sorption processes predominantly occur.

Download This Paper Open PDF in Browser Add Paper to My Library Share: Permalink Using these links will ensure access to this page indefinitely Copy URL Copy DOI

In recent years, climate change has disrupted the hydrological cycle, intensifying water management challenges around the world, leading to overabundance on some areas and scarcity in others. Groundwater as a dependable source of water has become a crucial part of water resource management. Unfortunately, in some places, natural arsenic (As) contamination in the groundwater threatens the safety of the user and has harmed human health. The prediction of groundwater As by using machine learning algorithms to improve efficiency is a critical issue. However, most of the existing machine learning models rely on geological and chemical monitoring data but fail to consider redox-dependent As variation induced by pumping and rainfall. Besides, comprehensive time series data for predicting temporal variations in As concentrations were rarely addressed by prior studies. The objectives of this study are to evaluate the trends and correlations between time series data for groundwater As concentrations, natural precipitation, and electrical power consumption of irrigation wells by dynamic factor analysis (DFA), and to develop an artificial neural network (ANN) for the prediction model of groundwater As variation. The results indicate that rainfall-induced recharge during the wet season significantly increases As concentrations, while groundwater pumping activities, particularly power consumption, are key factors influencing variations across seasons. The radius of influence of pumping wells plays a crucial role in altering redox conditions and As cycling in the groundwater. An ANN model, utilizing precipitation and pumping well power consumption, demonstrates high predictive accuracy with mean absolute error (MAE) values ranging from 0.001 to 0.012, mean squared error (MSE) values from 3×10-6 to 3×10-4, correlation coefficients (COR) from 0.23 to 0.85, and coefficients of determination (R²) from 0.96 to 0.99. These quantitative results underscore the model's effectiveness in reducing the cost of water quality analysis and mitigating the risk of As pollution in the food chain, providing valuable insights for water resource management and public health protection in As-affected regions.

This study aims at making a comprehensive assessment of the impact of land use and the hydrogeological properties on groundwater quality. First, factor analysis (FA) is applied to reveal the main pollutant sources and hydrogeological processes controlling the groundwater quality. FA identifies the four most important factors. Factor 1 (seawater salinization) is characterized by a medium loading of land use type of aquaculture. It is recognized that the high scores for factor 1 in coastal areas are due to over-pumping from aquafarms. Focused land use management is required to prevent saline-water intrusion in coastal aquifers. Factor 3 (nitrate pollution) shows high correlations with the land use type of fruit farming and the gravel thickness in unsaturated layers. High scores for factor 3 are also found in the proximal area of the Chuoshui River Alluvial Fan and the northeastern mountain area in the Pingtung Plain. Fruit farmers should be educated to reduce the application of fertilizers and promote the organic fruit farming. The impacts of land use and the hydrogeological properties on both Factor 2 (arsenic enrichment) and Factor 4 (reductive dissolution of Fe2+ and Mn2+) are negligible. Second, cluster analysis (CA) is performed on computed scores of the four main factors to separates 123 monitoring wells into cluster 1 (low polluted zone), cluster 2 (nitrate polluted zone) and cluster 3 (hybrid polluted zone). The results obtained from CA provide practical applications such as reduce agrichemical use in the areas of cluster 2 and enforce intensive monitoring in the prioritizing areas of cluster 3. This study successively uses the FA and CA to extract the meaningful information present by geographical visualization of scores for 4 main factors and 3 distinct clusters zones. The results are essential for formulating sound groundwater resource and land use management policies to ensure groundwater sustainability.

The contamination of groundwater aquifers by chlorinated solvents is a water-resource issue of great importance worldwide. Such contaminants are difficult to treat because they are often released as dense non-aqueous phase liquids (DNAPLs). Research shows that the application of remediation technologies to both the NAPL source and dissolved plume can lead to more efficient remediation, rather than to either alone. In these remediation efforts, analytical models that evaluate behavior and fate of contaminants do provide a better understanding of the performance of these remedial technologies. To the best of our knowledge, there exist no analytical model of simulating the plume migration of multiple contaminants with capabilities of accounting for both NAPL source and plume remediation simultaneously and different retardation for original chlorinated solvent contaminant and its degradation byproducts. In this study, we present a new analytical model for remediating both NAPL source and downgradient contaminant plume in groundwater at sites contaminated with chlorinated solvents and their degradation products with different retardation factors as well as considering both NAPL source and plume remediation simultaneously. A source model that accounts for the depletion of mass by the processes of dissolution or first-order decay reactions, corresponding with the removal or destruction of the source mass, is coupled to a plume reactive transport model. The source model is accounted for by relating source mass to the flux-averaged source discharge concentration through a power function. The developed analytical model considers 1-D advection, 3-D dispersion, first-order decay reactions and ingrowth as well as linear isothermal equilibrium sorption. The proposed analytical solution was derived through successive application of the Laplace transform in time and the double finite Fourier cosine transform regarding y and z. The correctness of the analytical model and its auxiliary FORTRAN computer program code are proved by showing excellent agreements between the simulated plume concentrations of all contaminants obtained from the derived analytical model and from a semi-analytical model available in the literature. Application of the proposed analytical solutions illustrates that the use of identical retardation factors for all contaminants may lead to underestimation or overestimation of the mobility of the contaminants, in cases when the retardation factors of the individual contaminants are greatly different from the identical retardation factor value adopted in all contaminants. From the experiments on six scenarios corresponding six remedial treatments, we found out that both the enhanced source decay and partial removal of source mass are main controlling factors at reducing the concentrations of all the contaminants, whereas plume decay leads to effective reduction in the concentrations of PCE, however, rather it causes unfavorable increases of the concentrations of the degradation byproducts. Ultimately, the developed model is used to better understand the impacts of various possible combinations of remedial efforts and management decisions on remediation of the subsurface contamination and quantify the benefit of a certain remediation decision.

Although efforts to analytically model multispecies reactive transport have already been reported, most of the studies utilized the advection-dispersion equations coupled with sequential decay chain reactions with a constant dispersion coefficient. Over the past 30 or 40 years, however, some studies have demonstrated that the dispersivity increases with the distance the solute travels, which is an outcome of variation in the hydraulic characteristics of the subsurface environment. The numerical modelling of multispecies reactive transport associated with a scale-dependent solute dispersion process has been discussed in the literature. In this study, an analytical model for plume migration of a chemical mixture comprised of the original pollutant and its degradation-related byproducts, subject to a scale-dependent dispersion process, is developed by taking advantage of the Laplace transform technique for the temporal variable and the generalized integral transform technique for the spatial variable. The correctness of the developed model is evaluated by comparisons with a numerical model in which the same set of governing equations are solved using the Laplace transform finite difference (LTFD) technique. The computational results obtained with the analytical and LTFD numerical models agree perfectly. To illustrate the impact of a scale-dependent dispersion process on the multispecies plume migration of the original contaminant and its degradation-related byproducts, this paper makes a comparison of the developed multispecies model with a model of the constant dispersion coefficient.

Stepwise reductive dechlorination of tetrachloroethene (PCE) to trichloroethene (TCE), three dichloroethylene (DCE) isomers (1,1-DCE, cis -1,2-DCE, and trans -1,2-DCE), vinyl chloride (VC), and ethane (ETH) may proceed under anaerobic conditions. However, most multispecies transport models for describing the plume migration of a chemical mixture comprising the original chlorinated solvent and its dechlorinated byproducts in the literature are often simplified to a sequential first-order reaction network which cannot account for the divergent reactions from TCE to the three DCE isomers and convergent reactions from the three DCE isomers to VC. In this study, general analytical solutions to multispecies transport equations with a complex reaction network were derived for a combination of semi-infinite and finite systems. The developed analytical solutions were robustly verified against a semi-analytical solution that can consider the same complex reaction network. The verification results indicated that the derived analytical solutions were accurate and robust. The general solutions derived in the present study were also used to investigate the effects of the outlet boundary conditions on solute transport involving a complex reaction network. The results showed that the analytical solution derived for infinite outlet BCs predicted lower concentrations of contaminants near the outlet boundary than those for finite outlet BCs. Moreover, a straight decay chain model may over-or under-estimate the DCE and VC con-centrations compared to the multiple branching isomer reaction model. The developed analytical solutions with a complex reaction network provide more realistic and efficient tools for assessing the movement of chlorinated solvents and their degradation-related byproducts in soil-water systems.

Groundwater resources are abundant and widely used in Taiwan's Lanyang Plain. However, in some places the groundwater arsenic (As) concentrations far exceed the World Health Organization's standards for drinking water quality. Measurements of the As concentrations in groundwater show considerable spatial variability, which means that the associated risk to human health would also vary from region to region. This study aims to adapt a back-propagation neural network (BPNN) method to carry out more reliable spatial mapping of the As concentrations in the groundwater for comparison with the geostatistical ordinary kriging (OK) method results. Cross validation is performed to evaluate the prediction performance by dividing the As monitoring data into three sets. The cross-validation results show that the average determination coefficients (R-2) for the As concentrations obtained with BPNN and OK are 0.55 and 0.49, whereas the average root mean square errors (RMSE) are 0.49 and 0.54, respectively. Given the better prediction performance of the BPNN, it is recommended as a more reliable tool for the spatial mapping of the groundwater As concentration. Subsequently, the As concentrations estimated obtained using the BPNN are applied to develop a spatial map illustrating the risk to human health associated with the ingestion of As-containing groundwater based on the noncarcinogenic hazard quotient (HQ) and carcinogenic target risk (TR) standards established by the U.S. Environmental Protection Agency. Such maps can be used to demarcate the areas where residents are at higher risk due to the ingestion of As-containing groundwater, and prioritize the areas where more intensive monitoring of groundwater quality is required. The spatial mapping of As concentrations from the BPNN was also used to demarcate the regions where the groundwater is suitable for farmland and fishponds based on the water quality standards for As for irrigation and aquaculture.

This study presents novel exact analytical solutions to a set of simultaneous three-dimensional advection-dispersion equations coupled with a sequential first-order degradation reaction network involving distinct retardation factor values among individual species. The analytical solutions to the coupled partial differential equation system subject to both the first - and third -type inlet source boundary conditions are derived by consecutive application of the three integral transformations in combination with sequential substitutions. The correctness of the developed analytical solutions is confirmed through numerical comparisons of three verification examples between our derived exact analytical solutions and a three-dimensional single-species analytical model in a semi-infinite domain, a two-dimensional multispecies exact analytical model in a finite domain and a three-dimensional multispecies semi-analytical model in a semi-infinite domain of the previous studies, respectively. The advantage of the derived analytical solutions is that its computational efficiency is 104 times the computational efficiency of two-dimensional multispecies analytical solutions for a finite domain when two solutions are simultaneously used to solve a two-dimensional multispecies transport problem. The developed analytical solutions are then used to evaluate the performance of the public domain BIOCHLOR model that simulates aquifer remediation by natural attenuation of dissolved multispecies at a chlorinated-solvent contaminated site provided by the Center for Subsurface Modeling Support (CSMoS) of USEPA. Results show that the BIOCHLOR model that was developed based on three-dimensional analytical model assuming a single retardation factor value for all dissolved species would not be suitable for simulating most multispecies plume migration at contaminated sites where each contaminant has its own retardation factor value. The effects of the inlet source boundary conditions on the multispecies plume migration are also investigated. The high computational efficiency of the developed analytical model in this current study renders it a very efficient tool for executing a probabilistic health risk assessment involving a large number of simulations required for problems of multispecies transport of dissolved chlorinated solvent in groundwater at a contaminated site with uncertainties in many different variables.

Lanyang Plain in northwestern Taiwan is an intensively productive agricultural area, most from the cultivation of crops and aquaculture. Groundwater fi

Although the average municipal water coverage in Taiwan is quite high, at 93.91%, only around half of the residents in the Pingtung Plain use tap water originating from the Taiwan Water Corporation to meet their needs. This means the exploitation of a substantial amount of groundwater as a source of water to meet drinking, agriculture, aquaculture, and industry requirements. Long-term groundwater quality surveys in Taiwan have revealed obvious contamination of the groundwater in several locations in the Pingtung Plain, with measured concentration levels of some groundwater quality parameters in excess of the permissible levels specified by the Taiwan Environmental Protection Administration. Clearly, establishing a sound plan for groundwater quality protection in this area is imperative for maximizing the protection of human health. The inappropriate use of hazardous chemicals and poor management of land use have allowed pollutants to permeate through unsaturated soil and ultimately reach the underlying shallow unconfined groundwater system. Thus, the quality of the water stored in shallow aquifers has been significantly affected by land use. This study is designed to characterize the relationship between groundwater quality and land use in the Pingtung Plain. This goal is achieved by the application of factor analysis to characterize the measured concentrations of 14 groundwater quality parameters sampled from 46 observation wells, the area percentages for nine land use categories in the neighborhood of these 46 observation wells, and the thicknesses of four unsaturated types of soil based on core samples obtained during the establishment of 46 observation wells. The results show that a four-factor model can explain 56% of the total variance. Factor 1 (seawater salinization), which includes the groundwater quality parameters of EC, SO42−, Cl−, Ca2+, Mg2+, Na+, and K+, shows a moderate correlation to land used for water conservation. Factor 2 (nitrate pollution), which includes the groundwater quality parameters of NO3−-N and HCO3−, shows a strong correlation to land used for fruit farming and a moderate correlation to the thickness of the gravel comprising unsaturated soil. Factor 3 (arsenic pollution), which is composed of groundwater quality parameters of total organic carbon (TOC) and As, is very weakly affected by land use. Factor 4 (reductive dissolution of Fe3+ and Mn2+), which involves Mn2+ and Fe3+, is weakly impacted by land use. Based on a geographic visualization of the scores for the four different factors and the patterns for land use, we can demarcate the areas where the groundwater in shallow unconfined aquifers is more vulnerable to being polluted by specific contaminants. We can then prioritize the areas where more intensive monitoring might be required, evaluate current land use practices, and adopt new measures to better prevent or control groundwater pollution.

Transport behaviors of contaminants through a heterogeneous formation consisting multiple layers are complicated because of the different physical and chemical properties for each individual layer. Few analytical solutions for single-species contaminant transport in a multi-layer heterogeneous formation have been reported in the literature. Some contaminants of concern such as radionuclide, nitrogen and chlorinated solvents can decay or degrade to form new successor products during their transport processes, thus making migration of these contaminants much complicated. Clearly, analytical models for multispecies transport coupled by a series of decay reactions in a multi-layer formation are useful tools for synchronous determination of the fate and transport of the predecessor and successor species of decaying or degradable contaminants. This study attempts to develop an analytical model for the multispecies reactive transport of degradable or decaying contaminants through a multi-layer heterogeneous formation. The derived analytical model is shown to be correct and accurate as the consistent results of comparisons between the derived analytical model and the numerical model. The developed analytical model will provide a more reliable predicting tool for real world application.