Antoine Collet
email: antoine.collet@orano.group

Olivier Regnault
email: olivier.regnault@orano.group

Alexandr Ozhogin

Assemgul Imantayeva

Loïc Garnier
email: loic.garnier@orano.group

Three-dimensional reactive transport simulation of Uranium in

Keywords: In situ recovery, Uranium, Reactive transport, History matching, Optimisation, HYTEC

come

Introduction

In situ recovery (ISR), also known as in situ leaching (ISL) is currently the leading uranium extraction technique. From sic either-or oxidising solutions, increasing uranium mobility within the permeable deposit, which is pumped out towards the processing plant. Before being reinjected, the leaching solutions are restored by the additions of some chemical, e.g., acid (IAEA, 2016).

The absence of mechanical extraction makes ISR by far the most cost-effective extraction technique [START_REF] Kidd | Uranium mining -what method works best[END_REF], and potentially more environmentally friendly [START_REF] Taylor | Review of Environmental Impacts of the Acid In-situ Leach Uranium Mining Process[END_REF].

Consequently, it is perfectly suited for large, deep and lowgrade sedimentary orebodies, where traditional open-pit and underground mining techniques are uneconomical [START_REF] Dahlkamp | Uranium Deposits of the World: Asia[END_REF][START_REF] Kyser | 13.19 -Uranium Ore Deposits[END_REF]. The highly permeable and well confined roll-front deposits found in the USA, Uzbekistan and Kazakhstan are the archetypes of these deposits and currently dominate the global production (OECD-NEA & IAEA, 2020). Although this technique is characteristic of uranium mining, re-cent studies highlighted a high potential for its use in extracting other commodities such as copper or gold [START_REF] Sinclair | In situ leaching of copper: Challenges and future prospects[END_REF][START_REF] Seredkin | In situ recovery, an alternative to conventional methods of mining: Exploration, resource estimation, environmental issues, project evaluation and economics[END_REF][START_REF] Filippov | The BIOMOre Project -A New Mining Concept for Extracting Metals from Deep Ore Deposits Using Biotechnology[END_REF][START_REF] Vargas | In situ recovery of copper sulfide ores: Alternative process schemes for bioleaching application[END_REF]. Hence ISR might be a fast-growing technol-30 ogy in the coming decades.

