The development of reactive transport soared during the 1990's, driven by the necessity to demonstrate the long-term efficiency of radioactive waste repositories (Bildstein et al. 2019; Claret, 2019; Cama et al. 2019, both this volume). The approach, based on a rigorous description of processes and their coupling, provides a basis of confidence to bridge the gap in time and space between knowledge gained in laboratory experiments and the dimension and lifetime of a high-level waste repository. Reactive transport codes progressively increased in complexity, as more processes were included to account e.g., for variably saturated flow, heat transport and more complex chemical processes (Steefel et al. 2015). In the meantime, improved algorithms and computer power opened the way for larger simulations.

Reactive transport simulation now embraces a large array of applications in the Earth and environmental sciences. Thus, in the field of contaminant hydrology, such studies address problem from the interface scale, to quantify migration and retention mechanisms, up to the scale of a watershed integrating multilayer reservoirs and their inherent spatial variability. Over the last decades, extensive monitoring and characterization studies allowed for better calibration of processes, parameters determination and more generally better quantitative understanding of sites (e.g., Hammond and Lichtner 2010; Li et al. 2010; Hammond et al. 2011; Arora et al. 2016) or on less emblematic sites (e.g., Mangeret et al. 2012; Estublier et al. 2013). As a result, better understanding of mechanisms, quantification of competing retention and migration processes were achieved. They allow for the estimation of plume migration, or the evaluation of remediation strategies.

In these applications, a unique application site is characterized, then a model is built in order to draw conclusions and propose recommendations. Apart from academic or agencies sponsored studies, the industry is also interested by this approach: estimation of reservoir properties in the oil and gas industry (Xiao and Jones 2006), evolution of gas storage sites (e.g., CO2 in Audigane et al. 2007, or H2 in Saínz García et al. 2017), evaluation of impacts of the mining industry (Metschies and Jenk 2011; Molson et al. 2012; Neuner and Fawcett 2015).

Another type of use is also relevant: industrial deployment of reactive transport within operation with a view to support operators in decision-making. Operating conditions are then very different: models should be built on a regular basis on numerous application cases, sometimes relying on incomplete sets of data. Additionally, this type of application only makes sense if it can be integrated in the regular workflow of operations. Hence, timing (to keep pace with the operations), and interfacing (with other aspects of production) are critical. Finally, such deployment can only be performed if the value of the simulations for the operation can be demonstrated.

An illustration of industrial use of reactive transport in operations is given: exploitation of uranium by in situ recovery (ISR). The technique is best suited for large, low grade, medium depth deposits. It consists in dissolving the metal in place by circulating solutions, using a large number of injection and production wells to cover the surface of the field. Relying on the circulation of fluids and their interaction with the deposit, the technique is particularly suitable for reactive transport analysis (Nguyen et al. 1983; Regnault et al. 2012). Although models can be built to understand the general behavior of such systems, when quantitative assessment of a unit production is wanted, the size of the exploitation and the strong spatial variability inherent to the formations require to build models specific to the local conditions (geology and constraints of exploitation) in each portion of the deposit.

After a brief description of the ISR technique, the use of reactive transport simulation is exemplified on the KATCO mine. The modeling approach is presented, including the data necessary to calibrate the model. The conditions to transform a mostly academic application into an operational tool are commented; some examples of use in operation are then proposed. Finally, in addition to exploitation, reactive transport modeling is also relevant for controlling environmental impacts. Although this application is more traditional, it shows that reactive transport modeling deployed in this context can create the conditions for an integrated approach to better identify post-operational constraints, or even to test alternative exploitation strategies limiting the need for active remediation a posteriori.


In situ recovery (ISR) mining is defined as the extraction of a metal from its host formation by the circulation of chemical solutions (lixiviants) using injection and extraction wells from the surface (Fig. 1). The solution is treated in surface plants to recover the metal. Reagents are then adjusted in the solution before recirculation in the ore formation (IAEA 2001). ISR was developed simultaneously by the USA and the USSR in the early 1960's for the recovery of uranium in secondary ore bodies in permeable aquifers. It has become one of the standard production methods, with a share exceeding 50% of world uranium production in 2014 (IAEA 2016). The technique is particularly suited for large ore bodies, at intermediate depth and low concentrations, like roll-front deposits (Dahlkamp 1993; Kyser 2014; Saunders et al. 2016). High permeability and presence of confining formations for the deposit are also desirable.

Uranium solubility increases markedly under oxidizing conditions, at low or high pH. Complexation with several anions further enhances solubility: particularly relevant for in situ context are SO42 and CO32−. Two technological choices are then available (Bhargava et al. 2015): acid (usually sulfuric acid) or alkaline (higher pH carbonate or bicarbonate solutions) leaching. In both cases, oxidizing agents can be added: air or oxygen, hydrogen peroxide, ferric iron and others.

The selection of technique is guided by regulatory aspects (mostly driven by remediation process opportunities), geological considerations (notably presence of mineral buffers, e.g., calcite) and economical optimization. However, acid leaching (Kazakhstan, Australia) is de facto the preferred technique, representing up to 96% of total uranium produced by ISR.

The fundamental advantage of ISR over other mining techniques is reduced costs: it is by far the most cost effective extraction technique (Kidd 2009). By bringing the reagents directly in the ore body, bulk material handling is reduced to the minimum. The absence of a mining fleet, milling capacity and heavy preparatory work (removal of the overburden for open-pit mines, development of access shafts for underground mining) drastically reduces upfront capital cost (Heili 2018). ISR mining plans are also highly flexible: new fields can be quickly developed or alternatively production rates can be reduced without substantial impact on production costs. Not only this progressive development allows for reduced return on investment time, it also provides for dynamic adjustment to the demand.

Reduced environmental footprint is another key advantage of ISR (IAEA 2005). Since liquid is moved rather than rock, disturbance at the surface is minimized (IAEA 2001). Hence, a major benefit is the absence of tailings; the quantity of solids generated by the exploitation is mostly limited to cuttings from the (extensive) well drilling. Other interesting features are limited water consumption, limited impact on landscape, and low impact on air quality (no milling operations). The major concern during operations is then the risk of contamination of groundwater by the leaching solution (Mudd 2001). Surface risks can be limited operations by the application of best practices; migration from the ore formation is limited when ISR is applied only in confined aquifers.

Closure operations are also facilitated since no legacy sites are created. The key point in post-mining is the restoration of water quality in the ore formation, and its a priori demonstration. Provided that good practices are followed and depending on the local geology and presence or absence of buffer minerals, monitored passive remediation or active restoration (e.g., pump and treat) can be chosen (IAEA 2016).

With the technical challenges and economic potential associated with ISR, this domain has triggered a lot of research interest (Märten et al. 2013; IAEA 2016; Le Beux et al. 2018): e.g., better assessment of the ore geometry, development of monitoring tools, increased valorization of production data, improvement of well engineering, economic and hydraulic optimization of mining plan. Simulation also plays an important role: hydrological software tools are used for well-field planning and operational control. Simulations, coupled with careful characterization of the formation, can help maximize the contact between the leaching fluid and the ore, reduce the risk of fluid migration out of the well-field, and micromanage injection rates including role reversal between injectors and producers.

