Retention and mobility of phosphogypsum constituents in carbonate aquifer rock materials

Abstract

The disposal of phosphogypsum, an acidic and metal-rich by-product of phosphate fertilizer production, represents a growing environmental challenge due to its potential to contaminate groundwater. Despite numerous reports of phosphogypsum leachate (PGL) impacting aquifers worldwide, the mechanisms governing the mobility and retention of contaminant elements remain poorly understood. To elucidate the underlying physical and chemical mechanisms governing the retention and mobility of the different components of the solution as it interacts with carbonate aquifer rock, we present a set of batch and column experiments using PGL and crushed aquifer rock. The results show that the flow of PGL through the rock induces non-uniform dissolution, creating preferential flow paths and reducing the rock–solution interactions and the retention of related elements. Moreover, rock dissolution causes a pH increase, affecting the speciation of the different elements and promoting precipitation. Element-specific retention was observed with some elements, e.g., Mo, Ge, Tl, and Rb, showing limited interaction and high mobility, raising concerns for groundwater contamination, especially for Tl, given its high toxicity. Other elements, including Al, Cr, B, Co, Ni, Cu, Zn, Cd, Cs, and U exhibited grain-size-dependent retention, with smaller grain sizes providing more surface area for sorption or precipitation. REEs were strongly immobilized across all conditions, indicating negligible mobility in acidic carbonate environments. The retention mechanisms of phosphogypsum-related elements include retention on the mineral surface and/or co-precipitation of newly formed minerals. Overall, these results underscore the need for site-specific assessments of PGL disposal in carbonate settings, given the distinct behavior of different elements and the dynamic nature of flow paths and pH conditions.

---

Environmental significance

Phosphogypsum leachates (PGL), generated by global phosphate fertilizer production, are acidic and metal-rich fluids that pose significant long-term threats to groundwater quality. This study provides detailed insight into the geochemical and transport processes governing the mobility of contaminants including lesser-studied elements such as thallium, germanium, molybdenum, and rare earth elements in carbonate aquifers. By examining interactions between PGL and porous subsurface media, the research reveals key mechanisms of sorption, reactivity, and flow path development that directly influence contaminant fate. These findings enhance our understanding of subsurface contamination processes and provide essential knowledge to inform remediation and risk assessment strategies near phosphogypsum waste sites globally, with implications for groundwater protection, ecological integrity, and public health.

1. Introduction

Phosphogypsum (PG), the primary byproduct of wet-phosphoric acid (WPA) production, results from manufacturing fertilizers from phosphate rock. The global annual production of PG reaches nearly 300 million tons.¹ A large proportion of this byproduct is stored in stacks, and estimates suggest that between 3 and 8 billion tons of PG piles are already distributed around the world.^2–9 Phosphogypsum is composed mainly of gypsum and residual phosphoric acid. However, it often contains impurities, including trace metals and radioactive elements derived from the phosphate rock feedstock. The principal constituents of phosphogypsum are fluoride, sulphate, calcium, phosphate, trace elements, and radionuclides. PG is typically formed as a slurry during the WPA process, and the process waters are recycled. This recycling of acidic process waters leads to the enrichment of trace elements and other potentially toxic substances within the PG. Moreover, using sulfuric acid during WPA retains an acidic character in the PG.

Large stockpiles of phosphogypsum, exposed to environmental weathering, generate phosphogypsum leachates (PGL). If left untreated, these leachates can contaminate surrounding water systems, contributing to the acidification and pollution of these waters. In dry stacking operations, rainwater and residual process waters can combine to form a hydraulic head, which increases the infiltration of acidic waters into underlying aquifers or nearby surface water bodies. This phenomenon presents significant environmental risks, as the contamination of groundwater systems and surrounding water bodies could have profound effects on both ecological systems and human health.^10,11 The long-term environmental consequences of these PG stacks are considerable, compromising the sustainability of groundwater, altering the properties of porous media, and potentially threatening biodiversity and public health.

The concentration of elements in phosphogypsum leachates is dictated mainly by the mineral composition of the phosphate ore processed.^12,13 Although trace elements are typically present in low quantities, their environmental toxicity can be severe. The concentration of these elements is influenced by the type of phosphate ore and the chemical processes, such as sulfuric acid treatment, used in phosphoric acid production. The release of PGL into the environment may occur due to the weathering of PG stacks or inadequate sealing of phosphogypsum ponds, which leads to leakage and infiltration into subsurface systems. This can result in significant ecological disturbances, with radionuclides and heavy metals accumulating in the food chain, posing long-term health risks.^3,4 Instances of groundwater contamination resulting from PG have been documented in numerous locations worldwide (Africa,¹⁴ Asia,¹⁵ North and South America,^16,17 and the Middle East¹⁸). For example, in Spain, Millán-Becerro et al.¹⁹ identified preferential flow paths of contaminants originating from phosphogypsum stacks with a discharge of high concentrations of pollutants (PO₄, F, U, As, and Cd). Small variations in pollutant concentrations were attributed to the interplay between dilution by rainwater recharge, which causes a slight decrease in concentration, and dissolution of the precipitated sulphate evaporitic salts, which caused an increase in the dissolved concentration in the leachates released from the PG. Outbakat et al.²⁰ reported that PG application alters soil physical properties (soil hydraulic properties, total porosity, and bulk density). Pliaka and Gaidajis⁵ indicated that the acidic PG leachate could accelerate water–rock interactions, releasing more metal(loid)s from the host rock into groundwater.

Groundwater contamination in various regions has been linked to PG stacks. In Tunisia, Melki and Gueddari²¹ observed that the Sfax-Agareb aquifer has been severely affected by PGL. Groundwater beneath PG stockpiles and downstream in the aquifer exhibited elevated concentrations of ortho-phosphate, fluoride, sulphate, acidity, and total dissolved solids compared to upstream water, with a marked increase in zinc and aluminium concentrations, indicating significant deterioration of water quality. In China, studies around phosphate production plants in Shifang City, Sichuan Province, revealed that dryland soils near PG piles contained elevated levels of heavy metals such as zinc, arsenic, copper, lead, cadmium, and mercury.²²

The mobility and retention of PGL-related elements, and, thus, the spatial distribution of the groundwater contamination, is governed by the chemistry of the PGL (e.g., pH, ionic strength, concentrations of different PG-related elements) and the composition of the aquifer rock.^7,13,19 The presence of geological carbonate formations, common aquifer rocks, underneath the ponds can facilitate relatively rapid infiltration of phosphogypsum into the groundwater system.^3,7,22–24 The interaction between the acidic PGL solution and carbonate rock leads to rock dissolution, which alters the rock permeability, porosity, and specific surface area. As a result of this non-uniform dissolution, preferential flow paths are formed, which reduce the rock–solution interaction and, thus, increase the mobility of the PG-related elements.²⁵

Groundwater is a crucial resource that sustains life, providing water to billions of people worldwide, supporting agricultural irrigation systems, and maintaining the health of various ecosystems.²⁶ Approximately 15.2% of the Earth's terrestrial surface comprises carbonate rocks like limestone and dolomite. Furthermore, about 9.2% of the global population relies on groundwater from aquifers within carbonate rock formations for their water supply.²⁷ This highlights the importance of carbonate aquifers in both sustaining human populations and supporting agricultural productivity. Therefore, prolonged infiltration of PGL from ponds situated on carbonate rocks is especially alarming, potentially resulting in significant regional groundwater contamination. An illustrative case of carbonate aquifer vulnerability was documented in Israel, where studies observed the spread of a PGL plume that contaminated a natural spring located more than 10 km from the contamination source.^18,28,29 This contamination resulted in the persistent elevation of PGL-related elements (e.g., Tl, U, Cd, and Ni) in the spring water, which severely impacted the local ecosystem and continued to increase in the spring water.^28–30

