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.

This study presents a comprehensive analysis of the principles, application procedures, domestic and international case studies, and application constraints of double-packer-based groundwater sampling in deep fractured rock boreholes, with the aim of precisely identifying the hydrochemical characteristics of the candidate sites for the disposal of high-level radioactive waste. The hydrochemical characteristics of deep groundwater are crucial factors that must be identified to ensure the long-term safety of the disposal system. To obtain reliable data on these factors, both precise groundwater sampling techniques and clearly defined quantitative criteria for the application procedures must be established and implemented. The results of this study highlight the technical requirements that must be considered at each stage of the groundwater sampling process, which include drilling methods, preliminary planning before the sampling, removal of drilling-disturbed fluid, stabilization of water quality, and maintenance of closed-system sampling conditions. In addition, application cases from leading countries in geological disposal, such as Sweden, Switzerland, and Japan, were reviewed and compared with domestic cases to analyze differences in technological maturity and the level of procedural refinement. Furthermore, key technical limitations were identified, including short-circuiting caused by fracture networks, pumping constraints in ultra-low permeability zones, reduced mechanical stability in highly fractured zones, and equipment performance limitations under highpressure deep subsurface conditions. Based on these findings, potential countermeasures and supplementary tasks were proposed to address these challenges and enhance domestic technological capabilities.

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.

Traditional methods for geological characterization often overlook or oversimplify the challenge of subsurface non-stationarity. This study introduces an innovative methodology that uses ancillary data, such as geological insights and geophysical exploration, to accurately delineate the spatial distribution of subsurface petrophysical properties in large, non-stationary geological fields. The approach leverages geodesic distance on an embedded manifold, with the level-set curve linking observed geological structures to intrinsic non-stationarity. Critical parameters rho and beta were identified, influencing the strength and dependence of estimates on secondary data. Comparative evaluations demonstrated that this method outperforms traditional kriging, particularly in representing complex subsurface structures. This enhanced accuracy is crucial for applications such as contaminant remediation and underground repository design. While focused on two-dimensional models, future work should explore three-dimensional applications across diverse geological structures. This research provides novel strategies for estimating non-stationary geologic media, advancing subsurface characterization.

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.

This study presents analytical solutions for describing contaminant storage and release from an aquitard with linear source depletion (LSD) boundary conditions. We investigated three scenarios for trichloroethylene (TCE) mass exchange before and after the LSD period in an aquifer bounded by an adjacent aquitard based on the LSD dynamics, a resistance coefficient, and the aquitard thickness. The developed analytical solutions showed good agreement with measured profiles and breakthrough curves from a previous study. In three scenarios, the factors delaying the onset of TCE release into the aquifer were a decrease in the resistance coefficient, an increase in LSD period and aquitard thickness. The changes in the duration, mass, and rate of TCE storage in the aquitard during LSD loading process affected the equilibrium of the aquifer-aquitard concentration gradient. After TCE loading, the period maintained above the maximum contaminant level was directly related to the three variables; the longest plume persistence occurred when TCE penetration distance at transition point from storage to release coincided with the aquitard thickness. Overall, the developed analytical solution aids in evaluating the risk of plume persistence, enhancing site management efficiency, and reducing remediation costs.

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.

In the study, a 2-D numerical model-delineating the high-level radioactive waste repository surrounded by hosting granite-rock was designed to evaluate pressure build-up, elevated temperature, and fracture flow caused by heat generation at the repository. During the early stage (0.1 years), build-up pressure at the repository induced radial groundwater flow. On the other hand, buoyancy-driven vertical flow was dominant during the late stage (2,000 years). The pressure build-up and elevated temperature varied significantly in bentonite, excavation-damaged zone and granite, resulting in different groundwater velocity. Fracture flow in both single and intersecting fractures were also influenced by the pressure and temperature, but their impact varied depending on the geometrical properties (e.g. orientation, elevation, lateral position) of the fracture. In complex discrete fracture network models, build-up pressure in the early stage mainly generated preferential fracture flow through few pathways, while in the late stage, convection-circulating flow induced by the buoyancy was distinct. In the early stage, higher density and connectivity of fracture network did not lead to an increase in average fracture flow, whereas the circulating flow in the late stage became more frequent and larger. Especially, long fractures served as both preferential flow pathways in the early stage and conduits for circulating flow in the late stage. The circulating flow was highly dependent on the spatial distribution of the long fractures, causing a significant uncertainty of average fracture flow in the late stage. The low angle of intersection between fractures resulted in low connectivity that decelerated the average fracture flow in both stages.

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