Along these developments, reactive transport modeling has been identified as a significant tool to quantify impacts, support the evaluation of remediation techniques, and help demonstrate their compatibility with regulatory targets (Kalka et al. 2006; Descostes et al. 2014; Johnson and Tutu 2016; Dangelmayr et al. 2017). Several authors also mention the use of reactive transport to help understand and quantify the processes at stake during exploitation (Kalka et al. 2006), although full 3D reactive transport at the scale of production units is still scarce (Regnault et al. 2014, 2017). The model developed in two applications is presented hereby. The quantitative results of the simulation can be used to improve the understanding of processes and their coupling: this offers in itself some opportunities to guide tentative optimizations. Furthermore, the simulation can also be used by engineering departments to directly test optimization ideas, improve production forecasting, and more generally help in planning: some examples illustrate how intensive use of reactive transport simulation can create value during and after exploitation.


Reactive transport model for the KATCO Mine

The KATCO mining company is a joint venture between Orano Mining (51%) and Kazakhstan national mining company Kazatomprom. It specializes in the in situ recovery of uranium, in roll-front type formations, in Southern Kazakhstan. With a 4,000 t uranium yearly capacity, it became in 2009 the world largest ISR mine. Its total production since 2006 amounts to 36,000 t uranium. KATCO operates on two distinct sites in the Chu-Saryssu basin: Tortkuduk and Muyunkum, 40 km apart, each with their own well-field (4,000 wells in activity at any given time) and separation plant.

The construction of reactive transport simulations at the block scale relies on a combination of hydrogeochemical description of the reservoir, and identification of main processes. Spatial variability is inherently high in this type of deposits; its huge impact on the shape of the recovery curve imposes a careful evaluation and integration in the simulation workflow. All simulations were performed using the HYTEC code developed at MINES ParisTech (van der Lee et al. 2003; Lagneau and van der Lee 2010).

Hydrogeochemical description.

The concession around the Muyunkum deposit is composed of several aquifers containing roll-front accumulations, permeable sandy formations from the Paleogene era (Petrov 1998). Particularly, the Uyuk aquifer is an Eocene sandy aquifer, mildly cemented, with a low content in calcite (< 0.1 wt%), approximately 400 m deep and 30 m thick. The aquifer presents a high spatial variability, both vertically and laterally, as a result of its sedimentation history: coarse to fine sands channels deposited under tidal influence, presence of clay intercalations stemming from continental flooding plains.

Hydrogeological characterizations show porosity and permeability in the sand ~ 20–25% and 10−4 m/s respectively (ben Simon 2011). With head gradients around 1‰, the natural regional flow velocity is around 5 m/y, from the oxidized to the reduced part of the aquifer through the uranium bearing roll-front.

Hydrogeochemical analyses show a clear discrimination between oxidized (pH ~ 7.6, Eh ~ 50 mV/SHE), and mineralized or reduced zones in the aquifer (pH ~ 7.8, Eh ~ −150 mV/SHE). The solid is a mildly cemented arkosic sandstone (Munara 2012; Robin et al. 2015a), with composition around 70% quartz, 10% feldspar, 5% micas, 5–10% clay minerals (mostly smectite, and kaolinite). Authigenic smectite is mostly composed of Beidellite, with a significant content in Fe(III) (Robin et al. 2015b). Other minerals include goethite or pyrite (< 1% in the oxidized or reduced part of the aquifer respectively), and low concentration of calcite (< 0.1%). Finally, uranium-bearing minerals are found in the roll-front, mostly uraninite and coffinite, at relatively low concentrations (~ 0.1%).

Initial porosity is around 20% in the sandstone. Porosity evolution remains small in the reservoir, limited to 1 to 2%, with the weak reactivity of most minerals. Actively reactive minerals (calcite, uraninite) represent a very small volume fraction. Likewise, precipitation of secondary minerals (mostly gypsum and possibly aluminum hydroxysulfates) has a limited impact on porosity if precipitation if spread over the reservoir, although chemical clogging can occur when precipitations are concentrated close to the production well.

Initial aquifer water composition was prescribed using base line monitoring: geochemistry is controlled by equilibrium with the host rock minerals, with total dissolved solids around 1 g/L at temperature around 20 °C. The injection solution composition evolves through the simulation. An initial acidification phase uses aquifer water acidified with sulfuric acid at typically 20 g/L. When pH drops below 2 at the producers, the block is connected to the plant. Fluid composition then stems from the multiple circulations of the fluid over the whole well field: the result is an accumulation of dissolved elements in the solution (only uranium I stripped at the plant). Typical solutions reach up to e.g., 0.5 g/L in Ca, or 1 g/L in Al. The model uses compositions in line with analyses at the plant, and acid content is adjusted daily according to the operator prescription: acid content is adjusted so that pH at the production well remains between 1 and 2, to ensure Fe(III) mobility and limit precipitation of secondary mineral; usually, the operator starts production with high acid content (typically 20 g/L H2SO4), then acid content is reduced for older, less producing, blocks (down to 5 g/L).

Chemical processes.

Uranium dissolution during acidic leaching is dominated by the oxidation of U(IV) minerals (e.g., uraninite UO2, coffinite USiO4) by Fe3+ in the bearing solution. The origin of ferric iron is multiple: local dissolution of gangue minerals (goethite, beidellite), recirculation of Fe3+ from the well field, or surface active regeneration of Fe3+ using peroxide (e.g., Heathgate operations in Southern Australia, Märten 2006) or other oxidants. Low pH (typically < 2) is required to allow for ferric iron mobility.

The governing reaction for uranium dissolution is then:


The reaction is controlled by the availability of ferric iron, solubility of uranium (in which complexation, notably with, SO42 can play an important role) and kinetics. In this case, it was verified that oxidant availability controls uranium dissolution, so that the specific (reduced uranium) mineral is not a determining factor. In other cases, several uranium-bearing mineral should be carefully considered, including possibly oxidized uranium minerals (e.g., schoepite in other geologic contexts).

Reactions with gangue minerals are also key, since they can control pH or local sources of ferric iron (Robin et al. 2016):


Sorption can also play an important role: clay surface sites can constitute an important sink term for protons. However, this effect is limited mostly to the initial acidification phase of the exploitation, when pH is driven from initial near neutral conditions to less than 2. Surface sites then remain mostly saturated in protons, independently of pH evolution during the whole exploitation phase. Their influence is again crucial post exploitation, when the remaining buffer capacity of the aquifer drives pH above 4.

Column tests are ubiquitous in ISR operations, a quick, easy means to determine the leachability of a sample ore. As already mentioned, direct use of the results to infer production curves should be exerted with care, since they are only representative of the conditions in which they were performed. Moreover, representativeness of a (small) sample compared to large heterogeneous ore bodies is limited. Finally, preservation and preparation of samples can significantly alter the behavior of the system; e.g., ben Simon et al. (2014) used a mixture of different samples from different parts of ore formation; moreover, oxidation between sample collection and test column led to a significant evolution: very low pH due to pyrite oxidation with air O2, oxidation of nearly half the uranium in the system. Nevertheless, these column tests offer perfect opportunities to test and calibrate the reactive paths and kinetics in controlled conditions. The chemical system can be further tuned using production data.

Hydrodynamic model.

