1. Introduction

2. Conceptual model and mathematical model

2.2. Governing equations

The governing equation involved in the model is as follows [26]:

(1) Water flow equation

The equation can be established from the conservation of mass:

$$\frac{\partial M_{\kappa }}{\partial t}=-\nabla F_{\kappa }+q_{\kappa }$$

(1)

For water phase:

$$F_w=x_{wl}{\rho }_l\overrightarrow{u_l}+x_{wg}{\rho }_g\overrightarrow{u_g},q_w=q_{wl}+q_{wg}$$

(2)

The mass flux of the multi-phase is according to Darcy's law:

$$\overrightarrow{u_{\beta }}=-k\frac{k_{r\beta }}{{\mu }_{\beta }}\left(\nabla P_{\beta }-{\rho }_{\beta }\overrightarrow{g}\right)$$

(3)

Combine with the mass conservation equation, the reactive transport model can be expressed as:

$$\frac{d\left(\phi S_l{c}_{jl}\right)}{dt}=-\nabla \left(c_{jl}\overrightarrow{u_l}-\tau \phi S_l{D}_l\nabla c_{jl}\right)+{\sum }_{i=1}^N{\nu }_{jn}{r}_n$$

(4)

Table 1. The mathematical model parameters.

Table 1. The mathematical model parameters.

symbol	meaning	symbol	meaning
${\mathit{M}}_{\mathit{\kappa }}$	The total mass of matter $\mathit{\kappa }$	$\mathit{S}$	saturation
$\mathit{t}$	Time	$\mathit{\phi }$	porosity
${\mathit{F}}_{\mathit{\kappa }}$	Flow rate of matter $\mathit{\kappa }$	$\mathit{k}$	Absolute permeability
${\mathit{q}}_{\mathit{\kappa }}$	Source and sink of matter $\mathit{\kappa }$	${\mathit{k}}_{\mathit{r}}$	Relative permeability
$\mathit{w}$	water	$\mathit{\mu }$	viscosity
$\mathit{c}$	concentration	$\mathit{l}$	Liquid phase
$\mathit{x}$	Mass fraction	$\mathit{g}$	Gas phase
$\mathit{\rho }$	density	$\mathit{D}$	Dispersion
$\mathit{u}$	Velocity of flow	$\mathit{r}$	reaction rate
$\mathit{\nu }$	stoichiometric number	$\mathit{j}$	species $\mathit{j}$
$\mathit{P}$	pressure	$\mathit{n}$	Reaction $\mathit{n}$

(2) Chemical reaction mathematical models

The main chemical reaction mathematical models in the simulation include equilibrium minerals and kinetic minerals. The equilibrium mineral model is:

$${\Omega }_m=K_m^{-1}\prod _{j=1}^{N_c}{c}_j^{{\nu }_{mj}}{\gamma }_j^{{\nu }_{mj}}$$

(5)

Where $K_m$ is the equilibrium constant.

The kinetic mineral model is:

$$r_n=\pm k_n{A}_n{\left|1-{\Omega }_n^{\theta }\right|}^{\eta }$$

(6)

Where $A_n$ is the specific surface area of the mineral; $k_n$ is the nth parallel mineral precipitation or dissolution reaction rate constant that depends on temperature;

Based on the Arrhenius equation, the correlation between $k_n$ and temperature can be written as:

$$k_n=k_{25}^{nu}exp\left[-\frac{E_a^{nu}}{R}\left(\frac{1}{T}-\frac{1}{298.15}\right)\right]+\sum _i\left\{k_{25}^{i}exp\left[-\frac{E_a^i}{R}\left(\frac{1}{T}-\frac{1}{298.15}\right)\right]\right\}\prod _j{a}_{ij}^{n_{ij}}$$

(7)

where $\prod _{\mathrm{j}}{a}_{ij}^{n_{ij}}$ describes the effect of specific ion activity on the i th parallel mineral precipitation and dissolution reaction;

(3) Models of mineral dissolution or precipitation