The economic potential and the technical challenges associated with ISR recently increased research interest. In particular, extensive literature has exposed the potential of reactive transport (RT) to investigate ISR [START_REF] Nguyen | Numerical simulation of uranium in-situ mining[END_REF][START_REF] Nos | Modèle Conceptuel d'une Exploitation d'uranium Par Récupération in Situ[END_REF] 35 [START_REF] Martens | In situ recovery of gold: Column leaching experiments and reactive transport modeling[END_REF][START_REF] Johnson | Reactive transport modeling at uranium in situ recovery sites: Uncertainties in uranium sorption on iron hydroxides I[END_REF][START_REF] Ben Simon | Kinetic reactive transport modelling of column tests for uranium In Situ Recovery (ISR) mining[END_REF][START_REF] Johnson | Predictive Reactive Transport Modeling at a Proposed Uranium In Situ Recovery Site with a General Data Collection Guide[END_REF][START_REF] Laurent | Numerical modelling of column exper-1135 iments to investigate in-situ bioleaching as an alternative mining technology[END_REF]. However, it is usually restrained to 1D laboratory experiments, such as column tests. In addition, the majority of ISR operations are conducted in post-Soviet countries (e.g., Kazakhstan), and 40 a large number of publications are rarely accessible and published in Russian. Hence, only a few recent papers illustrate the industrial usage of RT modelling in operations [START_REF] Regnault | 3D Reactive Transport simulations of Uranium In Situ Leaching : Forecast and Process Optimization[END_REF][START_REF] Noskov | Innovative Intellec-1165 tual Management Technology of Uranium Mining by the ISL Method[END_REF][START_REF] Lagneau | Predictive simulation and optimization of uranium in situ recovery using 3D reactive transport simulation at the block scale[END_REF][START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF]. In 2019, [START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF] revealed the first results of the operational deployment of RT simulations at the Katco mine in Kazakhstan using the HYTEC RT simulator (van der [START_REF] Van Der Lee | Module-oriented modeling of reactive transport with HYTEC[END_REF]. They showed that an RT simulation based workflow reproduced the production history of one block (61 wells), a process known as history matching. They assessed the robustness of the model by applying it to another block in the vicinity, yielding equally good results. This paper documents the robustness of the same model by presenting the results obtained from 39 production blocks (2394 wells) in two independent reservoirs. This article represents a significant step-forward, placing RT history matching on the scale of what has been achieved in the neighbouring industry of oil and gas: industrial applications of nonreactive transport simulation coupled with history matching to drive the operations, allowing to increase the amount of oil and gas that can be extracted from existing resources [START_REF] Oliver | Recent progress on reservoir history matching: A review[END_REF][START_REF] Rwechungura | Advanced History Match-1200 ing Techniques Reviewed[END_REF][START_REF] Udy | Review of Field Development Optimization of Waterflooding, EOR, and Well Placement Focusing on History Matching and Optimization Algorithms[END_REF].

In addition, the literature mentioned the main applications of a robust RT-based workflow as short-term planning [START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF], long-term planning [START_REF] Langanay | Uncertainty quantification for uranium production in mining exploitation by In Situ Recovery[END_REF] and environmental permitting through impact monitoring and postoperation remediation [START_REF] Coral | Microbial communities associated with uranium insitu recovery mining process are related to acid mine drainage assemblages[END_REF][START_REF] De Boissezon | Modeling uranium and 226Ra mobility during and after an acidic in situ recovery test (Dulaan Uul, Mongolia)[END_REF]. The present paper extends these results by presenting a real short-term well field optimisation case study that illustrates the geometallurgical approach and the economic benefits possible when managing the end of life of production areas with an RT-based workflow. This paper is organised as described below. After providing a brief description of the Katco uranium mine, the first section presents the ISR modelling process with HYTEC. It includes the description of the workflow inputs (3D geological models, geochemical model, and operational conditions), and the presentation of modelled areas. Section 3 describes the history matching and case study results. Finally, Section 4 discusses limitations and perspectives on the robustness of the workflow and the history matching process, the importance of operational conditions in forecasts, and finally, the industrial capacity of the workflow.

Materials and methods

Katco Mine

Created in 1996, Katco JV LLP is a joint venture mining 85 company between Kazatomprom (49%) and Orano Mining (51%), which are the Kazakh and French national nuclear fuel cycle companies, respectively. With a yearly production ranging from 3200 t to 4000 t since 2009, corresponding to 7% of the global uranium production in 2019 (OECD-NEA & IAEA, 2020), Katco 90 remains the largest in situ recovery (ISR) uranium mine worldwide [START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF]. The mine, located in the Shu Saryssu Bassin -Kazakhstan, exploits the Tortkuduk and Muyunkum deposits along a southwest-northeast axis over 40 km in length.

It currently operates this large territory in an acidic way thanks 95 to processing plants, each of which is connected to the deposit by a well field.

Well field description

A well field is a large network of wells, that connect the processing plant to the orebody. Figure 1 schematically illus-100 trates the organisation of one of the well fields sustaining the ISR process. The wells are arranged in hexagonal patterns with a radius of 40 m on average at Katco. At the hexagon vertices, injection wells (injectors) introduce the leaching solution to the orebody. Production wells (producers) are placed at the centre 105 of the hexagons and deliver the pregnant solution to the processing plant. This pattern, referred to as "seven spots", ensures good circulation of the solutions through mineralisation and hence results in more effective leaching [START_REF] Ward | Well Design and Construction for in Situ Leach Uranium Extraction[END_REF][START_REF] Shao-Chih | Well-Field Mechanics for In-Situ Mining[END_REF]. Each hexagon is called a production cell and is 110 considered a unit of extraction covering 4150 m 2 of the deposit.

Due to the reservoir thickness, vertical overlapping up to three levels of wells is frequently required to intersect mineralisation.

(Figure 1). In addition, the wells are frequently redrilled when the initial well is clogged or to locally accelerate the leaching 115 process.

The hydraulic network connecting the plant to the wells (Figure 1) is organised in three subsets. First, a primary network of pipelines, named primary lines, links the plant to pump houses called transit chambers. Then, a second level of pipelines, called secondary lines, links these chambers to well houses referred to as acidification units. Finally, the acidification units connect the secondary lines to each well individually through a third network of buried pipes. An acidification unit is usually linked to a dozen production cells composed of approximately 125 45 injectors and 15 producers. This group of cells is called a production block. The mine currently operates 80 production blocks at any given time.

Production management

Four key indicators are used to manage the production and 130 optimise the reservoir response. On the one hand, flow rates (injected and pumped) and added acid volumes are the two levers the operators use to adjust the production. On the other hand, the pH and uranium concentrations of pumped solutions allow us to calculate uranium production and acid consumption, from 135 which the economic performance of the mine is derived. Therefore, accurate and regular knowledge of these four process data is critical.

Each well and line are equipped with a flowmeter to monitor the flow rates, whose cumulative values are read every 24 h.

140

The pregnant solutions from the producers and production lines are sampled and analysed on the same daily basis, providing the uranium concentration and pH. A systematic correction workflow is applied to make the measured data consistent with the basic conservation principles. First, outliers are removed using 145 a combination of Hampel [START_REF] Hampel | The Influence Curve and Its Role in Robust Estimation[END_REF][START_REF] Liu | On-line outlier detection and data cleaning[END_REF]) and other median-based filters (see [START_REF] Gonzalo | Weighted Median Filters[END_REF]. Then, missing values and mass conservation are addressed using a dynamic mass balance reconciliation method performed at the acidification unit scale. The process exploits the measure redundancy 150 between wells and lines and is based on an unscented recursive nonlinear dynamic data reconciliation scheme [START_REF] Vachhani | Robust and reliable estimation via Unscented Recursive Nonlinear Dynamic Data Reconciliation[END_REF] which is an adaptation of unscented Kalman filters [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF]. This workflow is complemented by a systematic visual quality assurance before using the data.

155

Regarding composition of the pregnant leach solution, the transit of the pumped solutions through the processing plant before being reinjected induces a homogenisation of the concentrations [START_REF] Nos | Modèle Conceptuel d'une Exploitation d'uranium Par Récupération in Situ[END_REF]. Consequently, these concentrations are sampled at the outlet of the plant and considered identical 160 for the entire well field connected to the plant. Sampling was performed twice a day and the following major species were measured: Al 3+ , Ca 2+ , Cl -, H 2 SO 4 , Fe 2+ , Fe 3+ , K + , Mg 2+ , PO 3 - 4 , SO 2 - 4 , UO 2+ 2 . The final acidity of injected solutions is calculated using the baseline concentration from the plant, 165 complemented by the added acid in the acidification unit (one of the operational levers previously discussed).

Reactive transport workflow for ISR at Katco Mine

The RT simulations [START_REF] Lichtner | Reactive Transport in Porous Media[END_REF][START_REF] Steefel | Reactive transport modeling: 1230 An essential tool and a new research approach for the Earth sciences[END_REF][START_REF] Druhan | Reactive Transport in Natural and Engineered Systems[END_REF] et al., 2003). The HYTEC program provides a general mul-175 tipurpose flexible framework for solving complex hydrogeochemical problems [START_REF] Lagneau | HYTEC results of the MoMas reactive transport benchmark[END_REF][START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF]. The code was successfully employed in benchmark studies [START_REF] Carrayrou | Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case[END_REF] and various applications, e.g., cement degradation [START_REF] De Windt | Modeling the degradation of Portland cement pastes by biogenic organic acids[END_REF]Seigneur 180 et al., 2020), and geological storage of acid gases [START_REF] Sin | [END_REF]. The workflow developed using HYTEC for the ISR simulation has already been extensively described by [START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF]. Consequently only major elements of the workflow relevant to this study are presented below i.e., input 185 parameters and modelled areas. The designed RT-based workflow requires three types of input parameters [START_REF] Regnault | 3D Reactive Transport simulations of Uranium In Situ Leaching : Forecast and Process Optimization[END_REF][START_REF] Langanay | Uncertainty quantification for uranium production in mining exploitation by In Situ Recovery[END_REF]. First, a 3D geological block model, which contains the hydrodynamic parameters and the mineralogical description of the reservoir. Second, a geochem- (Olea, 1999, p.7-30).

3D geological block models

All 3D geological models used for the simulations were built from the same type of input data, and possess the same structure, to ensure a consistent representation of the Muyunkum sandstone roll-front [START_REF] Petrov | Epigenetic stratified-infiltration uranium deposits of Kazakhstan[END_REF]Dahlkamp, 2009, 241-253).