The ISR process relies on the recycling of the leaching solutions: after recovery at the producers, the solution is treated in the surface plant to strip uranium, then the reagents concentration (acidity, oxidizers) is adjusted before re-injection in the well field. Over time, the dissolved content of the solution increases. At such low pH, Ca, Mg or Al content can reach values as high as 0.5 to 1 g/L, and over 20 g/L for sulfate. Therefore, despite high flow rates imposed at the wells (several m3/h over 10 m screens), density differences between the leaching solution and the initial aquifer water are sufficient to create downwards migration: this was verified with observation wells in test production cells and confirmed by sensitivity analysis flow simulations (Bonnaud et al. 2014).

A gravity driven (saturated) flow solver is then used for the simulations. Density can be deduced from solution composition. The 3D model is built using the geometry of the well field (position of the wells, length and depth of the screens), on the whole thickness of the aquifer. Source terms are prescribed at each well, both in terms of flow-rate and solution composition: ideal conditions can be used, or they can be set to reproduce actual production data for history matching exercises. Boundary conditions can be imposed to allow for regional flow, although it has been demonstrated that the effect is limited.

Spatial variability.

The ore genesis process is responsible for a high spatial variability of the system. Sedimentation creates a first stage of heterogeneity, with sand channels of varying grain size and discontinuous clay interlayers. Flow and oxidation in the aquifer during the creation of the roll-front leads to a second stage of heterogeneity. Indeed, the roll-front is an interface between an oxidized upstream and a reduced downstream aquifer. The migration of the interface is therefore directly linked to the local reducing buffering capacity of the aquifer (organic matter, pyrite) and the flux of oxidizers i.e., the flow field (and in fine spatially variable permeability and porosity field). The sequential processes account for a very complex 3D shape for the roll-front. Far from a typical roll-front shape, the mineralized zone is a complex 3D structure, immerged in oxidized and reduced zones and intersected by low permeability clay barriers. Finally, uranium grade within the mineralized envelop is also variable, again a result of the complex history of the front migration.

The correlation length for uranium grade is typically a few 10's meters, similar to the size of the production cells (radius 42 m). Therefore, facies and uranium grade distribution play a huge role in the shape of the response function of the production cell. Particularly, smoother distributions tend to yield sharper production curves. On the contrary, high heterogeneity can present high concentration clusters that take longer to dissolve (mainly due to hydrodynamic control), resulting in very elongated residual production (Lagneau et al. 2018). Another aspect of the variability within the ore body is the diversity of response over the permit. Indeed, local differences in geometry of the mineralized envelop, grade distribution or clay content can lead to very different evolution of pH or uranium concentrations in the producers.

A consequence of the variability is that homogeneous models, or even 1D and 2D descriptions come short to being fully predictive. Reactive transport simulations must therefore tackle a 3D description and incorporate an accurate description of the geology, with a block model both in terms of geochemical facies, permeability, and ore grade.

Model description.

The workflow was applied on a particular block on the KATCO operation: Block A1. The block is composed of 8 adjacent production cells. Due to operation choices, mostly to enhance production in rich areas, several production cells contain 2 or more producers. In areas with two levels of mineralization separated by clay interlayers, multiple screens were also positioned. The model geometry then takes into account 18 producers and 43 injectors. Production data were used to individually constrain the flow rate and injection solution in each well. Both prescriptions evolve through the operations. Solution composition (particularly acid content) is a key lever for the operator. Injection and production flow-rates are regularly adjusted by the operator; they are also prone to modifications due to well evolution (clogging), maintenance (work-over), and flow allocation constrained at the mine scale by the plant specification.

One realization of the block model was chosen to constrain the variability. The implications of this simplification are discussed further.

The geochemical model was adjusted against production history, mostly pH, ferric iron balance and uranium concentration in the producers. During this process, an evaluation of the hierarchy of processes was performed. The aim was to determine the main controlling reactions; second order reactions were then removed in order to reduce the complexity of the system. The minimal mineralogical assemblage and associated kinetic parameters (from Palandri and Kharaka 2004; Robin et al. 2016) for each chemical unit are given Tables 1 and 2.

The simulation gives access to a wide array of results: e.g., Figure 2 shows the uraninite concentration in the system after 600 days of exploitation. The progression of uranium recovery is noticeable by the depletion within the limit of the block. The yellow cloud shows the extension of the acidification (pH< 1.8). This 3D representation is useful to identify and localize possible defects in the exploitation scheme: poor leaching, residual high pH areas, slow dissolving high grade clusters, …

Another option is to plot the concentrations at the production wells, for each producer or in the collector (i.e., a mixing of all producer wells in the block). This visualization is particularly relevant for the operator as it corresponds to the only information available on site to monitor the exploitation (daily analysis of the solution). Such representation is proposed in Figure 3: pH and dissolved uranium concentration in the collector of all the producers from Block A1. This was used to adjust the parameters of the model against production data.

As can be seen, a very good fit was obtained, even after various events during the exploitation. The adjustment is facilitated since most parameters are constrained by measures (permeability, mineral concentrations), databases (solubility constants, rate constants), geometry of the system (well field, block model), and operations (flow-rates, composition of the leaching solution). In the end, the adjustment was performed on poorly determined parameters: calcite concentration (too low in the system to allow for an accurate determination), beidellite reactive surface area, and ferric iron content of the beidellite. The adjustment is fairly easy to perform, since the three parameters have mostly independent impacts on the system. Indeed, calcite concentration is the first determining parameter (with porosity) on the time to obtain the pH drop. Beidellite (kinetically controlled) dissolution is the major long-term acid consumer, so that beidellite surface can be adjusted on the pH values or acid balance after one year. Finally, once beidellite kinetics is adjusted, its content in ferric iron controls the release of oxidizer within the block and the increase in uranium dissolution.

Interestingly, three parameters only are then sufficient to reproduce the complexity of the reaction of the system to multiple events: evolution due to the acidic and oxidizing solution, modification of the acid content in the injected solution, modification of flow-rates, and redrilling.


The model was verified by application to another block in the vicinity: Block A2. For this application, the model was obviously adapted: one realization of the block model over Block A2 area, position of the wells, exploitation scenario (flow-rates, composition of the leaching solution). However, the geochemical model calibrated from Block A1 was used without any adjustment. Figure 4 shows the simulation results for Block A2: the model is very close to the production data, providing a good validation of the model.

Interestingly, the two blocks behave very differently. Block A2 displays a sharper uranium peak, with high peak concentrations: these differences are mostly due to the local distribution of uranium in both blocks. Also, the two blocks were exploited sequentially, so that part of the acidic plume from Block A1 invaded block Block A2, resulting in lower initial pH.


The strength of the model relies on its roots in processes and geometric description of the system. As a result, very few adjustment parameters are needed. The operator currently strives to acquire these missing parameters in order to increase the predictive capacity of the model. Also, although chemistry evolves over the width of the mining concession (~ 20 km), the wavelength is sufficiently large that adjustment on a block can be used without modification in adjoining blocks. However, the adjustment needs to be performed again in farther areas, where chemical specification may be significantly different.