<p>A new numerical method was developed to accurately and efficiently compute a solution of the nonlinear Richards equation with a layered soil. In the proposed method, the Kirchhoff integral transformation was applied. However, in the Kirchhoff integral transformation approach, the transformed Kirchhoff head has dyadic characteristics at the material interface between different soil types. To avoid the dyadic characteristics at the material interface, a truncated Taylor series expansion was applied to the Kirchhoff head at the material interface and so the Kirchhoff head was replaced with a single pressure head value at the material interface. Accordingly, through the Taylor series expansion, a set of algebraic equations in the one-dimensional control volume finite difference discretized system formed a tridiagonal matrix system. Through a series of numerical experiments, the new method was compared to other numerical methods to determine its superiority. The results clearly demonstrated that the approach was not only more computationally efficient, but also more accurate and robust than other numerical methods. Computational performance was greatly enhanced with the proposed method, and which could be used to simulate complicated heterogeneous flow at a large-scale watershed or regional scale.</p> <p><strong>Acknowledgments</strong></p> <p>This work was supported by the basic research project (23-3411) of the Korea Institute of Geoscience and Mineral Resources (KIGAM) funded by the Ministry of Science and ICT.</p>

A series of numerical simulations were performed to investigate the heterogeneity effect resulting from spatial distribution of clast-matrix in conglomerates on CO2 migration under the reservoir conditions. Natural conglomerate cores drilled from geological CO2 storage site in South Korea as well as constructed cores comprising regularly distributed clast-matrix were numerically analyzed. Throughout the study, CO2 distribution, velocity of the CO2 front, breakthrough curves (BTCs), differential pressure (ΔP), gravity number (Ngv), and capillary number (Nc) were assessed to highlight the heterogeneity effect. In the regularly distributed clast-matrix cores, the CO2 plume with the corrugated shape of the front migrated faster than the matrix-only core. Even though the clast ratio was changed more than 10 %, there was not much difference in the arrival velocity of the CO2 plume (∼2.0×10−4 m/min). Nevertheless, there was a correlation between the clast ratio and the arrival velocity, clearly shown in the L-5 core having alternating layers of matrix and clasts. At natural conglomerate cores (JC-1 and JC-2), the migration pattern of the CO2 plume significantly differed. However, the clast-matrix ratio and the arrival velocity were strongly correlated with a large variation in the arrival velocity (∼2.6×10−3 m/min) despite the small change of the clast-matrix ratio. The variance of the cumulative arrival velocity was correlated with the clast-matrix ratio and the irregularity of the BTCs, which implied that the instability in the frontal displacement of the CO2 plume became greater where the clast ratio increased more. The ΔP was significantly affected by the distribution of clasts and matrix; ΔP increased as the CO2 passed narrow or heterogeneous pathways enclosed by the clasts and decreased as CO2 passed wide or diverging pathways. Interestingly, the highest ΔP was shown in the cores having the smallest clast ratio, implied that the heterogeneous distribution of clast-matrix could influence a degree of overpressure generation. The lowest Ngv below 1 and highest Ngv appear in JC-2 where the total clast ratio is more than 50 % and in H-1 consisting of matrix-only, respectively. The obtained ranges of Nc show that all the models are under capillary-dominated flow conditions relevant for CO2 sequestration. The magnitude of both Ngv and Nc correlates to the total clast ratio.