While extensive research has been conducted on the composition of the PG piles and the release of PG-related elements into the environment,^{11,23,31–35} studies that examine the mechanisms governing the mobility and retention of the various PG-related elements in aquifer material are generally scarce. No studies have specifically addressed these mechanisms in the context of carbonate rock aquifers. Given the growing body of evidence for groundwater contamination by PGL worldwide,^{14,15,17,21,30,36,37} it is crucial to investigate the mechanisms governing the mobility and retention of these elements in various environmental compartments. This understanding is critical to assessing and managing the long-term environmental risks of PGL release.

Therefore, the primary objective of this study was to elucidate the leading mechanisms governing the mobility and retention of PG-related elements as they interact with carbonate aquifer materials. We specifically focus this study on capturing interactions of various PGL elements with the aquifer rock, utilizing batch retention and column transport experiments with crushed carbonate rock of varying grain sizes. Additionally, we examine the changes in the hydraulic properties of the columns, which were influenced by the dissolution of the carbonate rock and the formation of CO₂ bubbles, thereby providing insights into the dynamic processes occurring in the aquifer.

2. Experimental

2.1. Sampling and sample handling

The phosphogypsum solution was collected using a bailer from a well (RPG-21) located near the PG ponds of the fertilizer plant in Mishor Rotem. The well was refreshed with three well volumes before the collection of the solution. The PG contamination was diluted ×3 using double distilled water (DDW; 18.2 > MΩ cm) to form a simulated PGL used in the batch and column experiments. The solution contained elevated concentrations of Na, Ca, Cl, P, and F (649, 817, 629, 2405, and 1797 mg L⁻¹, respectively). The concentrations of all elements present in the PGL are detailed in Table S1 in the SI. The solution pH, after dilution, was 3 (pH units). Zafit aquifer rock was sampled from an outcrop at coordinates 31.046810 N/35.083527 E. An XRD analysis of the Zafit rock showing its (dolomitic) mineral composition is provided at the SI (Fig. S1). The Zafit rock was crushed using a rock grinder and subsequently classified into three distinct grain-size fractions: 1–2 mm, 2–5 mm, and 5–10 mm. The crushed rocks were washed using DDW and dried at room temperature for seven days.

2.2. Batch adsorption experiments

Batch adsorption experiments were conducted by adding 10 g of rock (2–5 mm grains) to a flask containing 0.5 L PGL. The experiment was conducted in duplicate. The solution was allowed to be mixed with rock on a rotating table and sampled hourly/daily for 52 days. Control samples were taken before mixing the solution with the rock to determine the initial concentrations of the different elements. Solution sampling was conducted several seconds after manually shaking the experiment bottles, followed by sample filtration through a 0.22 µm filter. The sampled solutions were analysed for the different elements using inductively coupled plasma mass spectrometry (ICP-MS). At the end of the retention experiment, the solution was separated from the rock grains by filtration through a 1 mm mesh size net. Then, the powder of newly formed minerals was separated from the solution by filtration through 1.6 µm pore size microfiber filter paper (Whatman).

2.3. Column transport experiments

Two sets of column experiments were conducted. The first set examined the degree of interaction of the different elements with the aquifer rock using columns packed with aquifer rock with a mixture of grain sizes (1–2 mm, 2–5 mm, and 5–10 mm). The duration of each experiment was one week. The second set of experiments examined the changes in hydraulic properties that occurred in the column, using a red dye solution after flowing the PGL solution for one week.

2.3.1. Column experiments with different grain sizes. This set comprised three different column experiments, each with a different grain size. For each experiment, two polycarbonate columns, 19 cm long and 3 cm inner diameter, were packed with crushed rock from the Zafit formation and placed vertically. To eliminate gravity effects, the flow in the columns was from bottom to top via a multichannel peristaltic pump at a fixed flow rate of 0.6 ± 0.02 mL min⁻¹. The columns underwent a saturation phase by flowing DDW for 24 h before the experiment at a flow rate of 0.1 ± 0.01 mL min⁻¹. Then, the DDW was replaced with the PGL, and effluent collection at the column outlet began. After flowing the PGL for 60 h (32 ± 3 pore volumes (PV)), the solution was replaced with DDW for an additional 12 h (7 PV). The PV was determined by weighing the unsaturated column before the experiment started and after the column was fully saturated. Because the density of the aqueous solution is ∼1.00 g mL⁻¹ the PV can be directly calculated from the mass differences. The column porosity was calculated considering the PV and the total volume of the column derived from its dimensions. The porosity values for the three grain sizes were as follows: 1–2 mm PV = 66 mL and a porosity of 0.49, 2–5 mm PV = 71 mL and a porosity of 0.53, and 5–10 mm PV = 78 mL and a porosity of 0.58. Fig. S2, SI, shows the columns packed with grains of the Zafit aquifer rock in three sizes. At the end of the experiment, one of the replicates was used to examine the influence of PGL flow on the hydraulic properties of the column using a red dye solution (Section 2.3.2.). The other replicate was unpacked, and grains from the column inlet were further analysed by aqua regia digestion (Section 2.4.).

2.3.2. Red dye column experiment. To examine the changes in the flow regime of the column due to the interaction with the PGL, we compared the flow of red dye solution through the column that underwent PGL flow to a control column that experienced only DDW flow. Following the 24 h saturation of the columns with DDW, a red dye solution flowed through the columns for 8 h at a flow rate of 0.4 ± 0.02 mL min⁻¹. The red dye solution eluted from the column was collected at the column outlet, and the flow of the red dye through the column was monitored by imaging. The concentration of red dye in the effluents was measured using absorption spectroscopy at λ = 430 nm (UV-1600; Shimadzu Corp.).

2.4. Aqua regia digestion and sequential extraction

Aqua regia digestion procedures suitable for high-carbonate samples³⁸ were conducted to study the composition of the newly formed mineral powder extracted from the batch experiment and the retained minerals on the rock grains obtained from the inlet of the column. The retained mineral fraction was calculated as the difference between the elemental concentration measured in inlet-end column grains and the bulk concentration of the aquifer rock (Table S2). This inlet-derived retained fraction was then extrapolated to the full column mass (188 g) to obtain the solid retention efficiency, η, defined as the percentage of the input mass captured. To distinguish between surface sorption and mineral incorporation, and to quantify the accumulation of REEs in specific mineral phases, a sequential extraction procedure was conducted. This was performed on both the pristine bulk aquifer rock and the rock grains retrieved from the column inlet after the experiment. A modification of the procedure by Scheplitz et al. (2021), suitable for high-carbonate samples, was utilized to isolate the following phases: (1) adsorbed/exchangeable; (2) carbonates; (3) hydrous iron oxides; (4) crystalline iron oxides; and (5) sulfides/organic material (see Section S1 SI for details).

Table S3 summarize and compare the experimental parameters of the batch and column experiments.

2.5. ICP-MS