It was mentioned that spatial variability was key to control the shape of production curve. And yet, correct adjustments were possible using a single realization of the block model. Use of another realization would have given a similarly good adjustment, perhaps with marginally altered parameters, provided the total amount of uranium in the block was similar. There seems to be a contradiction between the two remarks: spatial variability is key but simulation results seem independent of the realization. This conflict disappears when looking as the relative dimension of the range of the geostatistical model for uranium grade (~ 50 m) and the block (250 –400 m). The block size is well above the range, so that stochastic variability is averaged over the size of the block. In the end, at the block scale, the structure of the distribution is more important than the actual distribution. This is not true at the size of a unit cell (radius 42 m): the simulation result for a specific producer (as opposed to the collector of the block) changes radically for different realizations. Langanay et al. (2018) provide some insights into the impact of geological uncertainty on the simulated uranium recovery.


The demonstration provides solid foundations for the reactive transport simulation of ISR operations. A direct benefit is a better understanding and quantification of the processes at stake over the lifetime of a block. Indeed, the information available to the operator is very limited: after the initial geologic reconstruction, the main flux of information during operations comes from the monitoring of flow-rates, composition of injected fluids, and chemical analysis of produced solutions for each well on a daily basis. This information is very useful for operation management: at the mine scale, the operator needs to arbitrate resource allocation (acid consumption and overall flow manageable by the plant). However, deconvolution of fluid composition is hard to perform, since it stems both from the fluid evolution along a flow tube and the mixing of all the flow tubes converging to a producer.

Once the model is validated, it can be used to evaluate the controlling processes in the system. For instance, kinetically controlled beidellite dissolution is the first acid consumer on the long term, over 80% of acid balance over the lifetime of a block. The weak catalyzing effect of protons, and conditions always far equilibrium at pH< 2, mean that dissolution rate is only weakly dependent on pH. Acid consumption (and operational costs) could therefore be reduced using faster injection rates. However, this should be compared to the impact of lower local production of ferric iron associated with beidellite dissolution. Another investigating concerns the dissolution regime of uraninite: depending on the position in the production cell, prescribed flow-rates and availability of reactants (Fe3+), uraninite dissolution can be controlled by kinetics or hydrodynamics. Here again, quantitative analysis can provide useful insights on a possible lever to optimize uranium production.

This kind of quantitative analysis using reactive transport codes, based on carefully devised model of geologic systems, is already performed in a large array of applications. However, a more intensive use of simulations can benefit the operations. Indeed, it is clear that a block behavior is heavily constrained by the local geology. A simulation calibrated on a specific block cannot therefore be used to illustrate a typical production curve. Going beyond generic recommendations, systematic simulations for each block are therefore needed. In this case, academic codes reach their limits.

The strength of a code developed in an academic environment is its versatility: continuous development offers a wide array of options to choose from, which allowed performing the demonstration. The downside of the versatility is a complex combination of options (and keywords in text-based input files), most of them probably not needed for the specific application. Also, although importing input parameters from the operator into the code is manageable for a single application, it can become cumbersome for routine use. Finally, systematic simulations are typically performed by engineering teams. Rapid turnover and dispersion between multiple tasks imply that training should be minimal.

Early in the development of the project, the need for an engineer tool (as opposed to a research code) was identified. A graphical user interface is obviously required, but not sufficient. Several specifications were identified: clear interface to minimize training, easy integration in the operator workflow, presentation of outputs to provide ready-to-use data. A GUI for HYTEC, specifically devised for ISR needs was then created: HYSR. The interface only proposes options relevant to the simulation of ISR exploitation. Moreover, the controls are organized to match the practice of engineers; they do not necessarily reflect numerical simulation thought-patterns. For instance, the interface does not ask for initial or boundary conditions. Rather, it imports a block model or well specifications, which are the standard data an ISR engineer is used to manipulate. Imported data, directly from the operator workflow, are converted by the preprocessor to launch the reactive transport simulation.

Output data management is also a key element to facilitate the use in operational conditions. 3D visualization of the simulation results is still possible, using ParaView (Ahrens et al. 2005). This allows for fast evaluation of leaching efficiency, particularly with application of preselected filters. HYSR also transforms the outputs into production curves: the curves show information at the producer wells or globally for the block, which directly relate to production data. Also, mass balances are performed at the scale of the block. Finally, the net present value (NPV) of the block is proposed: the calculation is based on a cost model integrating capital costs (construction of the block, connection to the plant), operation costs (power use, fraction load at the plant), reagent costs (acid consumed), and revenue (uranium produced). The combination of output formats offers different opportunities for the operator: identification of poorly productive areas to help propose corrective strategies, direct comparison to production data for validation or use in mining plan elaboration, or a robust means to compare alternative exploitation strategies.

Examples of application.

Mine planning for ISR exploitation consists in a predictive simulation of the sequence of well-field operations and construction works needed to reach the annual target of production. Planning also needs to respect regulatory constraints and financial limits imposed by the shareholders. The operator differentiates short- and long-term planning. Short-term planning (12–18 months) is the foundation for construction and revision of the company's budget for the current or following year. Long-term planning is a projection for the whole life expectancy of the mine. It provides a guide for long-term exploitation strategy and data to dimension infrastructures and associated CAPEX. In both cases, the quality of the mining plan depends on the reliability of estimated behavior of the system.

Short-term planning.

In the Katco operations, lifetime of exploitation blocks is between three and five years, with around 60 blocks running at any given time. Past experience of one-year projections shows that mature running blocks (i.e., in quasi stationary state after ramp-up) amount to roughly 80% of extracted uranium and 70% of acid consumption. The prediction of the end-of-life of mature blocks is therefore a major input data for short-term planning: it drives the gap between expected and target production, and therefore the sequence of new blocks needed to reach the target. Predictions of block end-of-life behavior also allows anticipating block closure, on criteria like expected uranium concentration or recovery rate. This information is needed to allocate resources (total flow capacity of the plant or main pipes) over the well-field and assess potential availability for new blocks. A correct assessment of the future behavior of currently operated blocks is therefore key for the robustness of the provisional budget allocated for new blocks development: earthworks on future exploited zones, wells digging, and development of hydraulic network.

The well field behavior is monitored and simulated at each technological block, for three main parameters: flow in the block, dissolved uranium concentration at the collector, and acid consumption. Traditionally in ISR operations, evolution is evaluated independently for the three parameters. Acid consumption is based on empirical coefficients (e.g., mass of acid consumed by mass of uranium produced), based on supposedly similar areas of the exploitation. Evolution of uranium concentration in the leach solution is based on analytical functions using parameters adjusted on production data: the extrapolation of these functions is used as predictor.

Such approach is acceptable as long as operating conditions are maintained: i.e., small evolution of the flow-rate and small variations of the acid prescription curve in the injected solution. Following the standard conditions, technological blocks display the characteristic bell curve (e.g., Figures 3 and 4 for blocks Block A1 and A2). However, this approach looses robustness when changing or complex operating conditions are applied. And yet, such modifications, planned (e.g., acid availability due to overall allocation decision) or provoked (e.g., additional wells, optimization decision) by the operator should be taken into account to allow for a precise short-term prediction.

The 3D reactive transport simulation is here particularly adapted. It offers the possibility to precisely adapt evolving conditions and evaluate their impact on uranium concentration and acid consumption. Better anticipation of alternative strategies (optimization, reaction to unforeseen constraints) can therefore be integrated in the short-term mining plan. Most common actions to increase uranium recovery are increasing acid content in the injected solution, and improving flow rates. Improved flow rates in a production cell can be obtained by work-over on ageing production wells or by additional wells drilling.