A 3D geological model consists of the reservoir discretised into a 3D grid of rectangular parallelepipedal voxels. The grid mesh size used (x, y, z) was (5 m, 5 m, 1 m) or (10 m, 10 m, 1 m). This grid held three geological properties. First, a two-class discrete lithological property discriminates the permeable (i.e., sandstone) from the impermeable facies (i.e., massive clays and carbonate layers). In HYTEC, these facies were associated with homogenous permeabilities of 1 × 10 -4 m/s and 1 × 10 -10 m/s for permeable and impermeable facies, respectively. For porosity, a unique value was used for each reservoir: 23% and 18% for Uyuk and Kanjugan, respectively (Figure 4). Second, the roll-front is represented by a discrete variable of three redox classes referred as OMR facies [START_REF] Petit | Application of Stochastic Simulations and Quantifying Uncertainties in the Drilling of Roll Front Uranium Deposits[END_REF][START_REF] Langanay | Uncertainty quantification for uranium production in mining exploitation by In Situ Recovery[END_REF] 2a). A homogenous dedicated mineralogy is defined for each facies. Third, the uraninite and coffinite accumulations are represented as a scalar field (Figure 2b). The input data available to build the models are the core description, core analysis (X-ray fluorescence) and drill hole logging (resistivity and gamma ray).

Although consistent input data were used to build all 3D geological models, the number of drill holes available to constrain the model varied depending on the area. Indeed, for each area of the mine, two datasets were available from the different drilling campaigns. The first available set corresponds to the exploration and development drill holes, which are required to characterise the deposit and to build the well field design (production cells and screen positioning). The usual drill spacing 230 observed is 100 m by 50 m ( 200 -250 drill holes per km 2 ).

The second is obtained from the drilling of technological wells, prior to the start of extraction. The 40 m diameter hexagon cell design provides much denser coverage (> 650 drill holes per km 2 ). Consequently, the geological model constraints used for a model mainly depend on the moment the model was built, in regard to the available data at that time. Figure 3 illustrates the spatial distribution of the drill holes used to build the block model covering the primary line A. The distribution of the ED wells is not homogenous and the technological well data included cover blocks A1 and A2. Furthermore, in this article, the models built only with exploration and development drill hole data are referred to as ED while those built on top of the technological wells are referred to as EDT.

In addition to the number of available drill holes, the modelling process differed as two geostatistical methods were employed. The first is a stochastic based method and relies on the one hand on plurigaussian simulations for the discrete variables [START_REF] Fontaine | Simulation of the Muyumkum uranium roll 1080 front deposit by using truncated plurigaussian method[END_REF][START_REF] Langlais | In the shade of the truncated gaussian simulation[END_REF][START_REF] Renard | 3D representations of a uranium roll-front deposit[END_REF], and, the other, relies on block simulation for the uranium grades scalar field. The models were built using ISATIS software combined with the R-geostats package (MINES ParisTech / ARMINES, 2021), based on the workflow described by [START_REF] Petit | Application of Stochastic Simulations and Quantifying Uncertainties in the Drilling of Roll Front Uranium Deposits[END_REF]. The second approach is deterministic and relies on simple kriging of indicators for the facies 255 and on punctual ordinary kriging for the uranium grades (Olea, 1999, 7-30). In that case, the models were built within the PE-TREL software (Schlumberger, 2018;[START_REF] Suslov | Allocation of paleochannels in sand-clay sediments of middle Jurassic age in petrel software complex[END_REF]. The main difference between the two approaches is that the first method preserves the geological variability while the 260 second method tends to smooth it.

The effects of drill hole density and geostatistical methods are investigated later in this study.

Geometry and operational conditions of the well field

The geometry and operational conditions used for the simu-265 lations were provided at the well scale. For each well, geometry was defined using 2 parameters: coordinates of the well head and the vertical position of the well screen. The flow rates and composition of leaching solutions were defined on a daily basis using the data mentioned in Section 2.1.2.

Geochemical model

The geochemical model employed is derived from previous studies that identified the reactive mineralogy of uraniumbearing sandstone in the Shu-Saryssu basin [START_REF] Yefteyeva | Research and development of methods for intensifying the process of in situ leaching of uranium from ores of the Chu-Saryssu depression (on the example of the Uvanas deposit[END_REF][START_REF] Petrov | Epigenetic stratified-infiltration uranium deposits of Kazakhstan[END_REF][START_REF] Ben Simon | Kinetic reactive transport modelling of column tests for uranium In Situ Recovery (ISR) mining[END_REF][START_REF] Robin | Occurrence of authigenic beidellite in the Eocene transitional sandy sediments of the Chu-Saryssu basin (South-Central Kazakhstan)[END_REF][START_REF] Robin | Dissolution of beidellite in acidic solutions: Ion exchange reactions and effect of crystal chemistry on smectite reactivity[END_REF].

275

The kinetic and mineral databases as well as the underlying chemical processes used for the simulations were described previously (Lagneau et al., 2019, and references therein), and are briefly recapitulated below.

The geochemical model relies on six minerals: cristobalite 280 alpha, kaolinite, smectites (beidellite), calcite, iron hydroxides (goethite), and uraninite.

On the one hand, the acidification of the aquifer is described in a straightforward manner through H 2 SO 4 injection and reactivity with calcite, beidellite and goethite using equations 1, 2 285 and 3, respectively [START_REF] Robin | Dissolution of beidellite in acidic solutions: Ion exchange reactions and effect of crystal chemistry on smectite reactivity[END_REF][START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF]. (2)

CaCO 3 (s) + 2H + ---⇀ ↽ ---Ca 2+ + CO 2 (aq) + H 2 O (1)
FeO(OH) + 3 H + ---⇀ ↽ ---Fe 3+ + 2 H 2 O (3)
On the other hand, the leaching of the U(IV) mineral phase (uraninite) is controlled by oxidative dissolution (Eq.4) by Fe 3+ , which is first derived from the reactions with gangue minerals described previously, as well as recirculation of Fe 3+ from the 290 well field. Gangue minerals play a major role since they both control the pH and provide the oxidiser.

UO 2 (s) + 2Fe 3+ ---⇀ ↽ ---UO 2+ 2 + 2Fe 2+ (4)

Modelled areas

As shown in Table 1, history matching work was performed on 39 production blocks from ten different primary lines, including 538 production cells composed of 2394 wells (1744 injectors and 650 producers). The primary lines are labelled History matching is the process of building one or more sets of models that account for the observed reservoir behaviour. It is a type of inverse problem which consists of an attempt to adjust some of the reservoir model parameters to better reproduce the state variables of interest observed in the past, taking into account the real operational conditions imposed on the reservoir. The greater the agreement among the forecasts, the more predictable the model [START_REF] Oliver | Recent progress on reservoir history matching: A review[END_REF]. In practice, we choose to adjust parameters for which we have a limited knowledge, and for which the state variables display sufficient sensitivity [START_REF] Yeh | Review of Parameter Identification Procedures in Groundwater Hydrology: The Inverse Problem[END_REF][START_REF] Carrera | Inverse problem in hydrogeology[END_REF].

Adjusted parameters and methodology

Three geochemical parameters were adjusted to fit the observed data: initial grades of beidellite (a single value for the three redox facies), grades of iron hydroxides in the oxidised facies and, finally, grades of calcite in the mineralised facies.

The choice of these three parameters was motivated both by internal studies that showed their spatial variability, and by the limited data available. All other parameters remained constant among the simulations and were described by Lagneau et al.

(2019) (as detailed in Section 2.2.3).

Parameters adjustment followed an iterative process. For

RC(γ, γ) = γ -γ |γ| LRD(γ, γ) = ln γ γ (5)
where γ and γ are the simulated and the observed data, respectively. In both cases a positive relative change indicates that γ > γ, i.e., the simulation overestimates γ and vice versa. The reason for using the LRD rather than another indicator (i.e.

the RC) is that it is a symmetric, additive and normed relative change metric [START_REF] Taylor | Review of Environmental Impacts of the Acid In-situ Leach Uranium Mining Process[END_REF]. Indeed, the RC is bounded to the interval [-1, +∞[, where γ being null or positive, and possibly infinitely superior to γ. With RC, the negative and positive deviations would not be treated in the same man-375 ner, e.g., a cell with an observed production of 10 t and a 5 t simulated cell has an RC of -50% while the opposite, a cell with an observed production of 5 t and a 10 t simulated cell has an RC of +100%. In the case of the LRD, however, the previous example produces symmetric values of -30.1% and +30.1% respectively.

In addition to the LRD, the MASE has been used to capture the instantaneous deviation and score the fit quality over the full simulation. Among the many measures proposed in the literature for evaluating the forecasting performances of different 385 methods [START_REF] Hyndman | Another look at measures of forecast accuracy[END_REF][START_REF] Kim | A new metric of absolute percentage error for intermittent demand forecasts[END_REF], the MASE, introduced in 2006 by Hyndman and Koehler, has proven to possess good mathematical properties, defined mean, finite variance and scale independence [START_REF] Franses | A note on the Mean Absolute Scaled Error[END_REF]. The MASE is now a well-established metric [START_REF] Makridakis | The M4 Competition: 100,000 time series and 61 forecasting methods[END_REF] as it overcomes most of the bias occurring when comparing method performances on different datasets. It is defined as the mean absolute error (MAE) for the predicted values, scaled by the MAE of a simple one-step "naive forecast method", which consists of shifting the values one step backward (one step equals one day in our case):

MASE = 1 n n i=1 |γ i -γ i | 1 n-1 n i=2 |γ i -γ i-1 | (6)
where γ i -γ i is the simulation error for the sample i, and n is the number of samples in the time series. This metric should be understood as follows: a MASE of 0 equals a perfect fit; a MASE of 1 equals a fit as good as the 400 one-step "naive forecast method"; and a MASE of n, n times worse, e.g MASE of 10 indicates a 10 times higher MAE. The U and acid concentration evolution being continuous and rather smooth, as shown in the results section (e.g., Figure 5), a MASE should be visually appreciated as an almost perfect match in the 405 interval [0, 2], as very good in ]2, 3], as good in ]3, 5], acceptable in ]5, 7] and poor beyond.