All samples were analysed via ICP-MS (Agilent 7700s) for the different element concentrations. Drift corrections were carried out using indium as an internal standard and by repeatedly analysing a calibration solution of 10 µg L⁻¹ concentration as a drift monitor throughout the analysis. Memory effects were avoided by additional manual cleaning using 5% HNO₃. The following masses were measured: Al27, P³¹, Ca40, Ca⁴⁴, Cr52, Fe⁵⁶, Co⁵⁹, Ni⁶⁰, Cu⁶³, Zn⁶⁶, Ge⁷¹, Br⁷⁹, Rb⁸⁵, Y⁸⁹, Mo⁹⁵, Cd¹¹¹, Cs¹³³, Ba¹³⁷, La¹³⁹, Ce¹⁴⁰, Ce¹⁴², Pr¹⁴¹, Nd¹⁴⁶, Nd¹⁵⁰, Sm¹⁴⁷, Sm¹⁴⁹, Eu¹⁵¹, Eu¹⁵³, Gd¹⁵⁵, and Gd¹⁵⁷, Tb¹⁵⁹, Dy¹⁶³, Dy¹⁶⁴, Ho¹⁶⁵, Er¹⁶⁶, Er¹⁷⁰, Tm¹⁶⁹, Yb¹⁷², Yb¹⁷⁴, Lu¹⁷⁵, Lu¹⁷⁶, Tl²⁰³, Tl²⁰⁵, and U²³⁸. All calibration curves were with R² > 0.98.

2.6. Visual MINTEQ

Speciation analysis and saturation indices of the various elements present in the PGL were performed using Visual MINTEQ (v3.1) at 25 °C, with the MINTEQ3.1 thermodynamic database and default complexation constants for all aqueous species. Model inputs included measured 3× diluted PGL composition (major ions and 39 trace elements from Table S1) with ionic strength 0.08–0.15 M (Davies equation) and pH range 3.0–5.7 matching experimental evolution; precipitation/dissolution was allowed for carbonates (dolomite, calcite), Ca-phosphates (hydroxyapatite, CaHPO₄), fluorides (CaF₂, MgF₂), REE minerals (REPO₄, REF₃), and metal hydroxides (Al, Cr, Cu, Zn, Cd, Ni) consistent with Zafit rock composition and observed secondary precipitates. Full input parameters and validation against experimental observations are detailed in Table S4 (SI). This analysis enabled the prediction of the species in which the elements exist under given environmental conditions, including pH, ionic strength, and the presence of natural organic matter. The analysis further identified the distribution of elements among free ions, complex forms, and potential precipitates, which are critical for evaluating their mobility and environmental behaviour.³⁹

3. Result and discussion

The interaction of the PGL while flowing through carbonate aquifer rock was examined in the context of the changes in the hydraulic properties of the matrix and the retention of the different elements present in the PGL. In general, as the acidic PGL interacts with carbonate rock, the resulting rock dissolution modifies the matrix's hydraulic properties, raises the solution pH, and drives the leaching of elements from the rock.^21,40,41 In parallel, the PGL-related elements are retained on the surface of the rock, and newly formed minerals are precipitated.^25,29 The changes in the retention of the PGL-related elements and changes in the hydraulic properties of the matrix are discussed in Sections 3.1. and 3.2., respectively.

3.1. Retention of PGL-related elements while interacting with carbonate rock

3.1.1. Batch experiment. The batch retention experiment, which lasted 52 days, reveals notable variations in the retention patterns of different elements (Fig. 1). Based on these observations, the elements can be classified into two distinct groups characterized by their unique retention behaviours.

Fig. 1 Batch retention experiment results (52 days duration) showing concentration profiles (C/C₀) for representative elements and solution pH changes. Panels (A–E) show Group I elements with limited retention (Ge, Ba, Rb, Mo, Ni); panels (F–J) show Group II elements with progressive retention (Al, Cr, B, Cu, Tl); panels (K–O) show Group III elements with rapid complete retention (representative REE shown: La); panel (P) shows solution pH evolution driven by carbonate rock dissolution. Each data point represents the average of two experimental replicates, with error bars indicating ±1 standard deviation.

The first group of elements is elements whose relative concentrations in solution (C/C₀) after 52 days are >1. This group includes Ba, Ni, Ge, Rb, and Mo (Fig. 1A–E). This limited interaction could result from the higher affinity of the elements toward the solution rather than the rock due to complexation with different ligands. The dominant ligands in this PGL solution are F⁻, Cl⁻, SO₄²⁻, and PO₄³⁻, which are present at concentrations of 1797, 629, 267, and 7535 mg L⁻¹, respectively (Table S1). The CO₃²⁻ could also serve as an important ligand, as its concentration in the solution changes according to the extent of the rock dissolution and the change in solution pH. For this batch experiment, the pH changed from 3 at the beginning of the experiment to 5.3 after 10 d and 5.7 after 52 d (Fig. 1P). This increase is attributed to the dissolution of the carbonate aquifer rock, which leads to the addition of carbonate ions to the solution (eqn (1)).

CaMg(CO₃)₂(s) + 4H⁺ → Ca²⁺ + Mg²⁺ + 2H₂O + 2CO₂(g) (1)

The second group comprises elements whose solution concentration after 52 days is lower than at the beginning of the experiment, with C/C₀ < 1. This group includes Al, B, Cr, Cu, Zn, Cd, Cs, Tl, U, and REE (represented by La; Fig. 1D–O). Within the second group, the retention of some elements, such as U, Al, and Cr is more robust, with a C/C₀ < 0.2 reached after several days. On the other hand, elements such as B and Tl exhibit more gradual retention. The REE group, represented in Fig. 1 by ‘La’, showed immediate and complete retention.

In the batch experiment, the emergence of newly formed minerals was detected visually, manifested as a fine powder settled at the bottom of the experimental bottles. This powder of newly precipitated minerals was extracted from the solution 52 days after the experiment and subsequently dissolved using aqua regia. The elemental composition that resulted from the dissolution of the newly formed mineral powder is presented in Fig. S3, SI.

High concentrations of P, Ca, and Mg, together with lower concentrations of B, Cd, Cr, Cu, Ni, U, and Zn were observed. Significant concentrations of P, Ca, and Mg are indicative of phosphorous mineral precipitation (e.g., hydroxyapatite formation, eqn (2)).

10Ca²⁺ + 6PO₄³⁻ + 2H₂O → Ca₁₀(PO₄)₆(OH)₂(s) + 12H⁺ (2)