An example of reactive transport simulation used in short-term planning is illustrated on block Block B (Fig. 5). The block life started with a classical production curve for the first 18 months, with a sharp uranium peak followed by a rapid decrease (phase A). Since the block had important estimated residual reserve, the operator tested several strategies to increase uranium yield. First, the operator decided to increase acid content in the injected solution (phase B), which stabilized uranium concentration at the collector. The operator then decided to improve production flow rate by repeated work-over on several production wells (to counter well clogging), and drilling of additional production wells (phase C). The impact was an important increase in uranium concentration at the collector, which stayed on a plateau for two years. Finally, when estimated optimal recovery rate was reached, the operator decided to reduce the amount of acid injected in the block (phase D), which resulted in a rapid decrease of uranium concentration.

The analytical function, with parameters adjusted on production data, yields correct results during the first phase of exploitation (Fig. 5, top), but fails when operating conditions deviate from the standard procedure (phase B, C).

On the other hand, the reactive transport simulation produces a correct estimation of the whole production life (with a unique set of parameters): Fig. 5, bottom. The link between injected acid and uranium concentration is correctly anticipated thanks to the geochemical model used: the increased acidity in phase B induces lower pH, faster kinetics for beidellite dissolution, and finally increased Fe3+ availability. Uranium oxidative dissolution in the ore body is then increased resulting in high uranium concentration in the producers. The opposite is observed when the operator decreases the acidity in phase D.

The model provides correct estimation of the impact of evolving velocity field within the ore body (phase C). Improved injection and production rates increase the acidification of the ore, including mineralized areas far from the producers. Also, production of dissolved uranium is improved by minimizing uranium-rich solution loss, particularly in the lower parts of the reservoir. Finally, additional production wells in high-grade areas of the block stabilize uranium concentration in the collector by adding high yield lines to the overall leached solution. The simulation correctly predicts the impact of this combination of factors. The 3D approach is here a necessity: the importance of uranium grade spatial variability is demonstrated here by the impact of local modifications (additional wells) on the overall behavior of the block.

Long-term planning.

The strength of reactive transport in short-term planning is the capacity to rapidly anticipate the impact of fluctuating operational conditions at the block scale. For long-term planning, the key is to produce simulations, on a large scale, to help propose a sequence of production all the way to the end of exploitation.

By construction, long-term planning is mostly constrained by the predictive behavior of future blocks. Geological and grade data on these areas is usually coarse: well-spacing on the exploration grid (typically 100 –400 m) is larger than at the technological stage (around 40 m). Moreover, the planning team cannot rely on the operational feedback on undeveloped areas: at best, a pilot operation can provide some indication on acid consumption or uranium yield on a small scale (a few cells) and on a limited operation time (several months). In these conditions, semi-analytical recovery functions are unsatisfactory: indeed, their parameterization relies on (unavailable) production data or analogy on similar blocks produced in the same area.

On the other hand, reactive transport simulation can be used even if the model is not completely constrained: a 3D geological model can be derived from (large grid) exploration data, uranium and gangue mineral behavior can be inferred from the geochemical model based on (limited number) mineralogical analyses. The simulation provides an interesting tool to illustrate and quantify the impact on production of certain features of the ore body. Thus, geometry has a huge impact, which can be calibrated by simulations: thickness of the reservoir, density of clay interlayers, morphology of the mineralization and localization in the reservoir. For instance, recovery is facilitated for compact ore bodies in confined low thickness reservoir. Conversely, dispersed mineralization in a thick reservoir leads to reactant dispersion, increased acid consumption, and lower recovery due to uranium loss.

Figure 6 illustrates the behavior of technological blocks with different positions of the ore in relation with the reservoir boundaries: basal mineralization in contact with the underlying clay formation, perched mineralization close to the top of the reservoir, and finally stacked mineralization with a complex vertical superimposition of roll-fronts. Uranium recovery is fast and complete for the basal mineralization (Fig. 7). The leaching solution, with higher density than reservoir water, is guided to the production wells by the bottom impermeable layer, through the mineralization. On the contrary, perched mineralization leads to slow, incomplete leaching. In this configuration, correct acidification of the ore is complicated by downwards flow due to density gradients between the leaching solution and the reservoir, resulting in incomplete dissolution of the uranium bearing minerals. Also, reservoir volume below the mineralization and the bottom impermeable barrier increases dispersion, so that a fraction of the dissolved uranium remains in the reservoir, and decreases uranium recovery rate. Stacked mineralization display intermediate behavior both in term of dynamic and recovery rate.


Reactive transport modeling can also be used as a powerful tool to assess the environmental footprint of an ISR mining site. As previously stated, the potential issues lie mostly within the aquifer targeted by ISR. Low environmental impact for all surface operations can be achieved by careful use of industrial best practices: clean well-field operation, well and plant decommissioning, removal of networks. Therefore, the remaining leaching solutions within the aquifer can be considered as the main source of chemicals of interest to be treated. Stakeholders are rightfully cautious about operators' capacity to guarantee aquifer quality after the end of exploitation (Mudd 2001; Saunders et al. 2016). Modeling can play a role to help demonstrate that post-exploitation environmental targets can be reached.

The first role of this demonstration is to comply with local regulations and the companies own targets. It can help plan a monitoring and restoration strategy, including an estimation of costs. For instance, Clay (2015) insists on the low amount of published data and the importance to include evolution and migration processes in the definition of rehabilitation targets. Finally, the demonstration is fundamental to obtain the social license, a prerequisite to any new project development. With time scales from 10 years to a century after exploitation and spatial range over 10 km (extension of the deposit), the use of reactive transport simulation is fundamental to quantify the demonstration. The demonstration is based on initial conditions data (baseline) including the inherent variability of roll-front systems, a good understanding and quantification of the physicochemical processes of evolution in the wake of the exploitation, and the identification of the key possible contaminants.

The objective of the remediation strategy revolves around the exploited aquifer, in relation with down-gradient outlets, and the demonstration that a potential water resource is not compromised. Pre-exploitation water conditions are usually poor quality (high contents in metal, uranium and its daughter elements, sometimes high salinity), resulting from the ore genesis mechanisms (Sodov et al. 2016). The physicochemical quality target to be reached (and demonstrated) depends on the local regulation: usually same class of water quality, or concentration range as close to the baseline as possible.

In ISR operations specifically, the buffering capacity of the aquifer is largely reduced by the circulation of the leaching solution. Also, recycling of the solution in the reservoir (with uranium only stripping at the plant) tends to increase the concentration in all the dissolved elements including metals, creating a saline plume. Radioactive daughter products of uranium decay (mostly 226Ra) can be mobilized by the leaching solution and accumulated in the reservoir. The actual hierarchy of chemicals of potential concern depends on the geology and the choice of ISR technology. For instance, alkaline ISR does not create pH singularities, however the high complexation capacity of carbonate uranyl increases the effective solubility of residual uranium. In acidic ISR, the main chemicals of potential concern are the pH plume, sulfate and residual uranium.