The single-well push-pull (SWPP) test has been extensively used in parameter estimation models. Numerous analytical solutions to the problem are available, all of which consider solute transport in a non-uniform flow field from steady radial flow created by the SWPP test. However, the non-uniform flow velocity is based on the Thiem equation and varies in a spatial manner, rather than a temporal manner. No analytical solution has been previously described for the solute transport equation of SWPP under fully transient flow. In this study, the generalized integral transform technique (GITT) was used to develop a new semi-analytical solution for solute transport in the SWPP test under fully transient flow. Four phases of the SWPP test were included: injection, chaser, rest, and extraction. With the proposed solution, the differences between a transient flow SWPP model solution and a piecewise steady-state flow SWPP model solution were non-negligible; such differences increased with decreases in the dimensionless parameter related to aquifer flow properties, which is proportional to transmissivity, aquifer thickness, and porosity, but inversely proportional to the storage coefficient and pumping rate. Additionally, long tails of BTCs are characteristic of type curves under the transient flow model. The long tails of BTCs under transient flow could result in overestimation of dispersivity if the SWPP model uses piecewise steady-state flow, rather than transient flow, for parameter estimation. The proposed semi-analytical solution provides a useful tool for curve-fitting during parameter estimation when the effects of transient flow are significant. The proposed solution can also serve as a performance comparison for testing numerical solutions to the SWPP test.

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.

Surfactant flushing with intermittent air injection, referred to as enhanced flushing, has been proposed at a site in Korea contaminated by military activity to overcome the difficulty of treatment caused by a layered geological structure. In this study, we developed a simple numerical model for exploring the effects of various physical and chemical processes associated with enhanced flushing on pollutant removal efficiency and applied it in a field-scale test. This simple numerical model considers only enhanced hydraulic conductivity rather than all of the interacting parameters associated with the complex chemical and physical processes related to air and surfactant behavior during enhanced flushing treatment. In the numerical experiment, the removal efficiency of residual non-aqueous phase liquid (NAPL) was approximately 12% greater with enhanced, rather than conventional, flushing because the hydraulic conductivity of the low-permeability layer was enhanced 5-fold, thus accelerating surfactant transport in the low-permeability layer and facilitating enhanced dissolution of residual NAPL. To test whether the enhanced flushing method is superior to conventional flushing, as observed in the field-scale test, successive soil flushing operations were simulated using the newly developed model, and the results were compared to field data. Overall, the simulation results aligned well with the field data.

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.

The multiple sorptive sites on illitic clays (e.g., frayed edge [FES], type II [TS], and planar sites [PS]) are related to variable 137Cs immobilization in subsurface environments. This study investigated the diverse Cs+ sorption using 10 illitic clays under various competing K+ (distilled water-10-1 mol L-1) and Cs+ concentrations (10-7-10-3 mol L-1). In addition, the multisite cation exchange model was simulated the best-fit sorption models compared to experimental datasets, and subsequently, optimized the sorption capacities of multiple sorptive sites on the illitic clays. The best-fit sorption model exhibited that diverse Cs+ sorption of 10 illitic clays was closely linked to the individual sorption capacities at the FES (1.76 x 10-5 to 1.12 x 10-4 eq kg-1), TS (1.59 x 10-3 to 9.76 x 10-3 eq kg -1), and PS (2.14 x 10-2 to 1.51 x 10-1 eq kg -1). The FES predominantly sorbed Cs+ at low aqueous-phase concentrations, whereas the TS and PS contributed to Cs+ sorption at relatively high concentrations. These sorption capabilities were correlated to illite contents and crystallinity of 10 illitic clays so that such parameters could be significant factors to predict the Cs+ sorption for illitic clays. Finally, the 1-D transport simulations showed significantly diverse Cs+ retardation (Rd approximate to 20-100) at low dissolved Cs+, implying that the various FES capacities of illitic clays could play an important role on Cs+ distribution in actual radioactive contamination sites.