Based on the geological condition, the uranium mineral is UO₂. The minerals include Quartz, K-feldspar, Na-feldspar, oligoclase, Na-smectite, Ca-smectite, illite, kaolinite, gypsum and anhydrite, hematite, muscovite, dolomite, siderite, and ankerite. These minerals are also considered in the model. Here we've listed some chemical reactions.

$$2\mathrm{N}\mathrm{a}\mathrm{A}\mathrm{l}\mathrm{S}{\mathrm{i}}_3{\mathrm{O}}_8+{\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4+9\mathrm{H}2\mathrm{O}\to 2{\mathrm{N}\mathrm{a}}^{+}+{\mathrm{S}\mathrm{O}}_4^{2-}+{\mathrm{A}\mathrm{l}}_2{\mathrm{S}\mathrm{i}}_2{\mathrm{O}}_5{\left(\mathrm{O}\mathrm{H}\right)}_4+4{\mathrm{H}}_4{\mathrm{S}\mathrm{i}\mathrm{O}}_4$$

(8)

$$2\mathrm{K}\mathrm{A}\mathrm{l}\mathrm{S}{\mathrm{i}}_3{\mathrm{O}}_8+{\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4+9{\mathrm{H}}_2\mathrm{O}\to 2{\mathrm{K}}^{+}+{\mathrm{S}\mathrm{O}}_4^{2-}+{\mathrm{A}\mathrm{l}}_2{\mathrm{S}\mathrm{i}}_2{\mathrm{O}}_5{\left(\mathrm{O}\mathrm{H}\right)}_4+4{\mathrm{H}}_4{\mathrm{S}\mathrm{i}\mathrm{O}}_4$$

(9)

$$2{\mathrm{K}}_{0.5}{\mathrm{M}\mathrm{g}}_{0.25}{\mathrm{A}\mathrm{l}}_{2.3}\mathrm{S}{\mathrm{i}}_{3.5}{\mathrm{O}}_{10}{\left(\mathrm{O}\mathrm{H}\right)}_2+{\frac{11}{10}\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4+\frac{63}{20}{\mathrm{H}}_2\mathrm{O}\to \frac{1}{2}{\mathrm{K}}^{+}+{\mathrm{M}\mathrm{g}}^{2+}+{\mathrm{S}\mathrm{O}}_4^{2-}+{\mathrm{A}\mathrm{l}}_2{\mathrm{S}\mathrm{i}}_2{\mathrm{O}}_5{\left(\mathrm{O}\mathrm{H}\right)}_4+\frac{6}{5}{\mathrm{H}}_4{\mathrm{S}\mathrm{i}\mathrm{O}}_4$$

(10)

$$\mathrm{C}\mathrm{a}{\mathrm{C}\mathrm{O}}_3+{\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4\to \mathrm{C}\mathrm{a}{\mathrm{S}\mathrm{O}}_4+{\mathrm{C}\mathrm{O}}_2+{\mathrm{H}}_2\mathrm{O}$$

(11)

$$\mathrm{C}\mathrm{a}\mathrm{M}\mathrm{g}{\left({\mathrm{C}\mathrm{O}}_3\right)}_2+{2\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4\to \mathrm{C}\mathrm{a}{\mathrm{S}\mathrm{O}}_4+{\mathrm{M}\mathrm{g}}^{2+}+{\mathrm{S}\mathrm{O}}_4^{2-}+2{\mathrm{C}\mathrm{O}}_2+{2\mathrm{H}}_2\mathrm{O}$$

(12)

$${\mathrm{H}}_4{\mathrm{S}\mathrm{i}\mathrm{O}}_4\to {\mathrm{S}\mathrm{i}\mathrm{O}}_2+2{\mathrm{H}}_2\mathrm{O}$$

(13)

$${\mathrm{F}\mathrm{e}}_2{\mathrm{O}}_3+{3\mathrm{H}}_2{\mathrm{S}\mathrm{O}}_4\to {2\mathrm{F}\mathrm{e}}^{3+}+{3{\mathrm{S}\mathrm{O}}_4^{2-}+3\mathrm{H}}_2\mathrm{O}$$

(14)

(4) Relative permeability and capillary pressure calculation model

The calculation model of relative permeability and capillary pressure adopts the Van Genuchten-Mualem model.

3. Numerical model

3.5. Model accuracy and verification

The calibrated accuracy of the concentration field uses the linear correlation coefficient to measure. The correlation coefficients of the three wells ranged from 0.9433 to 0.9990. The simulation accuracy of this study is above 90%. The calibrated accuracy of the concentration field uses the linear correlation coefficient to measure. The expression of the linear correlation coefficient is as follows [35]:

$$r=\frac{{\sum }_{i=1}^n\left(C_i-\overline{C_i}\right)\left(b_i-\overline{b_i}\right)}{\sqrt{{{\sum }_{i=1}^n\left(C_i-\overline{C_i}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(b_i-\overline{b_i}\right)}^2}}$$

(15)

Where, $C_i$ is the simulated concentration, and $\overline{C_i}$ is the mean value of the simulated value; $b_i$ Observed concentration, $\overline{b_i}$ is the mean value of the observed value. Near the simulation area, there are three typical observation holes: W₁, W₂, and W₃. Three observation holes were sampled respectively. The pH value is for the field test, and SO₄^2- was sent to the mining area laboratory for analysis and test.

Figure 8, figure 9, and figure 10 show the calibrate points of each test index concentration of each observation well. From the figures, it can be seen that the simulated values of most points are similar to the observed values. The variety trend of the simulated values is consistent with the observed values, indicating that the model can reflect the actual mining and can be used to predict the ISL.

Figure 5. Comparison between observed and simulated values of at the observation well W₁.

Figure 5. Comparison between observed and simulated values of at the observation well W₁.

Figure 6. Comparison between observed and simulated values of at the observation well W₂.

Figure 6. Comparison between observed and simulated values of at the observation well W₂.

Figure 7. Comparison between calculated values and observed values at the observation well W_3.

Figure 7. Comparison between calculated values and observed values at the observation well W_3.

4. Influence of in-site-leaching of uranium on groundwater quality

4.2. Simulation schemes design and the principle of pollution control

(1) Simulation schemes design

Due to the analysis the impact of pumping ratio and non-uniform injection ratio, five schemes have been designed for simulation (scheme 1 to scheme 5). The water level contours and the streamlines are extracted from these models. The contour map of the water level is the fifth layer grid representing the production zone. The profile of the streamlines has been generated at the central axis. After these numerical simulations, the uranium concentration, pH value, and sulfate radical concentration at the observation site W₁ were selected for comparative analysis. Table 4 displays the simulation schemes.

Table 4. The simulation schemes.

Table 4. The simulation schemes.

	Schemes(Simulation time=1 year)
	1	2	3	4	5
Pumping ratio	0	0.01	0.02	0	0
Non-uniform injection ratio	0	0	0	0.05	0.1
Inner injection rate(m3/d)	176.00	174.24	172.48	184.82	194.09
Outer injection rate(m3/d)	176.00	174.24	172.48	167.18	158.89
Pumping rate(m3/d)	264.00	265.76	267.52	264.00	264.00

(2)The principle of determining the pollution control

After determining the above scheme simulation, the better pollution control scheme selection principle of this study is as follows.

(a) Principle of pollution control based on water level

Based on the numerical simulation results of several schemes, the water table map of the mine layer is extracted respectively. If a scheme can have a higher water head distribution outside the mining area than inside, we believe that this scheme is more optimal.

(b) Principle of pollution control based on streamline

The streamline of the vertical section in Q₂-Q₅-W₁-W₂-W₃-W₄-W₅ is analyzed by comparison. The position of Q₂-Q₅-W₁-W₂-W₃-W₄-W₅ is shown in Fig.3. If the characteristics of the streamline back to the mining area are more obvious, this scheme is the best scheme.

(c) Comparison of concentrations at the observation well.

The concentration simulation value of the W₁ observation well, which is most significantly affected by the mining area, is selected for comparison. The lower the concentration, the more optimal solution we consider.

(3) Sensitivity analysis

Through sensitivity analysis, the parameter with the most significant influence is obtained, which is the priority control parameter. Sensitivity coefficient means how much the dependent variable factor is affected by an independent variable factor. The calculation method is a partial differentiation of the dependent variable concerning the independent variable:

$$S_{ij}=\frac{\partial y_i}{\partial x_j}$$

(18)

Where,$S_{ij}$ is the sensitivity coefficient of the $i$th dependent variable to the $j$th independent variable. In this study, the difference quotient is used instead of the derivative, and the sensitivity coefficient is calculated by the difference quotient.

$$S_{ij}=\frac{y_i\left(x_j+∆x_j\right)-y_i\left(x_j\right)}{∆x_j}$$

(19)

where $y_i$ represents the concentration value of the $i$th index, $x_j$ represents the $\mathrm{j}$th dependent variable, $∆x_j$ represents the increment of the $j$th dependent variable. In this study，$y_i$ was the concentration when the penetration curve was stable.

5. Results and discussion

5.1. The variation of uranium ore,UO₂²⁺ ,SO₄^2-,H⁺ for base scheme(scheme 1)

(1) Variation of uranium ore content and migration of uranium

When the acid solution is injected into the aquifer, the uranium minerals dissolve in the groundwater, which is the main mechanism of in-situ leaching. In Figure 8 a, the blue area is the uranium dissolution zone, and the red area is the uranium precipitation zone. The maximum volume fraction of a uranium increment may exceed 0.006%, but cannot exceed 0.008%. Uranium ore volume fraction can be reduced by more than 0.01%. The results show that the uranium ore is mainly dissolved. According to the simulation results, the uranium ore is mainly dissolved near the water injection well, which is the main source of uranium in the mining area. Near the injection well, the volume fraction of uranium is reduced to 0.01%. Near pumping wells, the volume fraction of uranium can be increased to 0.006%. This is because of the groundwater flows, the acid in the leaching solution is gradually consumed, which will cause the pH value to be decreased so that the previously dissolved uranium will be precipitated again. Outside the mine area, there is no uranium in the aquifer.

From the distribution of uranium concentration, the uranium concentration near the pumping well is higher. The highest concentration of UO₂²⁺ occurs between injection wells and pumping wells. The maximum concentration of UO₂²⁺ can exceed 0.004 mol/kg. (Figure 8 b)

Figure 8. (a) Spatial distribution of volume fraction variation of uranium ore body, (b) the distribution of the concentration of UO₂^2+.

Figure 8. (a) Spatial distribution of volume fraction variation of uranium ore body, (b) the distribution of the concentration of UO₂^2+.

(2) SO₄^2- Spatial distribution characteristics

In the acid-leaching process, the amount of sulfuric acid is large. The concentration of SO₄^2- is a key indicator that needs to be described. In this simulation, the concentration of SO₄^2- is relatively high, with the highest concentration exceeding 0.13mol/Kg. The migration distance can be more than 150m. This is mainly because SO₄^2- is a large amount of injected ions in mining areas, and its consumption is limited. Although some SO₄^2- interacts with minerals, there is still a large amount of redundancy due to the large amount of injection. The concentration spatial distribution of SO₄^2- is shown in Figure 9.

(3) H⁺ spatial distribution characteristics

An important factor affecting ground leaching is H⁺ concentration. Figure 10 shows the change in pH after one year of mining. At the inside of mine, the pH value is relatively low. The pH value increases gradually with the distance from the mining area. Figure 10 also shows that if the lower pH value range is extended to about 80 meters. Figure 11 shows the concentration spatial distribution of H⁺ after one year of mining. The migration amplitude of H⁺ concentration is not as large as SO₄^2-.This may be due to H⁺ being chemically active, which participates in more reactions and increases consumption during the migration process.

6. Conclusions

References

1. Saunders:, J.A.; Pivetz, B.E.; Voorhies, N.; Wilkin, R.T. Potential aquifer vulnerability in regions down-gradient from uranium in situ recovery (ISR) sites. J Environ Manage 2016, 183, 67–83. [Google Scholar] [CrossRef] [PubMed]
2. Bhargava, S.K.; Ram, R.; Pownceby, M.; Grocott, S.; Ring, B.; Tardio, J.; Jones, L. A review of acid leaching of uraninite. Hydrometallurgy 2015, 151, 10–24. [Google Scholar] [CrossRef]
3. Wellmer, F.W.; Becker-Platen, J.D. Sustainable development and the exploitation of mineral and energy resources: a review. International Journal of Earth Sciences 2002, 91, 723–745. [Google Scholar]
4. Poole, T.; Bruneton, P.; Fairclough, M.; Schnell, H.; Valter, O. Uranium Resources as Co- and By-products of Polymetallic Base, Rare Earth and Precious Metal Ore Deposits; International Atomic Energy Agency: 2018.
5. Tulsidas, H.; Fairclough, M. World Distribution of Uranium Deposits (UDEPO). 2018.
6. Edwards, C.R.; Oliver, A.J. Uranium processing: A review of current methods and technology. JOM: the journal of the Minerals, Metals & Materials Society 2000, 52, 12–20. [Google Scholar]
7. Mudd, G.M. Critical review of acid in situ leach uranium mining: 2. Soviet Block and Asia. Environ. Geol. 2001, 41, 404–416. [Google Scholar] [CrossRef]
8. Zhang, C.; Tan, K.; Xie, T.; Tan, Y.; Kong, L.J.P. Flow Microbalance Simulation of Pumping and Injection Unit in In Situ Leaching Uranium Mining Area. Processes 2021, 9, 1288. [Google Scholar] [CrossRef]
9. Zhang, C.; Xie, T.T.; Tan, K.X.; Yao, Y.X.; Wang, Y.A.; Li, C.G.; Li, Y.M.; Zhang, Y.; Wang, H. Hydrodynamic Simulation of the Influence of Injection Flowrate Regulation on In-Situ Leaching Range. Minerals 2022, 12, 15. [Google Scholar] [CrossRef]
10. Collet, A.; Regnault, O.; Ozhogin, A.; Imantayeva, A.; Garnier, L. Three-dimensional reactive transport simulation of Uranium in situ recovery: Large-scale well field applications in Shu Saryssu Bassin, Tortkuduk deposit (Kazakhstan). Hydrometallurgy 2022, 211, 105873. [Google Scholar] [CrossRef]
11. Mibus, J.; Sachs, S.; Pfingsten, W.; Nebelung, C.; Bernhard, G. Migration of uranium(IV)/(VI) in the presence of humic acids in quartz sand: a laboratory column study. J Contam Hydrol 2007, 89, 199–217. [Google Scholar] [CrossRef]
12. Kurmanseiit, M.B.; Tungatarova, M.S.; Kaltayev, A.; Royer, J.J. Reactive Transport Modeling during Uranium In Situ Leaching (ISL): The Effects of Ore Composition on Mining Recovery. Minerals 2022, 12, 21. [Google Scholar] [CrossRef]
13. Yabusaki, S.B.; Fang, Y.; Waichler, S.R. Building conceptual models of field-scale uranium reactive transport in a dynamic vadose zone-aquifer-river system. Water Resources Research 2008, 44. [Google Scholar] [CrossRef]
14. Kohler, M.; Curtis, G.P.; Kent, D.B.; Davis, J.A. Experimental Investigation and Modeling of Uranium (VI) Transport Under Variable Chemical Conditions. Water Resources Research 1996, 32, 3539–3551. [Google Scholar] [CrossRef]
15. Zhang, X.; Jiao, C.; Wang, J.; Liu, Q.; Li, R.; Yang, P.; Zhang, M. Removal of uranium(VI) from aqueous solutions by magnetic Schiff base: Kinetic and thermodynamic investigation. Chemical Engineering Journal, 2012; 198-199, 412–419. [Google Scholar]
16. Bear, J.; Jacobs, M. On the movement of water bodies injected into aquifers. Journal of Hydrology 1965, 3, 37–57. [Google Scholar] [CrossRef]
17. Xu Shaohui, Z.X. Hydraulic interception technology and numerical simulation of groundwater oil pollution control. JOUR NAL OF HYDRAULIC ENGINEERING 1999, 71–76. [Google Scholar]
18. Javandel, I.; Tsang, C.F. Capture-zone type curves-A tool for aquifer cleanup. Ground Water 1986, 24, 616–625. [Google Scholar] [CrossRef]
19. Faybishenko, B.A.; Javandel, I.; Witherspoon, P.A. Hydrodynamics of the capture zone of a partially penetrating well in a confied aquifer. Water Resources Research 1995, 31, 859–866. [Google Scholar] [CrossRef]
20. Bair, E.S.; Roadcap, G.S. Comparison of flow models used to delineate capture zones of wells.1.leaky-confined fractured-car bonate aquifer. Ground Water 1992, 30, 199–211. [Google Scholar] [CrossRef]
21. Wang, J.C.; Booker, J.R.; Carter, J.P. Analysis of the remediation of a contaminated aquifer by a multi-well system. Computers and Geotechnics 1999, 25, 171–189. [Google Scholar] [CrossRef]
22. Hudak, P.F. Numerical modeling assessment of three-gate structures for capturing contaminated groundwater. Environ. Geol. 2008, 55, 1311–1317. [Google Scholar] [CrossRef]
23. Zheng, H.Q.Y. Site conditions experiments for field test on in-situ leaching and hydrodynamic simulation of uranium in Bayan-Uul. University of Technology(Nanchang),East China. 2017.
24. Chang, Y. Simulation study on groundwater dynamics control of leaching range of in-situ uranium well field. University of South China 2020. [Google Scholar]
25. Chen, Q.Q. Numerical Simulation of In-situ Leaching of Uranium by Acid Method Based on Reactive Solute Transport Theory. University of Technology(Nanchang),East China. 2020.
26. Xu, T.; Spycher, N.; Sonnenthal, E.; Zheng, L.; Pruess, K. TOUGHREACT_V2_Users_Guide. 2012.
27. Xu, T.; Apps, J.A.; Pruess, K. Reactive geochemical transport simulation to study mineral trapping for CO2 disposal in deep arenaceous formations. Journal of Geophysical Research: Solid Earth, 2003; 108. [Google Scholar]
28. Pruess, K.; Oldenburg, C.M.; Moridis, G.J. TOUGH2 User's Guide Version 2. United States: N. p. Web. 1999.
29. Xu, T.; Spycher, N.; Sonnenthal, E.; Zhang, G.; Zheng, L.; Pruess, K. TOUGHREACT Version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions. Computers geosciences 2011, 37, 763–774. [Google Scholar] [CrossRef]
30. Dobson, P.F.; Salah, S.; Spycher, N.; Sonnenthal, E.L. Simulation of water-rock interaction in the yellowstone geothermal system using TOUGHREACT. Geothermics 2004, 33, 493–502. [Google Scholar] [CrossRef]
31. Sonnenthal, E.; Ito, A.; Spycher, N.; Yui, M.; Apps, J.; Sugita, Y.; Conrad, M.; Kawakami, S. Approaches to modeling coupled thermal, hydrological, and chemical processes in the drift scale heater test at Yucca Mountain. International Journal of Rock Mechanics and Mining Sciences 2005, 42, 698–719. [Google Scholar] [CrossRef]
32. Spycher, N.F.; Sonnenthal, E.L.; Apps, J.A. Fluid flow and reactive transport around potential nuclear waste emplacement tunnels at Yucca Mountain, Nevada. J Contam Hydrol 2003, 62-63, 653–673. [Google Scholar] [CrossRef] [PubMed]
33. Aradottir, E.S.P.; Sonnenthal, E.L.; Bjornsson, G.; Jonsson, H. Multidimensional reactive transport modeling of CO2 mineral sequestration in basalts at the Hellisheidi geothermal field, Iceland. International Journal of Greenhouse Gas Control 2012, 9, 24–40. [Google Scholar] [CrossRef]
34. Audigane, P.; Gaus, I.; Czernichowski-Lauriol, I.; Pruess, K.; Xu, T.F. Two-dimensional reactive transport modeling of CO2 injection in a saline Aquifer at the Sleipner site, North Sea. Am. J. Sci. 2007, 307, 974–1008. [Google Scholar] [CrossRef]
35. Zheng, C.M.; Bennett, G.D. Applied Contaminant Transport Modeling(Second Edition); 2002.