The use of reactive transport modeling for such purpose is closer to the type of simulations mentioned in the introduction. However, the industrial context brings several specificities. The operator can start testing remediation strategies in the early stage of development of the project, with a view to optimizing the process while maintaining the targets. Here, pilots can play an interesting role: central in the technical and economic demonstration of the viability of the project, pilots developed prior to the operations can be monitored so that experience on the natural evolution of the system is obtained very early in the mining development. Another opportunity is to use the model in a global optimization perspective. Traditionally, mining operators decouple exploitation technique and remediation. However, particularly in ISR, some industrial choices during exploitation can have important impacts on the remediation costs: notably management of the hydraulic balance and confinement of the saline plume during exploitation. A global model, integrating exploitation and post-exploitation phase allows for a global evaluation, taking into account the specific source term resulting from production choices and their evolution under several remediation strategies.

Geochemical mechanisms involved in post acidic ISR mining

Chemicals of Concern.

According to best environmental assessment practices, the chemicals of potential concerns (COPCs) for an ISR operation are defined according to the pre-mining groundwater conditions, the type of ISR process and the post-mining groundwater quality targets agreed upon by the mining company and the regulators. However, it is possible to assess the changes in chemical composition of the aquifer after mining operations, which will help in the remediation solutions to be deployed. Depending on the buffering properties of the aquifer and the water composition target, one may count on to attenuate the residual solution. Indeed, concentrations of the main chemicals of concern will tend to decrease due to dilution on the one hand, and geochemical reactions on the other hand. These contaminants will undergo, natural attenuation while migrating down gradient with regional groundwater.

Three main chemicals of potential concern are addressed specially in acidic ISR context: sulfate, acidity (pH) and uranium. At the end of the exploitation, residual SO4-enriched solutions will have to be managed. In the mined part of the aquifer, the acidic plume remains in a mostly pH-buffers depleted aquifer. Finally, depending on the industrial optimization targets, residual uranium concentration can still be above the initial baseline. Moreover, acidification of the ore body also induces the increase in other elements concentrations. Concurrently to the dissolution of the uranium minerals, the oxidizing solutions lead to the dissolution of redox sensitive metals (transition metals) and radioactive decay products of uranium (mainly 226Ra). Lastly, dissolution and exchange reactions at the surface of clay minerals participate to the increase of the salinity with Al, Ca, Mg, Na, K concentrations reaching locally the g/L range. The pregnant solutions should therefore be considered as a saline and oxidizing plume.

It is important to pinpoint the distribution of the different chemicals of potential concern according to the aqueous and the solid compartments. Historically, the solid fraction, e.g., neoformed minerals or elements sorbed on the surface of minerals exhibiting high sorption properties, were most of the time disregarded. This retroactively explains why most of remediation operations based on pump & treat strategies failed or were costly as only the aqueous fraction of the chemicals of concern is targeted. Injection of fresh waters tends to dissolve newly-formed minerals such as SO4 bearing minerals (gypsum, alunite, jarosite…) known to be less soluble in the physico-chemical conditions prevailing during the production stage (high SO4 concentrations). The desorption of metals from the surface of clay minerals also depends on the concentrations of other competitive cations which differ drastically between the fresh waters injected during the pump and treat process and the pregnant solutions remaining after the acidification stage. This is illustrated in Figure 8, where the distribution of uranium between solid, aqueous and sorbed species is simulated for a pilot at the end of production. Aqueous uranium represents only less than 20% of the total uranium stock while exchangeable uranium accounts for more than 40%. This stock is therefore available and may be desorbed when fresh waters are renewed during the pump and treat operations. This result clearly highlights the importance of the reactive transport modeling as a tool to help the mining operator in building a remediation strategy.

The geochemical modeling strategy should therefore be seen as assessing the buffering capacity of the aquifer submitted to a saline, oxidizing, and acidic plume. It is important to carefully analyze the main geochemical mechanisms involved in the fate of the different chemicals of potential concern. These data will help in the definition of the remediation strategy, from passive remediation, e.g., natural attenuation, to active remediation operations such as pump and treat or bioremediation (Descostes et al. 2014).

A special attention is put on the development of the geochemical model. Indeed, even if the area of interest is the same as for the production optimization, the key-geochemical reactions involved in the environmental impact can significantly differ. The model previously presented satisfies the constraints of the operator and focuses mainly on the oxidation of uranium and the acid consumption. In other words, geochemical reactions that may be considered as secondary for the production of uranium may become relevant when trying to assess the migration of chemicals of concern at larger time and space scale.

During the exploitation stage, an emphasis is made on the competition between kinetics of major reactions and flow rates in the order of several cubic meters per hour. Meanwhile, when considering longer duration for the environmental post-mining remediation, migration of chemicals of concern is mainly driven by the very low natural flow of the aquifer (~ m/y). Less extreme ranges of pH, compared to production phase, also reduce the impact of kinetics.

One typical example is the consideration of sorption on clay minerals. Mineralogical characterization demonstrated the presence of clay minerals (smectite and kaolinite) within the host formation (Robin et al. 2015a). Their micromorphology and spatial distribution, localized mainly in the porosity of sandstone, allow them to be also involved in the regulation of the water chemistry. Indeed, they exhibit significant cationic sorption properties (Robin et al. 2015b, 2017) and low solubility in acidic conditions (Robin et al. 2016). During the acidification stage, smectite minerals occurring in the porosity sorb protons (Fig. 9); the amount of sorbed protons remains negligible when comparing the total amount of acid used during the lifetime of a technological block (several years). However, once the mining operations cease, proton exchange from sorption sites contributes to maintain moderately acidic conditions, between pH 4-5 depending on the slow renewal of cations from circulating fresh waters. This is well illustrated in Figure 10, where the pH measured in a production well was measured during and after an ISR test. In such case, the acidification was performed for approximately 200 days. The environmental survey performed after the ISR test indicates an increase of pH after several years, but values remain buffered at pH 4. Such trend can only be simulated if sorption is considered on clay minerals (de Boissezon et al. 2017). Finally, sorption, in these physico-chemical conditions, is also a key mechanism governing the mobility of U, 226Ra and other metals (Reinoso-Maset and Ly 2016; Robin et al. 2017).

Another difference between the geochemical models developed for production and remediation, lies in the number of chemical elements to be considered and their respective range of concentration. In the production geochemical model, only major elements are taken into account with of course the reactivity of uranium. Concurrently, in the geochemical model devoted to the assessment of the environmental footprint, even if SO4 and in a lesser extent acidity can be seen obviously as major elements, the fate of residual uranium and other redox sensitive metals, has to be modeled as trace elements, or even ultra-trace elements in the special case of 238U radioactive decay products (mainly 226Ra). Hence, the long term evolution modeling of major elements constitutes a necessary preliminary step before attempting to model the fate of the trace elements.

Lastly, the major difference between the two modeling approaches is the time and space scales of interest. During the production stage, modeling focuses on short durations and on very restricted areas delimited by the localization of the injectors and producing wells (block of production scale). In such approach, local heterogeneities in the permeability distribution or in uraninite distribution may induce significant discrepancies in the predicted production curve. This is why a 3D description of the permeability and the mineralogy is unavoidable. This requirement is less clear for the assessment of the environmental footprint: it depends on the objectives of the simulation, more precisely on the timescale of interest.