End of block life management

The capacity to predict future reservoir performance is a key aspect of history matching and constitutes the true usefulness 410 of an RT-based workflow. The ability to forecast the reservoir response to operational conditions with increased confidence indeed allows operational experiments to be performed to optimise the well field management. The end of block life management is one of these optimisation applications.

Results

History matching

Uranium production

Figures 5 and6 show the observed uranium production curves against the simulation results for the eight primary lines (see 2.2.4) and four blocks of interest respectively. Interestingly, all primary lines and blocks presented in this paper produce a large 450 range of uranium extraction curves. This is due to several factors. First, the blocks of a primary line are usually exploited sequentially, e.g., for primary line E in Figure 5, each peak corresponds to a block start. This plays a major role in the production profile at the line scale but also at the block one, scales of interest: primary line, block and cell. Since a single block was modelled for both primary lines J and I (Table 1), their results were considered at the block scale and not at the primary scale. The LRD of the total uranium produced follows a quasi-normal distribution approximately centred at zero for all 500 scales, which implies that the underestimations and overestimations globally compensate. Additionally, the smaller the scale is, the higher the distance of the distribution mean and median from zero and the higher the standard deviation of the misfit ( 2.3±14.0%, 8.6±28.7%, 25.8±69.9% for primary lines, blocks good. It increases to 5.8 ± 5.2 at the block scale and more than doubles (9.7 ± 9.3) at the cell scale with a standard deviation 4.5 times higher than at the primary line scale.

Acid consumption

Figure 8 shows the history matching results for the pumped solution pH and the corresponding back calculated acid recovery for the blocks G2, E3, H2 and H5. For uranium, and for the same reasons (see Section 3.1.1), the profiles vary substantially from one block to another. However, unlike uranium, the acidity fit results in alternate intervals when the data are correctly matched and intervals with important LRDs. This fit is particularly acute during the first year of production (80.1±62.1%) because the model relies on the acidic dissolution of the smectites but ignores the ion-exchange reactions that play a major role in acid consumption [START_REF] Robin | Dissolution of beidellite in acidic solutions: Ion exchange reactions and effect of crystal chemistry on smectite reactivity[END_REF][START_REF] Robin | Fate of dioctahedral smectites in uranium roll front deposits exploited by acidic In Situ Recovery (ISR) solutions[END_REF]. This choice is mainly motivated by technical reasons, as incorporating sorption processes would be very CPU-intensive.

The results on the 39 blocks show that the intervals can be classified into three classes (Figure 8): The mean MASE is also 1.21 times larger at the primary line scale than at the block scale because the blocks of a primary line operate sequentially. As acidification is conducted at the block scale, the primary line scale aggregates intervals with potentially different behaviours. Consequently, separation between behaviours cannot be achieved at the primary line scale.

590

Interestingly, the deviations are more centred at zero and exhibit a lower standard deviation for acid consumption than for uranium, both at the block scale (3.8 ± 29.1% vs. 8.6 ± 28.7%) and at the well scale (-2.7 ± 46.6% vs. 25.8 ± 69.9%). Additionally, the average mean MASE is 2.23 times lower at the 595 block scale and 2.15 times lower at the well scale with a lower standard deviation. This result suggests a better performance of the model for predicting acid consumption than for predicting uranium output. This better performance is easily explained by the fact that the beidellite content is manually adjusted by the user, while the uraninite distribution is the result of the geostatistical process with high uncertainties. It also confirms that the local variations in beidellite content are much lower than those in uraninite content. Regarding uranium production, ED assets show a more offcentre average relative log deviation with a higher standard de-615 viation at every scale, and the MASE is systematically higher for ED with a larger standard deviation (more extreme over and underestimations). The accuracy of the uranium recovery fore-cast model is directly linked to the constraint density used to build the geological model. Therefore, the fact that a lower con-620 straint on the 3D block model results in higher misfits proves that a major part of the misfit is due to the incorrectness of the 3D block model, and not the RT-based workflow itself.

For the acid recovered through the pumped solutions, LRD and MASE are only slightly better when the well-density con-625 straint is low (ED) which might seem contradictory. This result is, however, explained by the direct control of the acid consumption by the beidellite grade, which is not a geostatistical outcome, but imposed by the user.

Comparing the stochastic and the deterministic methods for 630 EDT, the first method produces better results at the primary line and block scale for uranium recovery and slightly better results at the cell scale. Therefore, this method is recommended for systematic use in all studies on operating blocks, i.e., already drilled (EDT).

Optimisation case study -Primary line A end of life

As introduced in Section 2.4, the true usefulness of history matching and of a robust model lies in the capacity to predict the response of a reservoir, with increased confidence, while considering various operational conditions. The next section presents the results of the optimised end of life of primary line A, based on the history matching results previously presented.

First, the optimisation results from 2019 are described in detail.

Second, the predictions made for the optimised well-field based on the history-matched model are compared to 16 months of observed production data.

Optimisation results

The history matching results showed that 37.1% of the initial estimated uranium reserves were still present under the outline of the technological blocks (day 3248 in Figure 12) before the optimisation of primary line A. Following the implementation of the procedure described in Section 2.4, the optimised design was successively improved.

Figure 13 shows the overall operational conditions tested and the production predictions for four different scenarios. Only the first two years of predictions were considered because the uncertainties in the operational parameters were considered too strong after this period. The first of these scenarios, called the reference case, relied on a well field design established empirically. In that case, the cells retained for the end of life were selected based on the previously observed production curves (high uranium concentration) and low estimated recovery (back calculated from the initial reserve of the cell minus the observed production) recorded as is and simulated in HYTEC. The three other scenarios were built based on the procedure using the RTbased workflow, which is described in detail in Section 2.4.2.

The "Int. 1" and "Int. 2" cases represent intermediate test scenarios where additional wells have been proposed. The final scenario corresponds to the optimum scenario reached in the limited time available for the study.

The operational conditions used to generate the predictions vary very little from one scenario to another. For example, the cumulative flow and the total acid tonnage varied by 8% and 6%, respectively. The final case shows the lowest acid usage and a water imbalance. This imbalance is due to injectorproducer inversions, temporary injection through the producer to fully invade its vicinity and lower the pH, which were not 695 considered in the previous cases.

In contrast to operational conditions, production predictions increase substantially as scenarios are optimised (P base < P int1-2 < P f inal , with P the uranium production). Between the base case and the final scenario retained, HYTEC predicted a 28% in- The results show that the predicted and observed production 735 trends match very well. After 484 days, or 16 months, predictions only underestimate total uranium production by 9.2%.

Using actual operating conditions the trends are even better than expected, and the underestimation drops to 1.6%, which is excellent. First, this result once again proves the robustness and accuracy of the model to reproduce the reservoir response.

It is indeed capable of reproducing a complex production history spanning more than 3500 days, including operational constraints and well field design that vary substantially over time, which no analytical approach is capable of reproducing. Sec-

Discussion

Model robustness and history matching

The results presented in Section 3 have shown the capacity 755 of the workflow to correctly reproduce production histories, and to predict the response of the reservoir while considering the operational constraints imposed on it. This robustness is strong at the scale of the production block and primary lines, which makes the tool usable for all mine planning operations and con-760 stitutes a major breakthrough in the exploitation of uranium by ISR.

However, manual calibration performed globally on the initial grades of calcite, iron hydroxides and clays, is not sufficient to explain all the variability in production data, particularly at However, in roll-front deposits, two in-situ chemical elements are gamma emitters, namely, uranium (U) and radium (Ra), and are not evenly spatially distributed [START_REF] Boulesteix | Ilmenites and their alteration products, sinkholes for uranium and radium in roll-front deposits after the example of South Tortkuduk (Kazakhstan)[END_REF]. This results in Ra/U ratios varying sometimes up to an order of magnitude, with complex 795 distributions. The few available measurements increase the difficulty of mapping these disequilibria, and in practice, a single correction ratio value is used across several primary lines.

It logically has the effect of providing a good estimate of the overall grade but taints the local values with error. This explanation is completely consistent with the results obtained in Section 3. 1.1, i.e., errors at the cell scale follow a normal distribution centred at zero, and compensate for each other, resulting in better matches at block and primary line scales.

In the absence of additional data (drill holes) to constrain the geostatistical processes (which represent unaffordable costs for the operator), the remaining alternative to improve 3D geological models is the automation of the parameter adjustment by a computer: mineral (e.g., uraninite and beidellite) grade scalar fields, as well as hydrological parameters such as per-810 meability and porosity. This approach, known as automated history matching, has been widely studied and applied to reservoirs since the 1960s (de [START_REF] De Marsily | Four decades of inverse problems in hydrogeology[END_REF], in hydrogeology [START_REF] Carrera | Inverse problem in hydrogeology[END_REF][START_REF] Zhou | Inverse methods in hydrogeology: Evolution and recent trends[END_REF], and in the petroleum industry [START_REF] Oliver | Recent progress on reservoir history matching: A review[END_REF][START_REF] Rwechungura | Advanced History Match-1200 ing Techniques Reviewed[END_REF].

815

Automated history matching allows us to adjust scalar or discrete fields up to the voxel scale, improving the forecast of models with higher confidences. Therefore, the development of large-scale reactive-transport inversion frameworks is highly desirable, and it is probably the area with the highest potential 820 for workflow improvements.

Operational conditions

Although a robust workflow and an efficient history matching process are key aspects that help solve the problem of predicting future reservoir performance, a good estimate of future

RAM used

Fig. 15. Runtime, number of voxels in the 3D blocks, number of days simulated and RAM used in the history matching simulations presented.

quently an important research topic, which should be developed alongside automated history matching.

Industrial capacity

One of the fundamental aspects of a numerical model in- ) is very important and should not be neglected while adding complexity to the workflow. An interesting approach to overcome these speed and memory issues is the usage of surrogate modelling, i.e., a statistical model of the model, also called a meta-model or proxy model [START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF]Wang and 885 Shan, 2006). This approach has already been widely used in the oil and gas industry at industrial scale [START_REF] Feraille | Prediction under Uncertainty on a Mature Field[END_REF][START_REF] Feraille | An Optimization Strategy Based on the Maximization of 1070 Matching-Targets' Probability for Unevaluated Results[END_REF]Costa et al., 2014, among others), with limitations especially for automated history matching [START_REF] Zubarev | Pros and Cons of Applying Proxy-models as a Substitute for Full Reservoir Simulations[END_REF].