The retention of Group II elements (Al, Cr, Cu, etc.) is pH-dependent (e.g., eqn (3) and (4), explaining the decline in solution concentration as pH increases from 3 to 5.7 (Fig. 1).

M²⁺ + 2H₂O ↔ M(OH)₂(s) + 2H⁺; K_sp = [M²⁺][OH⁻]² (3)

where M²⁺ represents transition metals (Zn²⁺, Cd²⁺, Cu²⁺, Ni²⁺).

U⁶⁺ + 4H₂O → UO₂(OH)₂(s) + 6H⁺ at pH > 5 (4)

Ahmad et al.⁴² reported a decrease in heavy metal mobility through adsorption and pH-driven precipitation. As CaCO₃ dissolves, it raises the solution pH, prompting the formation of stable carbonate or hydroxide phases that immobilize metals such as Cd, Zn, and Pb. Chromium's fate hinges on its oxidation state, with Cr³⁺ adsorbing or precipitating more readily than Cr⁶⁺, which must first be reduced. At higher pH, Cu²⁺ may form hydroxides or carbonates, and Fe³⁺ typically hydrolyses into Fe-hydroxides that interact with CaCO₃ surfaces.

The saturation indexes of potential minerals were calculated using Visual MINTEQ 3.1 software in the relevant pH range (pH 3–5). The saturation indexes of Ca and Mg minerals, the only minerals that displayed saturation indexes exceeding 1, are presented in Fig. S4, and Table S5, SI. The figure highlights the increasing saturation indices of these minerals as the solution pH rises, suggesting the possible precipitation of Ca–P minerals (CaHPO₄, CaHPO₄: 2H₂O, and hydroxyapatite), Ca–F minerals (fluorite; eqn (5)), and Mg–F mineral (MgF₂; eqn (6)).

Ca²⁺ + 2F⁻ → CaF₂(s) (5)

Mg²⁺ + 2F⁻ → MgF₂(s) (6)

The precipitation of fluorinated minerals commences at pH 3, whereas the precipitation of phosphorous minerals occurs with an increase in solution pH. The precipitation of these minerals could lead to the co-precipitation of different elements in the mineral lattice, as previously reported.⁴³ The detection of trace metals in the dissolved fraction suggests their co-precipitation within the newly formed minerals, as the saturation indexes of these elements with phosphorous and fluorine (e.g., Cd₃(PO₄)₂ and CdF₂) were lower than 1 (Table S6, SI). The saturation indexes of REE minerals in solution and the concentration of REE in the dissolved fraction of the mineral powder extracted from the batch solution after 52 days of the experiment are shown in Fig. S5 and S6, SI, respectively. Fig. S5, SI, depicts the potential precipitation of REE-F₃ and a majority of REE-PO₄ minerals within the pH range of 3–5, showcasing an increase in saturation index for these minerals with higher solution pH.

3.1.2. Column experiments of compound. Three sets of column transport experiments were conducted to examine the mobility and retention of elements found in the PGL and better frame the mechanisms that are related to surface area interactions. These experiments used columns filled with crushed aquifer rock of varying grain sizes, namely 1–2 mm, 2–5 mm, and 5–10 mm. The representative elements were categorized into three groups based on their level of interaction with the rock, which was determined by analysing the recovery of each element in the respective sets of column experiments (Fig. 2). This classification offers valuable insights into the behaviour of the elements in relation to their recovery, which is determined by the fraction of the total mass of the element eluted from the column compared to the total mass that initially entered the column.

Fig. 2 Recoveries of selected elements present in the PGL solution following column transport experiments with different grain sizes. Error bars represent the range between two experimental replicates (shown as ±standard deviation).

The first group includes elements that exhibit limited interaction with the rock, with recoveries >60% across all three-grain sizes. The second group comprises elements whose recoveries are grain-size dependent. For this group of elements, the recoveries are higher as the grain size diameter is higher due to the smaller reactive surface area of the column. The third group consists of elements completely retained on the aquifer rock. While batch experiments (52 days) approached equilibrium retention, column duration (60 h ≈ 2.5 days contact) represents early transient transport; Group I high mobility (>60%, Fig. S7) persisted across timescales, but slower Group II sorption may increase over field-relevant multi-year residence times requiring site-specific verification.

3.1.3. Group I – limited interaction: Ge, Rb, Mo, and Tl. Ge, Rb, Mo, and Tl demonstrated comparable mobility while flowing through columns packed with carbonate rock of different grain sizes (breakthrough curves are shown in Fig. S7). The breakthrough curve (BTC) of Ge, Rb, and Mo, shown in Fig. S7A–C, displays a steep increase in relative concentration to C/C₀ = ∼0.8, followed by a plateau. The BTC of Tl also increases sharply to C/C₀ = ∼0.8, but then Tl relative concentration gradually increases to C/C₀ = ∼1.1. The occurrence of C/C₀ > 1 could be attributed to the leaching of Tl from the source rock, which has a Tl concentration of 0.03 mg kg⁻¹ (Table S2), due to rock dissolution or cation exchange. Similar mobility and retention in columns with different reactive surface areas could be attributed to the limited interaction of these elements with the rock grains, as elaborated in Section 3.1.1.

Zhao et al.⁴⁴ demonstrate how carbonate minerals within mine tailings (e.g., calcite and dolomite) buffer the pH of percolating fluids. Zhao et al. outline a “dissolution–transport–re-precipitation” cycle that governs how metals move through carbonate tailings. Initially, acidic conditions drive metal release via carbonate and sulphide dissolution. Over time, as the pH climbs and secondary mineral precipitates accumulate, metals become immobilized again, settling out of solution onto reactive surfaces or forming stable mineral phases.

This interplay shapes the spatial extent of metal migration and ultimately dictates contaminant concentrations in drainage waters downstream of the tailings. Zhao et al.⁴⁴ reported that alkali and post-transition metals like Tl and Rb often remain relatively mobile but can be removed when many sorption sites exist. Visual MINTEQ speciation analysis (Table S4) reveals that Mo in the acidic PGL exists predominantly as molybdate oxyanions (MoO₄²⁻), which form stable aqueous complexes resistant to sorption on carbonate minerals^45,46 at pH < 5. Germanium exists primarily as undissociated H₃GeO₃⁰ and H₂GeO₃⁻ species under the experimental pH conditions.⁴⁷

The high solubility of Mo oxyanions and Ge silicic acid analogues, combined with their weak affinity for carbonate surfaces (compared to Group II transition metals), explains their rapid breakthrough and high recoveries in column experiments⁴⁷ (Fig. S7). The negligible grain-size dependence of Mo and Ge recovery (>60% across all grain sizes) confirms that element retention on carbonate rock is not the controlling mechanism for their environmental transport.^45–48

3.1.4. Group II – grain-size dependent behaviour: Al, Cr, B, Co, Ni, Cu, Zn, Cd, Cs, and U. The second group of elements consists of those whose recoveries were influenced by the grain size of the crushed rock used to pack the different columns (Fig. 3). Variations in grain size affect the reactive surface area of the columns and subsequently govern the extent of interaction between the rock and the solution. In columns with similar volumes, those packed with smaller grains exhibit a larger reactive surface area than those packed with larger grains. As a result, elements that interact with the rock experience increased retention processes when flowing through columns packed with smaller grains. Measuring retained elements in the column inlet confirms Group II grain-dependence: retention efficiency η declines by 20–96% from fine to coarse grains (e.g., Ni 47 → 31%, U 26 → 1%, Al 8 → 3%; Fig. S9). Finer grains (higher surface area) capture more despite heterogeneity; coarser grains show bypass via preferential paths (Fig. 4). Quantifies how non-uniform dissolution reduces overall retention efficiency in coarser media, as inferred from recoveries (Fig. 2). The recoveries of all selected elements in this group of elements increased as the reactive surface area of the column decreased (Fig. 3). Moreover, the mechanisms that control the mobility of each element differ and depend on its chemical properties (e.g., solubility, concentration, and speciation).

Fig. 3 Breakthrough curve (BTC) pattern classifications (P1–P5) of various elements as a function of grain size. Five distinct BTC patterns are identified based on concentration profiles: P1 (A), rapid rise followed by gradual decline correlated with pH increase; P2 (B), rapid rise to plateau without pH correlation; P3 (C), gradual increase independent of pH; P4 (D), rise-plateau-decline-plateau pattern inversely correlated to pH; P5 (E), complete retention with near-zero effluent concentration. These patterns indicate different retention mechanisms (sorption, precipitation, cation exchange) and their dependence on pH, surface area (grain size), and element speciation. Black symbols, 5–10 mm grains; Green symbols, 2–5 mm grains; Red symbols, 1–2 mm grains. A replicate of (F–O) is provided as Fig. S8.

Fig. 4 Red dye tracer experiment results showing the effect of phosphogypsum leachate (PGL) exposure on column flow paths. (A) Breakthrough curves (BTC) of red dye in two columns: PGL column (red line) was pre-exposed to PGL solution for 60 h creating preferential flow paths through selective carbonate dissolution and mineral precipitation; DDW control column (black line) received only deionized water flow. The PGL column shows: (1) faster red dye breakthrough (earlier C/C₀ rise), (2) lower peak concentration (indicative of broader tracer distribution), and (3) longer tailing (slower drainage from heterogeneous flow regions), confirming development of preferential flow paths. (B–E) Spatial distribution snapshots of red dye solution in the PGL column at four time intervals marked on the BTC curve in panel (A), showing non-uniform dye penetration and preferential pathways compared to expected symmetric flow in homogeneous media.

By examining the BTC patterns of the various elements in columns with different grain sizes, it is possible to compare the leading processes that govern transport and retention. It is noted that analysing a BTC requires caution, as the element BTC represents its concentration at the column outlet throughout the experiment, and processes occurring within the column cannot always be detected. The concentration of each element in the effluent solution can change due to the magnitude of the retention on the rock, the dissolution of the rock, the cation exchange between the rock and solution, and the release of previously retained elements. These processes depend on the hydraulic and chemical properties of the column (i.e., the solution pH, the reactive surface area, and the pore flow speed), which change as the PGL solution flows through the column (Section 3.2.).

The different BTCs were divided into five general patterns, reflecting five processes (P1–P5; Fig. 3A–E). The first pattern, referred to as ‘P1’, is characterized by a rapid rise in relative concentrations during the early stages of the BTC, followed by a gradual decline (Fig. 3A). This BTC pattern was measured for Al, Cr, B, Co, and Ni while flowing through columns packed with the largest grains (Fig. 3F–J, respectively; black symbols). This BTC shape exhibits an inverse trend to the solution pH, indicating a decrease in the mobility of these elements as the solution pH increases. This pattern can be attributed to two processes associated with pH: (1) a change in the speciation of the elements as a function of pH, which alters the mobility of the new species, and (2) reduced rock dissolution throughout the experiment, which leads to a decreased contribution of the element from the rock and, in turn, reduced concentrations of the element in the solution. During the first PVs, when the pH is relatively low, the rock dissolves, and the element concentrations in the effluents increase. Later, as the contribution of rock dissolution gradually decreases, the concentrations of the different elements also decrease.

The second observed BTC pattern, referred to as ‘P2’, is characterized by a rapid increase in concentration from a background level, followed by a plateau where the concentration stabilizes (Fig. 3B). This plateau could occur after a delay in elution (Cr, Fig. 3G; green symbols), immediately after the solution enters the Column (B, Fig. 3H; green and red symbols), or after a gradual increase in relative concentrations (Co and Ni, Fig. 3I and J, respectively; green symbols). These variations suggest that the BTC plateau pattern does not exhibit a direct response to changes in solution pH. The occurrence of a plateau BTC pattern for a particular element signifies that the element concentrations in the column have reached a state of saturation, resulting in the formation of mineral precipitates.

The third observed BTC pattern, referred to as ‘P3’, is characterized by a gradual increase in relative concentrations as the solution flows through the column (Fig. 3C). The gradual increase in concentration may commence either upon the solution entering the column (e.g., Cu, as shown in Fig. 3K with green symbols) or be delayed (e.g., Cs, as shown in Fig. 3N with red symbols). The elution variations suggest that the mobility of these elements is not linked directly to the solution pH, similar to the ‘P2’ pattern. The gradual increase in relative concentrations is attributed to the occupation of the retention sites of the rock that are specific to the retention of the specific element species present in the solution.

The fourth observed BTC pattern, referred to as ‘P4’, is characterized by an elution pattern exhibiting an initial increase in relative concentrations, followed by a plateau phase, and subsequently, a decrease in relative concentrations leading to another plateau (Fig. 3D). This BTC pattern was observed for Co, Ni, and Cd in columns packed with the smallest grain size (Fig. 3I, J, and M, respectively; red symbols). This pattern shows the opposite trend from the pH pattern. The inverse shapes of the BTC and pH curves imply that the mobility of these elements depends on the changes in the solution pH, similar to the ‘P1’ pattern. The influence of pH variation on the mobility of elements is similar to the mechanisms described for the ‘P1’ BTC pattern.

The fifth observed BTC pattern, referred to as ‘P5’, displays a flat BTC representing complete retention of the elements in the column (Fig. 3E). This BTC pattern was observed only for columns packed with the smallest grains for Al, Cr, Cu, Zn, and U (Fig. 3F, G, K, L, and O, respectively, red symbols). Zhao et al.⁴⁴ emphasize that each heavy metal exhibits distinct solubility and speciation patterns as the carbonate tailings evolve geochemically. Under acidic microenvironments, carbonate minerals dissolve quickly, releasing metals that were either sorbed or structurally incorporated. As pH increases, buffering by carbonate minerals, and as secondary minerals (e.g., iron oxyhydroxides) form, most metals become less mobile. Zinc and cadmium follow similar pathways, first dissolving under acidic conditions and then re-precipitating or adsorbing at higher pH. Copper exhibits a strong affinity for sorption onto iron oxyhydroxides, while nickel can remain soluble longer due to stable aqueous complexes but ultimately partitions into secondary mineral phases. Aluminium re-precipitates as Al-hydroxides as pH climbs, limiting its migration, while boron commonly remains soluble unless extreme alkalinity fosters adsorption. Column experiments were conducted in open-system conditions with continuous atmospheric contact at the inlet and outlet, maintaining aerobic conditions throughout the 7 days duration. Although redox potentials (Eh) were not directly measured, the open-system design and the availability of atmospheric oxygen suggest predominantly oxidizing conditions during the experiments. For redox-sensitive elements (Cr, U), speciation in the PGL was determined by Visual MINTEQ assuming oxidizing conditions (pe ≈ 13 at pH 3–5), predicting Cr⁶⁺ as the dominant oxidized Cr species and U⁶⁺ as the stable uranium species. The observed grain-size-dependent retention of Cr and complete retention of U in the smallest grain-size columns (Fig. 3, P1 and P5 patterns) may reflect both the chemical properties of these oxidized species and their interaction with the carbonate rock. In future studies, controlled redox experiments (reduced vs. oxidized conditions) would help isolate redox effects on element mobility and retention.^49,50

3.1.5. Group III – complete retention: rare earth elements. Complete retention of the REE group was observed across all three column experiments, regardless of the grain size employed (data not shown). This complete retention of the REEs demonstrated that either the complexation of the REE was insufficient to overcome the binding affinity of REEs to the surface sites of the carbonate rock and/or REEs precipitated in the column.⁵¹

An aqua regia digestion of the aquifer rock was conducted to investigate the REE concentration before and after the column experiment. The total concentration of each REE in the aquifer rock before the retention of REEs in the column experiment, measured using aqua regia digestion, was 0.17–5.5 mg kg⁻¹ (Table S2). After completion of the column experiments, the column was disassembled, and analysis of the inlet samples revealed REE concentrations of 57–73 and 49–60 mg kg⁻¹ for the two replicates, demonstrating consistent and robust retention of REEs on the rock surface (Fig. S10). In addition, the presence of REEs was also detected in the powder of new precipitated minerals obtained from the batch experiment (Section 3.1.1.; Fig. S2).

To elucidate the specific retention mechanisms, the fractionation of REEs in post-experiment column grains was compared with that in pristine aquifer rock using sequential extraction (Fig. S11). The pristine rock exhibited low background REE concentrations (0.17–5.5 mg kg⁻¹) distributed naturally across the mineral phases. In contrast, the rock grains from the column inlet showed a marked enrichment in REE concentrations, specifically within the hydrous iron oxide (Step 3) and carbonate (Step 2) fractions. Significant concentrations were also associated with crystalline iron oxides (Step 4). Conversely, the concentration of REEs in the exchangeable/adsorbed fraction (Step 1) remained negligible (<1 µg L⁻¹) in both pre- and post-experiment samples.

This comparison provides quantitative evidence that REE retention is not governed by reversible surface adsorption (Step 1). Instead, the dominant mechanisms are scavenging by iron oxides—consistent with the significant iron content of the Zafit rock (approx. 9500 mg kg⁻¹) -and coprecipitation into the carbonate lattice. This aligns with the ‘dissolution-transport-reprecipitation’ cycle, in which REEs are immobilized in secondary mineral phases formed during the neutralization of the acidic PGL. Robust REE retention while flowing through carbonate rock and carbonate-rich soils was previously reported e.g.,^28,52,53 Amiel et al.,⁵¹ showed that the retention of REE increases as the carbonate content of the matrix increases.

3.2. Effect of rock dissolution on the changes in hydraulic properties

The flow of the PGL solution, characterized by acidic pH and high metal content (Table S1), caused changes in the hydraulic properties of the column due to dissolution and precipitation reactions. In general, carbonate dissolution caused by acidification expanded the pore volume of the carbonate rock, while precipitation of newly formed minerals, resulting from the super-saturation of chemical constituents, filled the pore space.²⁵ The interplay between these contrasting processes leads to the development of preferential flow paths.^25,54

The formation of preferential flow paths in the column experiment was examined by simultaneously flowing red dye solution through two columns packed with 2–5 mm carbonate rock. The first column experienced the flow of PGL solution before the red dye experiment (‘PGL’ column). The second was a control column, which experienced the flow of DDW (‘DDW’ column). Fig. 4 shows the BTC of the red dye in these columns. The red dye solution in the ‘PGL’ column emerged faster, reached lower C/C₀ values in the ‘peak’, and had a longer ‘tail’ than the control column. Fig. 4B–E shows the distribution of the red dye solution in the columns in four time-steps throughout the experiment. These figures show that the red dye solution progresses faster in the ‘PGL’ column due to forming preferential flow paths, relative to the more homogeneous flow in the control column (‘DDW’). The higher mobility of the red dye in the ‘PGL’ column, compared to the ‘DDW’ column, implies that the dissolution and precipitation processes occurring in the ‘PGL’ column increase not only the porosity of the rock but also its permeability.^54,55

The formation of a preferential flow path, resulting from the blockage of pores by minerals, reduces the reactive pore volume of the column, increasing the pore flow velocity and rendering the solution–rock interaction less efficient.^25,56 These processes reduce both the retention of elements on the rock and the dissolution of the rock itself. In parallel, and because of the greater heterogeneity of the PGL column, a small portion of the solution flows more slowly through the less mobile regions, creating a longer ‘tail’ in the ‘PGL’ column BTC than the ‘DDW’ column, resulting in higher solution–rock interaction efficiency.⁵⁴ Inlet solid retention (η, Fig. S9) quantifies heterogeneity: coarser grains retain ∼50% less PGL elements (Fig. S9; e.g., Ni 47 → 31%). Preferential paths (Fig. 4) bypass reactive zones, reducing contact despite some stagnant retention (longer tails). Overall efficiency drops, amplifying mobility risks.

3.3. Implications for field-scale applications

Laboratory experiments using crushed rock provide controlled conditions ideal for mechanistic understanding, but differ from natural in situ carbonate aquifers. First, the increased surface area of crushed material (three discrete grain-size fractions) enhances element retention compared to intact rock containing fractures and solution-enlarged pathways. Second, preferential flow paths in natural aquifers are typically localized within fracture networks and dissolution channels, creating greater heterogeneity than observed in our uniformly packed columns. Third, natural carbonate aquifers host diverse biological communities (bacteria and archaea) and organic matter that can influence element speciation and precipitation, factors not considered in this sterile system. Fourth, the pH-induced dissolution patterns in our study assume acidic conditions persist; in natural systems, pH buffering by surrounding limestone may occur more rapidly. Therefore, the observed element mobilization represents a maximum mobility scenario under continuous acid exposure, with enhanced surface area. Field predictions should account for aquifer mineralogy (limestone vs. dolomite), natural pH-buffering capacity, groundwater quality, and flow-regime characteristics.^25,57,58

4. Conclusions

This study employed batch and column experiments to investigate the mobility and retention of various elements derived from PGL ponds when interacting with carbonate aquifer rock. Through these experiments, we examined the underlying physical and chemical processes that govern the behaviour of PGL elements. The passage of PGL solution through the columns resulted in rock dissolution and mineral precipitation, which subsequently influenced the hydraulic characteristics of the systems. The dissolution of the carbonate rock by the acidic PGL solution simultaneously increased the solution pH and caused the blockage of pores by mineral precipitation. The formation of pore blockages induced the emergence of preferential flow paths, which were observed using a red dye solution in the column experiment. The preferential flow led to a decrease in the available reactive surface area. This, in turn, reduces the interaction of the solution with the rock, which results in reduced leaching.

The different elements in the PGL solution showed high variation in their retention behaviour. Some elements (Mo, Ge, Tl, and Rb) exhibited limited interaction with the rock, resulting in high recoveries in the column and batch experiments and similar BTC shapes and recoveries as the reactive surface area of the solution changed. These elements were not (or were minimally) affected by changes in the solution pH. The mobility of these elements is expected to be relatively high in cases where low-pH industrial effluent infiltrates through a groundwater system composed of carbonate rock. The notable mobility of Tl in these cases raises concerns regarding its potential impact on human health, considering the substantial toxicity associated with this element e.g.,.⁵⁰ Another group of elements showed grain size dependence behaviour, resulting in decreased recoveries as the grain size was smaller, yielding higher reactive surface area. The mechanisms that led to this grain size-dependent behaviour vary within this group of elements. The mechanism of each element in each experiment was analysed by the correlations of its BTC shape with the physical and chemical processes occurring in the column. For example, some retention patterns correlated to the change in solution pH, while others were more gradual, independent of the solution acidity. The retention of the REE group exhibited the highest level of robustness compared to all other elements in the solution. These elements were entirely retained in all column experiments. Thus, REEs are not expected to be mobile in carbonate aquifer rock exposed to acidic solutions. Overall, these results underscore the need for site-specific assessments of PGL disposal in carbonate settings, given the distinct behaviour of different elements and the dynamic nature of flow paths and pH conditions.

Author contributions

Nitai Amiel: conceptualization, methodology, investigation, validation, visualization, writing – original draft; Han Zhang: investigation, validation, data curation; Ishai Dror: resources, validation, supervision, funding acquisition, writing – review & editing; Brian Berkowitz: resources, supervision, funding acquisition, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: the modified sequential extraction procedure (Section S1), elemental concentrations in the experimental phosphogypsum leachate (PGL) and bulk aquifer rock composition (Tables S1 and S2), and detailed batch and column experimental conditions (Table S3). It further contains Visual MINTEQ inputs, parameters, and calculated saturation indices for PGL-related elements and mineral phases over pH 3–5 (Tables S4–S6; Fig. S4–S5), XRD confirmation of Zafit aquifer rock mineralogy (Fig. S1), images of column setups for different grain sizes (Fig. S2), concentrations of elements in newly precipitated minerals from batch experiments (Fig. S3 and S6), breakthrough curves for elements exhibiting limited rock interaction across grain sizes (Fig. S7), solid-phase retention efficiency calculations and trends (Section S4; Fig. S9), and measured elemental concentrations on aquifer rock before and after PGL column exposure, including sequential extraction results documenting REE fractionation (Section S5; Fig. S10 and S11). See DOI: https://doi.org/10.1039/d5em00254k.

Acknowledgements

The authors gratefully acknowledge the support of a research grant from the Crystal Family Foundation, Mr Stephen Gross, the Emerald Foundation, the P. & A. Guggenheim-Ascarelli Foundation, and the DeWoskin/Roskin Foundation.

References

1. M. Contreras, R. Pérez-López, M. J. Gázquez, V. Morales-Flórez, A. Santos, L. Esquivias and J. P. Bolívar, Fractionation and fluxes of metals and radionuclides during the recycling process of phosphogypsum wastes applied to mineral CO₂ sequestration, Waste Manage., 2015, 45, 412–419.
2. E. Bilal, et al., Phosphogypsum circular economy considerations: A critical review from more than 65 storage sites worldwide, J. Cleaner Prod., 2023, 414, 137561.
3. Y. Chernysh, O. Yakhnenko, V. Chubur and H. Roubík, Phosphogypsum recycling: A review of environmental issues, current trends, and prospects, Appl. Sci., 2021, 11, 1575.
4. L. Hermann, F. Kraus and R. Hermann, Phosphorus processing—potentials for higher efficiency, Sustainability, 2018, 10, 1482.
5. M. Pliaka and G. Gaidajis, Potential uses of phosphogypsum: a review, J. Environ. Sci. Health, Part A, 2022, 57, 746–763.
6. R. El Zrelli, L. Rabaoui, H. Abda, N. Daghbouj, R. Pérez-López, S. Castet, T. Aigouy, N. Bejaoui and P. Courjault-Radé, Characterization of the role of phosphogypsum foam in the transport of metals and radionuclides in the Southern Mediterranean Sea, J. Hazard. Mater., 2019, 363, 258–267.
7. R. Pérez-López, F. Macías, C. R. Cánovas, A. M. Sarmiento and S. M. Pérez-Moreno, Pollutant flows from a phosphogypsum disposal area to an estuarine environment: an insight from geochemical signatures, Sci. Total Environ., 2016, 553, 42–51.
8. R. Pérez-López, J. M. Nieto and J. D. de la Rosa, Environmental tracers for elucidating the weathering process in a phosphogypsum disposal site: implications for restoration, J. Hydrol., 2015, 529, 1313–1323.
9. G. Zhou, Z. Zheng, R. Liu, C. Liu and R. Tao, Multi-stage precipitation for the eco-friendly treatment of phosphogypsum leachates using hybrid alkaline reagents, J. Water Process Eng., 2023, 53, 103626.
10. M. S. Al-Hwaiti and J. F. Ranville, Distribution of potentially toxic metal and radionuclide contamination in soils related to phosphogypsum waste stockpiling in the Eshidiya Mine, Jordan, Geochem.: Explor. Environ. Anal., 2010, 10, 419–433.
11. R. Pérez-López, J. M. Nieto, I. López-Coto, J. L. Aguado, J. P. Bolívar and M. Santisteban, Dynamics of contaminants in phosphogypsum of the fertilizer industry of Huelva (SW Spain): from phosphate rock ore to the environment, Appl. Geochem., 2010, 25, 705–715.
12. F. Macías, J. M. Nieto and A. M. Sarmiento, Trace elements from phosphate rock in east-central Europe: Environmental impact assessment, Environ. Sci. Technol., 2017, 51, 3021–3027.
13. P. M. Rutherford, M. J. Dudas and R. A. Samek, Environmental impacts of phosphogypsum, Sci. Total Environ., 1994, 149, 1–38.
14. F. Telahigue, B. Agoubi, F. Souid and A. Kharroubi, Groundwater chemistry and radon-222 distribution in Jerba Island, Tunisia, J. Environ. Radioactiv., 2018, 182, 74–84.
15. S. J. P. Cañete, L. J. H. Palad, E. B. Enriquez, T. Y. Garcia and T. Yulo-Nazarea, Leachable 226Ra in Philippine phosphogypsum and its implication in groundwater contamination in Isabel, Leyte, Philippines, Environ. Monit. Assess., 2008, 142, 337–344.
16. A. Miller and J. W. Sutcliffe Jr, Effects of phosphogypsum on soil properties and plant growth, Soil Sci., 1982, 133, 278–283.
17. A. J. G. Santos, P. S. C. Silva, B. P. Mazzilli and D. I. T. Fávaro, Radiological characterisation of disposed phosphogypsum in Brazil: evaluation of the occupational exposure and environmental impact, Radiat. Prot. Dosimetry, 2006, 121, 179–185.
18. A. Burg and I. Gavrieli, Groundwater contamination and water-rock interaction during leakage of industrial wastewater into a carbonate aquifer in an arid zone, Israel, Procedia Earth Planet. Sci., 2013, 7, 101–104.
19. R. Millán-Becerro, R. Pérez-López, C. R. Cánovas, F. Macías and R. León, Phosphogypsum weathering and implications for pollutant discharge into an estuary, J. Hydrol., 2023, 617, 128943.
20. M. B. Outbakat, K. El Mejahed, M. El Gharous, K. El Omari and A. Beniaich, Effect of phosphogypsum on soil physical properties in Moroccan salt-affected soils, Sustainability, 2022, 14, 13087.
21. S. Melki and M. Gueddari, Impact of phosphogypsum leachate on groundwater quality in Sfax-Agareb, Tunisia, J. Environ. Manage., 2018, 217, 614–622.
22. L. F. O. Silva, M. L. S. Oliveira and T. J. Crissien, A review on the environmental impact of phosphogypsum and potential health impacts through the release of nanoparticles, Chemosphere, 2022, 286, 131513.
23. H. Tayibi, M. Choura, F. A. López, F. J. Alguacil and A. López-Delgado, Environmental impact and management of phosphogypsum, J. Environ. Manage., 2009, 90, 2377–2386.
24. E. Saadaoui, N. Ghazel, C. Ben Romdhane and N. Massoudi, Phosphogypsum: potential uses and problems – a review, Int. J. Environ. Stud., 2017, 74, 558–567.
25. X. Zhang, J. Guo, Q. Hu, X. Gao, C. Li, M. Luo and Y. Wang, Effects of Fe-rich acid mine drainage on percolation features and pore structure in carbonate rocks, J. Hydrol., 2020, 591, 125571.
26. T. Gleeson, Y. Wada, M. F. P. Bierkens and L. P. H. van Beek, Water balance of global aquifers revealed by groundwater footprint, Nature, 2012, 488, 197–200.
27. I. Vadillo and L. Ojeda, Carbonate aquifers threatened by contamination of hazardous anthropic activities: challenges, Curr. Opin. Environ. Sci. Health, 2022, 26, 100336.
28. N. Amiel, I. Dror, A. Zurieli, Y. Livshitz, G. Reshef and B. Berkowitz, Selected technology-critical elements as indicators of anthropogenic groundwater contamination, Environ. Pollut., 2021, 284, 117156.
29. A. Burg and J. Guttman, Mitigation of downstream propagation of contaminated water in a carbonate aquifer – the northeastern Negev desert, Israel, Sci. Total Environ., 2019, 654, 550–562.
30. A. Burg, Ein Boqeq – Chemical Composition after Fifteen Years of Contamination, Geological Survey of Israel, Report GSI/03/2011, pp. 4–7.
31. J. M. Abril, R. García-Tenorio, S. M. Enamorado, M. D. Hurtado, L. Andreu and A. Delgado, The cumulative effect of three decades of phosphogypsum amendments in reclaimed marsh soils from SW Spain: ²²⁶Ra, ²³⁸U, and Cd contents in soils and tomato fruit, Sci. Total Environ., 2008, 403, 80–88.
32. M. Al-Hwaiti, V. Carney, J. F. Ranville, and P. E. Ross, Heavy metal assessment of phosphogypsum waste stockpile material from Jordan, in Proc. 22^nd American Society of Mining and Reclamation Annual National Conference, 2005, pp. 19–23.
33. M. S. Al-Masri, Y. Amin, S. Ibrahim and F. Al-Bich, Distribution of some trace metals in Syrian phosphogypsum, Appl. Geochem., 2004, 19, 747–753.
34. J. M. Arocena, P. M. Rutherford and M. J. Dudas, Heterogeneous distribution of trace elements and fluorine in phosphogypsum by-product, Sci. Total Environ., 1995, 162, 149–160.
35. L. A. Attar, K. Shamali, B. A. Ghany and S. Kanakri, Case study: heavy metals and fluoride contents in the materials of Syrian phosphate industry and in the vicinity of phosphogypsum piles, Environ. Technol., 2012, 33, 143–152.
36. A. May and J. W. Sweeney, Assessment of environmental impacts associated with phosphogypsum in Florida (USA), Rep. Invest. U.S. Bur. Mines, No RI 8693, U.S. Dept. Int., Washington, 1982.
37. E. V. Molski and M. N. Kotchneva, Groundwater contamination in Voskresiensk industrial region (Russia), in Future Groundwater Resources at Risk, 1994, pp. 169–175.
38. J. Scheplitz, S. Koopmann and H. Fröllje, et al., Optimization and assessment of a sequential extraction procedure for calcium carbonate rocks, Environ. Monit. Assess., 2021, 193, 577.
39. J. P. Gustafsson, 'Visual MINTEQ 3.0 user guide', KTH, Department of Land and Water Resources, Stockholm, Sweden, 2011, p. 550.
40. H. Coetzee, M. Kotoane, M. Atanasova and F. Roelofse, 'Interactions between dolomite and acid mine drainage in the Witwatersrand—Results of field and laboratory studies and the implications for natural attenuation in the West Rand Goldfield', in Reliable Mine Water Technology, IMWA, Golden, 2013, pp. 307–312.
41. G. Kaufmann, F. Gabrovšek and D. Romanov, Dissolution and precipitation of fractures in soluble rock, Hydrol. Earth Syst. Sci., 2016, 1–30.
42. K. Ahmad, I. A. Bhatti, M. Muneer, M. Iqbal and Z. Iqbal, Removal of heavy metals (Zn, Cr, Pb, Cd, Cu, and Fe) in aqueous media by calcium carbonate as an adsorbent, Int. J. Chem. Biochem. Sci., 2012, 2, 48–53.
43. K. E. Wright, T. Hartmann and Y. Fujita, Inducing mineral precipitation in groundwater by addition of phosphate, Geochem. Trans., 2011, 12, 1–13.
44. R. Zhao, X. Wu, G. Zhu, X. Zhang, F. Liu and W. Mu, Revealing the release and migration mechanism of heavy metals in typical carbonate tailings, East China, J. Hazard. Mater., 2024, 464, 132978.
45. S. Goldberg, Mechanism of molybdenum adsorption on soils and soil minerals, in Soil Mineralogy with Environmental Applications, Soil Science Society of America, 2009, pp. 345–372.
46. Y. Arai, Molybdate adsorption by birnessite: effects of pH and ionic strength, Appl. Clay Sci., 2015, 115, 1–8.
47. S. Koopmann, H. Prommer, A. Siade and T. Pichler, Molybdenum mobility during managed aquifer recharge in dolomitic aquifers, Environ. Sci. Technol., 2023, 57, 7298–7308.
48. V. A. Pokrovskii and H. C. Helgeson, Thermodynamic properties of aqueous germanium species to 1000 °C and 3 kbar, Geochim. Cosmochim. Acta, 1998, 62, 3751–3766.
49. D. Rai, B. M. Sass, and R. C. Moore, Remediation Processes of Hexavalent Chromium from Groundwater: A Review, IWA Publishing, 2004.
50. M. E. Romero-González, L. E. Fletcher and C. J. Williams, Uranium(VI) attenuation in a carbonate-bearing oxic alluvial aquifer, Chem. Geol., 2019, 524, 1–12.
51. N. Amiel, I. Dror and B. Berkowitz, Mobility of rare earth elements in coastal aquifer materials under fresh and brackish water conditions, ACS Environ. Au., 2024, 4, 186–195.
52. A. Brewer, I. Dror and B. Berkowitz, Electronic waste as a source of rare earth element pollution: Leaching, transport in porous media, and the effects of nanoparticles, Chemosphere, 2022, 287, 132217.
53. K. P. Saripalli, P. D. Meyer, D. H. Bacon and V. L. Freedman, Changes in hydrologic properties of aquifer media due to chemical reactions: a review, Crit. Rev. Environ. Sci. Technol., 2001, 31, 311–349.
54. C. N. Fredd and H. S. Fogler, Influence of transport and reaction on wormhole formation in porous media, fluid mechanics and transport phenomena, AIChE J., 1998, 44, 1933–1949.
55. C. Noiriel, L. Luquot, B. Madé, P. Raimbault, P. Gouze and J. van Der Lee, Changes in reactive surface area during limestone dissolution: an experimental and modelling study, Chem. Geol., 2009, 265, 160–170.
56. L. Osorio-Rico, A. Santamaria, and S. Galván-Arzate, Thallium toxicity: General issues, neurological symptoms, and neurotoxic mechanisms, in Neurotoxicity of Metals, ed. M. Aschner and L. Costa, Adv. Neurobiol., Springer, Cham, 2017, vol. 18, ch. 17.
57. P. G. Cook, A Guide to Regional Groundwater Flow in Fractured Rock Aquifers, Geoscience Australia, 2nd edn, 2003.
58. G. Medici, L. J. West, P. J. Chapman and S. A. Banwart, Prediction of contaminant transport in fractured carbonate aquifer types: a case study of the Permian Manesian Limestone Group (NE England, UK), Environ. Sci. Pollut. Res., 2019, 26(24), 24863–24884.

---