The evaluation of the extension of the acidic plume in the vicinity of a technological block during the production phase (e.g., due to hydraulic imbalance) is controlled by flow, which derives from spatial heterogeneity. The environmental piezometer network localized in the vicinity of each technological block allows detecting very shortly such intrusions. For this application, the modeling strategy needs to rely on a 3D approach, with careful analysis of the block model; i.e., production modeling can be used directly in the optimization of the environmental footprint (see Fig. 2).

On the other hand, longer durations and distances are investigated for post-mining simulation, up to the kilometer scale or more. In these conditions, the vertical extension of the aquifer is negligible. Also, the knowledge of geometry outside the exploited area is reduced so that permeability and facies distribution are less constrained. Both reasons lower the need for 3D modeling. Furthermore, although this is not a justification in itself, the geochemical complexity of the system leads to CPU intensive models, so that 2D is easier to handle: e.g., a faster 2D approach can help test different remediation scenarios and assess their efficiency. This 2D/3D modeling strategy is consistent with a far field/near field modeling devoted to different aims and is essential when discussing environmental impact of mining activities with stakeholders (Fig. 11).

Main geochemical mechanisms involved in environmental impact assessment

Sulfate may be considered as the main chemical of concern regarding its chemistry. Remaining aqueous sulfate within the technological block porosity is expected in a first stage to be diluted in the aquifer but also to precipitate in SO4-bearing minerals (gypsum, alunite, jarosite…) in the vicinity of the mined technological blocks.

Precipitation is favored by the general increase of concentrations in cationic species (Ca2+, Al3+ and other cations) inherent in the acidification of the ore body. Indeed, as previously indicated, calcium is released in solution consecutively to the dissolution of calcite but also to exchange reactions on clay minerals according to the overall reactions:


where X≡Ca2+ and X≡H+ arbitrarily stand for the Ca2+ and H+ sorbed at the surface of smectite type minerals (see Robin et al. 2015 and Reinoso-Maset and Ly 2014 for more details). Partial dissolution of alumino-silicate minerals provides sufficient aqueous aluminium to favor the precipitation of alunite type minerals. These reactions are expected to occur shortly after the closure of a technological block. Natural attenuation may be effective mainly thanks to aquifer dispersion, while SO4-bearing minerals precipitation lowers aqueous SO4 concentration during and right after the end of production. On the other hand, high solubility SO4-bearing minerals are slowly dissolved with the renewal of the pore volume by fresh water, therefore maintaining sulfate at relatively high concentration.

Additionally to the dilution of sulfate downstream the geological block, its immobilization is expected according to the prevailing reducing conditions typical of the classical geochemical architecture of a roll front. Taking into account the regional hydraulic gradient, mitigation of the residual solutions will occur mainly downstream, where reducing conditions are found. These conditions are confirmed and explained by the presence of native sulfate-reducing and metal-reducing bacteria involved in the biogeochemical reactions regulating the water composition (Coral et al. 2018). (Bio)reduction of sulfate into sulfide is known to be favorable to the formation of insoluble metallic sulfide minerals such as pyrite (FeS2(s)), therefore participating to the decrease of SO4 concentration.

Acidic conditions inherited from the mining operations will slightly evolve to mid-acidic conditions (around pH 4-5) in a first stage. Such evolution is linked to the slow dissolution of clays minerals and exchange reactions at the surface of clay minerals as previously illustrated (Fig. 10). For longer duration (> several years), as the plume of chemicals of concern migrates downstream, additional dissolution of carbonate minerals will participate to increase the pH and to favor a slow return to pH value close to pristine conditions. This is illustrated in Figure 12, where predictive modeling allowed us to assess the pH potential evolution at the kilometer scale in the vicinity of several technological blocks. Note that the modeling was performed before mining operations started in this area; here, the simulation evaluates the natural evolution towards mid-acidic conditions and locally a hypothetical return to pristine pH values. Validation of such mechanisms was obtained from laboratory experiments and field observations as shown in Figure 10 and in the literature (Yazikov and Zabaznov 2002; Kayukov 2005; Jeuken et al. 2009; Schmitt et al. 2013; Dong et al. 2016).

Residual aqueous uranium is expected to be immobilized by precipitation into newly-formed minerals, stable in geochemical conditions imposed by the ISR process, and by sorption on clay minerals (see dedicated references for sorption above). These reactions are expected to occur shortly after the closure of a technological block as shown in Figure 13. For longer duration (e.g., in the range of years), a reestablishment of reducing conditions is reasonably expected in agreement with the natural redox conditions occurring downstream the ore body as mentioned for the discussion related to the fate of sulfate. These conditions will favor the reduction of uranium and its stabilization under the form of highly insoluble uranium (IV) minerals like uraninite (UO2(s)). The same geochemical mechanisms are proposed to operate for other redox sensitive metals released during the ISR process. Reducing conditions are favorable to immobilize uranium and most metallic elements as their respective solubility is usually 103 to 104 times lower than in oxidizing conditions (Schmitt et al. 2013).

An emphasis is proposed here on the mobility of 226Ra as one of the main radioactive decay products of 238U. According to the amount of dissolved uranium, using a 238U/226Ra ratio of 1, one may expect identical aqueous concentrations, expressed in Bq/L. Field measurements show lower 226Ra activities than expected (see Fig.14), indicating that more than 90% of the 226Ra is trapped, even in acidic conditions. Through its specific activity and its relatively long half-life (1620 y), 226Ra must be considered as ultra-trace element with typical concentrations in the range of ppb to ppt. Therefore, its mobility in natural context is principally governed by sorption at the surface of ancillary minerals, natural organic matter and coprecipitation in SO4-bearing minerals (see among others references dedicated to the mobility of 226Ra in mining context: Lestini et al. 2013; Sajih et al. 2014; Reinoso-Maset and Ly 2016; Robin et al. 2017; Bordelet et al. 2018, Lestini et al. 2019). No proper Ra mineral phase has ever been found in natural environments.

The mobility of 226Ra was successfully simulated in acidic ISR context taking into account three main 226Ra hosting minerals: uraninite, which according to its age may be in secular equilibrium with 226Ra (typically if the mineralization is older than 2 My), barite (BaSO4), and montmorillonite. Such ancillary minerals were observed during and after the acidification stage illustrating their stability (Fig. 15). Also, their dispersion within the pore space allows maximum contact and reaction with circulating fluids. These minerals are known to be efficient traps for 226Ra, especially barite (Curti et al. 2010; Brandt et al. 2015). Even at low concentration, with typical barite grade around 1 to 3 mg/kg, this mineral may contain up to 1011 Bq/L of 226Ra, i.e., several order of magnitude higher than concentrations observed during acidic ISR.

Reactive transport simulations of 226Ra were able to reproduce environmental survey data, only when considering both 226Ra coprecipitation in barite and sorption on clay minerals:


where X ≡ Ra2+ and X ≡ H+ arbitrarily stand for 226Ra2+ and H+ sorbed at the surface of smectite type minerals (see Robin et al. 2017 for more details).

Implications for an industrial remediation solution strategy.