Comments about alternative simulation approaches

As explained in the section 2.1.2, the operator monitors the 920 well field behaviour watching three main parameters: flow in the field, dissolved uranium concentration at the block scale, and acid consumption. [START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF] illustrated the limit of traditional ISR operations forecast, where evolution is evaluated independently for the three parameters:

925

• Flow in the field is ignored.

• 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 solu-930 tion is based on analytical functions using parameters adjusted on production data: the extrapolation of these functions is used as a predictor.

Such approach is acceptable only if operating conditions remain constant which is not the case in practice as illustrated with 935 primary lines and block production curves (Figures 5 and6).

Another alternative is pure hydrological simulation which has been reported for well field optimisation [START_REF] Krause | Hydrological characterization and optimization of in-situ recovery[END_REF] and focuses on the flow within the reservoir, taking into account hydrological parameters only. In the end, only the 3D reactive transport approach gathers 950 the three pillars of ISR operations, which are (i) reservoir geometry and composition (ii) well field geometry (iii) operational conditions. This is the only method that correctly considers the leaching within the reservoir and provides a real metal balance and dissolution and recovery rates. Instead, an accurate quantitative evaluation is needed and can-990 not be provided by empirical approaches. Second, large datasets of concentrations, flow rates, and operational conditions are available and must be considered when making these decisions.

While qualitative approaches are still possible, numerical models such as HYTEC are required to fill in the gaps and make 995 full use of data. Third, the time is ripe for geometallurgical approaches, when economic and environmental aspects are the main drivers in operation [START_REF] Dominy | 1065 Geometallurgy -A Route to More Resilient Mine Operations[END_REF][START_REF] Michaux | How to Set up and Develop a Geometallurgical Program FINAL[END_REF].

The workflow presented here fulfils the needs by coupling geology, geostatistics and metallurgical response of the reservoir.

Fig. 1 .

 1 Fig. 1. Schematic view of the Katco Uranium In situ recovery (ISR) mine and the ISR process.

  allow us to model and 170 understand the coupled hydrodynamic and geochemical processes that occur within the aquifer during ISR operations. All RT simulations presented in this study were performed using the code HYTEC developed at Mines ParisTech (van der Lee

Figure 1 ,

 1 Figure 1, O=yellow, M=red and R=grey.

Fig. 2 .

 2 Fig. 2. Examples of 3D block models.

Fig. 3 .

 3 Fig. 3. Spatial distribution of the drill holes used to build the 3D model of primary line A. The red dots are the exploration and development drill holes (ED) while the blue dots are the drill holes added after the drilling campaign in blocks A1 and A2 (mostly technological wells). The light-grey polygon shows the total extent of the primary line A, while the dark grey polygon shows the extent of blocks A1 and A2.

  270

Fig. 4 .

 4 Fig. 4. Location of the simulated technological blocks. The primary line I, located north, is composed of four isolated blocks that are not displayed on the map.

340

  the simulation of a new block, the grades used for the closest calibrated block were established as a starting point and used in a first run. The results were compared to the observed data at the block scale, and the pH was adjusted by manually increasing or decreasing the beidellite grade, repeating the process un-345 til a satisfactory fit was obtained. The same principle was then applied to the calcite grade to fit the initial decrease in pH (∼ 30 first days) and finally to the iron hydroxide grade affecting the rate of uranium dissolution during the first year.2.3.3. Fit accuracy metrics 350As explained by[START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF], the necessity for the operator to produce reliable mining plans in both short and long terms prompts interest in block performance projections at different temporal and spatial scales. However, a long-term mining plan requires a good overall estimation of uranium produc-355 tion and acid consumption at the primary line scale i.e., cumulatively over long times; short-term planning demands a correct estimation on a monthly to weekly basis and at the block or cell scale. A good production model is therefore expected to provide both. Thus, the qualification of predictive models involves 360 the adaptation of metrics allowing the intercomparison of different variables (acid and uranium concentrations), observed at various temporal (yearly and daily) and spatial (primary line, block, and well) scales. The two metrics applied in this paper are the log relative deviation (LRD) and the mean absolute 365 scaled error (MASE). The LRD has been used to assess the final deviation of the cumulated uranium production and acid consumption, as an alternative to the classic relative change (RC). They are defined in equation 5 below:

  415 2.4.1. Concept As a production area ages, uranium reserves are gradually depleted and therefore the uranium concentration of pumped solutions globally decreases. Below a certain concentration cut-off, production cells are systematically stopped, until the 420 production block itself is closed. The spatial distribution of uranium is very uneven, and while many cells are depleted, some cells nevertheless still offer high potential and can be exploited. Managing the end of block life consists, on the one hand, of determining which areas must remain active and, on the other 425 hand, of optimising the design with well re-drillings. 2.4.2. Primary line A end of life In 2019, Collet et al. introduced an HYTEC-based optimisation methodology for the final stage of uranium recovery in old production blocks. This paper presents its implementation 430 in a real case study based on primary line A. The methodology follows an iterative process. As a first step, a careful history matching is performed to adjust the model parameters (see Section 2.3). This step reinforces the predictivity of the model and provides access to the location of the remaining uranium in the 435 deposit. The location of uranium at the end of extraction determines the choice of production cells to retain for the next step. As a second step, a forecast simulation estimates the metallurgical response of each cell, which does not depend solely on uranium grades. Based on these data, the design is sequentially optimised by well re-drillings and different acidification strategies. The final outcome is the implementation of the most profitable design and acidification strategies in the well field.

  455as the blocks that started later benefit from the acidification of the first blocks to be operated. Second, as studied by[START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF] and[START_REF] Langanay | Uncertainty quantification for uranium production in mining exploitation by In Situ Recovery[END_REF], the local distribution of uraninite and the geometry of the reservoir (Figure4) vary from one area to another and it has a great impact on arrival and 460 heigh of production peaks, and shapes of descending slopes.Third, these diverse production profiles are explained by the modification of the operational conditions with time, i.e., the acidity of the leaching solutions and the flow rates. Although some blocks such as E2 exhibit an intense initial acidification 465 followed by a constant decrease, other assets, as illustrated by blocks B1 and F1, undergo multiple acidification campaigns, sometimes carried out with cell additions and closures. In addition, some areas, particularly for the Kanjugan deposit, are prone to clogging, requiring many workovers and cleaning actions to maintain a sufficient flow rate. This phenomenon is characterised by sawtooth wave production profiles, each peak matching a well cleaning campaign. It is clearly observed at primary line H in Figure 5, and specifically at block H5 in Figure 6. 475 Despite the very different behaviours of blocks and primary lines, the model systematically manages to reproduce the trend and is usually very close to the production data, for both the instantaneous and cumulative extraction. This result indicates the robustness of the workflow to account for the geological and 480 geochemical spatial variability and its capacity to accurately predict the response of the reservoir with respect to the various operational conditions imposed by the operator. According to Lagneau et al. (2019), the misfit between the model and the observed data is largely due to the uncertainties 485 in the geology and, more specifically, in the uranium distribution. This explanation has been sustained by the study of the effects of the geological uncertainties on the simulated uranium recovery performed by Langanay et al. (2021). Additionally, Lagneau et al. (2019) argue that the block size is well above 490 the range of constraints for the geostatistical model, and the uncertainty is globally averaged at the block scale. The results of the statistical study of the misfit shown in Figure 7 validate these assumptions. The univariate statistics of the LRD and MASE for uranium and acid production examined three spatial 495

505Fig. 5 .Fig. 6 .

 56 Fig.5. History matching results of uranium production at the primary line scale. For visualisation, the daily values are smoothed by median filtering with a 3-day window. LRD is the logarithmic relative deviation, and MASE is the mean absolute scaled error.

530•Fig. 7 .Fig. 8 .

 78 Fig. 7. Log relative difference (fact vs. simulated) and mean absolute scaled error (MASE) of the uranium and acid production predictions per scale of interest. The boxes show the interquartile range (Q3-Q1) of the dataset while the whiskers extend to show the p5-p95 of the distribution.The red squares represent the mean values. µ is the average, m is the median, σ is the standard deviation and n is the number of samples.

3. 1

 1 Figures 9a -10b, show results of the statistical analyses (LRD and MASE), respectively, of the uranium and acid deviations with respect to the 3D block model building method (stochastic vs. deterministic) and to the density of wells used 610

  Spatial variability of adjusted geochemical parametersAs indicated in Section 2.3.2, the adjustment necessary to obtain the results previously presented only included three global parameters: initial grades of beidellite (a single value for the three redox facies), grades of iron hydroxides in the oxidised 640 facies and, finally, grades of calcite in the mineralised facies.

Figure 11

 11 Figure11shows the univariate distribution and spatial variability (from one block to another) of these parameters. Two remarkable observations were made. First, the parameters are spatially continuous i.e., they slightly vary between two neigh-645

Fig. 9 .

 9 Fig.9. Same as for Figure7but for uranium production only, grouped by 3D geological model modelling method and well constraint density.

Fig. 10 .

 10 Fig. 10. Same as for Figure 9 but for the acid recovered.

  Calcite grades (wt%).

  Iron hydroxides (goethite) grades (wt%).

  Distribution (box-and-whisker plot) of the parameters among the blocks.

Fig. 11 .

 11 Fig. 11. Spatial and univariate distributions of the geochemical parameters used to calibrate the models (adjusted at the technological block scale). In d) the boxes show the interquartile range (Q3-Q1) of the dataset while the whiskers extend to show the p5-p95 of the distribution. The red dashed lines and squares are the median and the mean of the datasets, respectively; µ is the average, m is the median, σ is the standard deviation and n is the number of samples.