Depending on remediation targets, a dedicated strategy has to be developed by the mining operator. To assess its efficiency and have preliminary associated costs, reactive transport modeling appears also as an interesting complementary tool. Among the remediation strategies, one may distinguish on the one hand passive solutions (usually referred to as natural attenuation) and on the other hand active solutions (pump and treat, bioremediation…).

Natural attenuation relies mainly on the buffering properties of the aquifer to attenuate the oxidative, saline and mid-acidic plume inherited by the acidic ISR mining operations. This remediation strategy was deployed beyond the framework of acidic ISR context, mostly for contaminated aquifers (Christensen et al. 2004, Jørgensen et al. 2010). Among its advantages, one may cite the cost as it relies mainly on an environmental monitoring associated to predictive reactive transport modeling. Natural attenuation in acidic context is however scarcely described even if some examples are available in the scientific literature (Yazikov and Zabaznov 2002; Fyodorov 2002; Bakarzhiyev et al. 2004; IAEA 2005; Dong et al. 2016). Uncertainties on duration and spatial extension appear as the main disadvantages of such remediation strategy. Duration can be shortened using additional pumping operations, which led to the concept of enhanced natural attenuation (Park et al. 2007). In mildly cemented arkosic sandstone where smectite type minerals are found, pH mitigation may drive the duration of natural attenuation as previously discussed. Sulfate appears also as the other chemicals of concern, which governs the overall duration of mitigation, since precipitation of gypsum like minerals is limited in time and space (see further Fig.16): indeed, fairly soluble SO4 minerals do not constitute a stable sink for sulfate. Concurrently to dispersion, reducing conditions usually observed downstream the natural hydraulic gradient of a roll front constitute the major natural mechanisms involved in the attenuation of SO4.

Pump and treat like solutions encompass several technical processes, which are mostly mining site specific. It may be summarized according to several steps. First, residual solutions are pumped out and then treated on dedicated water treatment stations on surface (usually by reverse osmosis). Once the targeted chemical quality is reached, waters are reinjected in the leached part of the aquifer. Additional steps or sequences can be adapted with the injection of a reducing agent to prevent the oxidation of sensitive redox elements, and/or a mixing step with fresh waters (IAEA 2001, 2005). Return of experience and associated reactive transport simulations both highlight the limited efficiency of pump and treat based solutions (Anastasi and Williams 1984; Hall 2009). Indeed, these solutions focus only on the aqueous fraction of the chemicals of concern, which may be in specific case a minor part of the targeted elements (see Fig. 8). To avoid the dissolution of the residual trapped target elements, the remediation strategies should therefore also consider the overall geochemical equilibriums, in the reservoir including the history of production.

All these aspects should benefit from reactive transport modeling, with a special emphasis on the link between production and remediation stages. This is illustrated in Figure 16 where simulation of the environmental footprint was investigated, integrating the initial phase of production of several technological blocks. Simulation was achieved on a 20 × 20 km2 area according to a sequential acidification of the technological blocks. Precipitation of gypsum is predicted on the edges of each technological block, which contributes to limit the spreading of SO4.

Recently, alternative remediation strategies have been proposed. They were developed initially and generally out of the scope of ISR. They rely on the catalysis of natural geochemical reactions to reduce the uncertainty on the duration needed for natural attenuation (Jove Colon et al. 2000). Such approach requires a full understanding of the geochemical equilibriums prevailing before the exploitation of an ISR mine. Amongst the main geochemical processes naturally regulating the aqueous concentrations of major and trace elements in the aquifer before ISR mining operations, microbial activity plays an essential role in establishing the contrasting redox conditions through the roll front. This is especially the case for the reducing conditions observed within the mineralized zone and downstream, where sulfato-reducing bacteria are present (Coral et al. 2018). Taking into account the low solubility of uranium and other redox sensitive metals in reducing conditions, a return to reducing conditions after mining operations can be targeted using relevant bacteria. Hence, depending on the use of native bacteria or exogenous bacteria, reducing conditions can potentially be reached in shorter durations than the ones related to natural attenuation. Concurrently, acidity and SO4 are also targeted. Such alternative remediation solutions fall into a bioremediation strategy with on the one hand, biostimulation and on the other hand, bioaugmentation. Frequently used for organic compounds, several industrial tests have been launched on uranium-polluted aquifers and (conventional) uranium mining sites (Wu et al. 2006a,b; Groudev et al. 2008, 2010; Yabusaki et al. 2014; Long et al. 2015). These bioremediation techniques are still under study, mostly in laboratory or small-scale pilots, with no industrial scale demonstration on ISR sites. Once more, reactive transport modeling can be used in the understanding and interpretation of bioremediation tests from the laboratory scale, in column experiments for instance, up to the field test (see for instance Yabuzaki et al. 2007).

Industrial use of reactive transport modeling for post-mining.

All these examples can be summarized in a general remediation strategy adapted to acidic ISR where reactive transport modeling is a key tool to assess the efficiency of the retained remediation solution. The specific remediation strategy (among others natural attenuation, pump and treat and bioremediation) depends on the local geology, mining technology, and environmental targets defined by local regulation.

  1. Prior to the closure of any geological block, a predictive modeling of the remediation solution is provided. Such model should merge, both the production and remediation stages in order to take into account the production history. Economical aspects are considered according to the design of the environmental survey network and the duration of monitoring.

  2. Deployment of a dedicated environmental survey, and monitoring of the efficiency of the remediation solution over a reasonable period of time to evaluate if predictive modeling agrees with observations.

  3. If there is a good consistency between the modeling and the environmental monitoring then maintain the environmental survey with an adapted chronic (decreasing frequency).

  4. If no agreement between modeling and observations is found, then the model should be updated (including using new data from the monitoring phase), another monitoring phase should be launched, and eventually the remediation solution should be adapted.


Several applications above illustrated how an academic, versatile, reactive transport code like HYTEC can be used in an industrial context. The fast pace imposed by operations, need for quick results, imposed the use of a specific graphical interface, specially devised for this application: HYSR. The code is presently tested, in operation conditions, on the Katco ISR mine.

An important lesson from this project is the benefit of co-construction. After the initial identification of an industrial problem and the possibility to address it with reactive transport tools, a close cooperation was built between an academic team, with the knowledge of the tool and expertise to optimize it, and an industrial team, with the knowledge of industrial context and issues and the driver to improve their simulation capacity. Also, access to an extensive mass of data was key in the early stage of the project to build, calibrate, and optimize the model. As the project moved from R&D department to technical department and finally operations, it became clear the code had to evolve from a versatile, keyword-based input files to a more intuitive interface. GUI is of course essential, but a key was to calibrate the interface to the needs of the engineers: the GUI is basically devised so that ISR engineers interact with the code using the frame of mind of their trade: inputs organize various data of the project and proposed exploitations scenarios, outputs are displayed following the need of production engineers including their financial conversion for comparison purpose.


This work was funded by the industrial chair “in situ recovery of uranium”, a program supported by the French National Research Agency and Orano. The authors also wish to thank Katco and the Mongolian project team, for fruitful discussions and access to years of operation data. Finally, the authors wish to thank many people involved in the simulation project, particularly Gwenaële Petit, Hélène de Boissezon, Tatiana Okhulkova, and Marie Mazurier, and more generally in the ISR project, Valérie Langlais, Nicolas Fiet and Anthony Le Beux.

Open access: Article available to all readers online.