700

  crease in uranium production over the first two years. Application of the cost model, which includes CAPEX (drilling, acidification unit,. . . ) and OPEX (electricity, acid, plant reagents,. . . ), as well as incomes, results in an increase in the net present value (NPV) of 35%. In addition, the final overall uranium concen-705 tration of the pumped solutions is increased by 23%, indicating that the trend will continue and that the gap widens beyond the two considered years.

Figure 12 Fig. 13 .

 1213 Figure12shows the spatial distribution of uranium and the pH plume at key time points in a simulation. The simulation 710

Fig. 14 .

 14 Fig. 14. Uranium production predictions from 2019 compared with observed production (March 2021) and predictions updated with observed operational conditions (March 2021).

765

  the technological cell scale. Our assumption is that a correct reproduction at this scale would require a much finer description of the 3D geological model, i.e., spatial variability of the parameters, a finer history matching (up to voxel scale adjustment), and the consideration of more parameters, including the uranium grades. The need for a finer description is supported by two main reasons.First, internal studies have shown that initial grades of calcite, iron hydroxides and clays, as well as hydrodynamic parameters (permeability and porosity), might vary strongly and sharply even at the cell scale. This variability is also observed vertically at drill hole level. Not surprisingly, for instance, the use of an average value of beidellite at the block scale fails to adequately reproduce local acid consumption at the cell scale.Second, the uranium grades, which are one of the main parameters controlling uranium production, have not been adjusted at all in the model calibration process, which in itself is a strength of the current workflow. Uranium grades are indeed only the result of geostatistical processes. However, in addition to the uncertainties inherent in estimation methods and to the uranium 785 local variability below the range of drill hole spacing, at least one major known source of bias affects the evaluation of the spatial distribution of uraninite. This bias is because the estimate of the grade at the boreholes, which is used to constrain the 3D geological model, is calculated mainly from a logging gamma signal (radioactivity).

825

  production also requires a good forecast of the operational conditions imposed to the reservoir. On the one hand, the geometry of the well field and the amounts of injected acid are the result of human decisions, and can therefore be determined precisely through sustained planning. The flow rates, on the other 830 hand, are much more difficult to predict because they are largely controlled by coupled chemo-mechanical clogging phenomena[START_REF] Vergnaud | Mécanismes de colmatage dans une exploitation d'uranium par récupération in situ[END_REF] and the operator can only perform curative actions whose effectiveness is temporary. Consequently, two blocks with equivalent reserves and geometries may have very 835 variable production profiles depending on their level of clogging. Understanding and ultimately predicting these phenomena is essential to improve predictions. The wrong forecast of flow rates is responsible for most of the error in the production forecast at the block scale for the case study presented in Sec-840 tion 3.2.2 (9.2% versus 1.6% for forecasted and observed operational conditions, respectively). This result was even stronger at the cell scale. Accumulated over blocks, this error might have an important impact on the reliability of mining plans, and might affect the result of local optimisations. This area is conse-

  tegrated in a regular workflow of operations is the speed with which the operator can obtain results on which to base his or her decisions. Regarding the presented workflow, due to its robustness and the low variability of adjusted geochemical parameters from one block to a neighbouring block, the model calibration only required two to three runs on average. This calibration is not a limitation for a systematic use of the workflow for mine planning or optimisation studies (several hundred runs per year).However, the calculation time and CPU capacity may become factors limiting several lines of research previously discussed. Figure15shows the runtime, the number of voxels of the 3D models, the number of days of production simulated, and finally the RAM used to generate the simulation results (Section 3). Although the number of production days simulated is generally high -∼ 2250 days on average, as the study focused on areas reaching their end of life, an average of 4 days, and 1000 Go of RAM were required to run a simulation. If these constraints do not prohibit a single run, history matching might need hundreds to thousands of calls of the model with increased memory needs[START_REF] Oliver | Recent progress on reservoir history matching: A review[END_REF].Moreover, the tendency is to increase the size of the grid and the period of time covered. First, simulating multiple blocks at once takes into account their interactions and removes the bias for border cells. Second, environmental monitoring and footprint assessment consider acid migration with the natural flow 875 of the aquifer during the rehabilitation, up to a hundred years after mine closure, and for the far field domain, extending beyond the contour of the production blocks[START_REF] De Boissezon | Remediation of a uranium acidic In Situ Recovery mine: From laboratory to field experiments using reactive transport modelling[END_REF][START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF][START_REF] De Boissezon | Modeling uranium and 226Ra mobility during and after an acidic in situ recovery test (Dulaan Uul, Mongolia)[END_REF]. Consequently, code performance optimisation (memory usage, parallelisation, 880 etc.

  to other deposits and commoditiesStarted in 2006[START_REF] Nos | Modèle Conceptuel d'une Exploitation d'uranium Par Récupération in Situ[END_REF], the project at the origin of the presented workflow represents 15 years of effort and development. These developments involved firstly the adaptation of reactive transport codes to the ISR case in porous media, sec-895 ondly, the creation of a geochemical model, and, finally, the creation of pre-and post-processing interfaces, which actually constitutes the main part of the work[START_REF] Lagneau | Industrial Deployment of Reactive Transport Simulation: An Application to Uranium In situ Recovery[END_REF].Our experience shows that the use of this workflow for a new roll-front uranium mine, such as found in Kazakhstan, Mon-900 golia or Uzbekistan, represents relatively little effort insofar as the geochemical model and the operating scheme do not change significantly. However, it requires much more effort when the geology changes, e.g., for high-grade uranium deposits of Athabasca Basin in Canada. In the latter case, the min-905 eralogy and uranium proportions greatly differ, which means a new geochemical model must be established. In addition, the extraction is performed in fractured media, which means the flow transport must rely on different modelling approaches such as double porosity that must be incorporated in reactive-910 transport tools if not existing. Similarly, mining another metal (e.g., copper) requires creating and testing a new geochemical model and potentially adapting the tools if the geological context greatly differs, which can take up to several years. Therefore, the creation of an ISR modelling workflow for a mine site 915 greatly depends on the existence or not of reactive transport model and pre-and post processing infrastructure for a similar deposit.

940

  However, the use of this approach is limited in operational conditions because (i) homogeneous global values in porosity and permeability are used, data acquisition being too costly (ii) the validation of the hydrodynamic parameters can only be done by means of tracer tests which are spotty, or by having the well-945 head pressures, which is not a standard practice in the ISR industry because of its cost. Moreover, elements such as uranium or acidity cannot be used as tracers as they are largely impacted by the chemical dimension.

  This paper presents one of the first large-scale well field applications of reactive transport (RT) simulation in the context of uranium in situ recovery (ISR). The RT-based workflow was successfully used on 39 production blocks containing 2394 960 wells of the Katco ISR Mine well field. The model was found to be very robust at the block and primary line scales, accurately reproducing uranium production and acid consumption over time scales of up to 12 years, with dramatic time-varying operational conditions (well field design, flow rates, and leach-965 ing solution acidity). The robustness of the model is based on the adjustment of only three global parameters of the 3D geological models: initial grades in clays (beidellite), initial grades in calcite, and initial grades in iron hydroxides (goethite). In addition, these calibration parameters are spatially continuous and consistent with mineralogical internal studies. The statistical study of the model deviations has shown that the errors at the production well scale compensate for each other, therefore indicating that most of the remaining forecast errors are attributable to the 3D geological model and not to reactive transport simulations. The optimisation case study validated the predictive nature of the RT-based workflow on a real example of an area end of life management, showing less than 10% deviation between predictions and observed data 16 months later. It also revealed 980 the significance of gains provided by the workflow in operations: +28% increase in uranium production and +35% increase in economic gains just over the first two years of end of block life management.Overall, this paper proved the capacity of an RT-based work-985 flow to be used industrially and how it became a game changer on several accounts. First, hydrogeologists and mining engineers are increasingly expected to make difficult decisions for which qualitative understanding of reservoirs is not sufficient.

  is the culmination of a fifteen-year research effort by scientists from Orano Mining and Centre de Géosciences Mines ParisTech. It also represents 20 years of careful data collection and preparation by operators and en-1005 gineers from Katco JV LLP. Therefore, this study should be recognised as collective work, and all people involved should be acknowledged. Additionally, the authors would like to thank Gwenaële Petit, R&D engineer at Orano Mining for her constant technical support, and Dr. Irina Sin, Pr. Hervé Chauris 1010 from Centre de Géosciences Mines ParisTech, and Dr. Valérie Langlais from Orano Mining for their valuable advice, comments and support while preparing this manuscript. The authors also would like to thank assistant editor Sadegh Safarzadeh, editor-in-chief Gamini Senanayake, and three anonymous re-1015 viewers for their valuable comments and suggestions that helped to improve the manuscript. Finally, KATCO JV LLP is thanked for providing its permission to share and present the results of this work.

Table 1

 1 Characteristics of the modelled primary lines.

	Size