Mechanistic models supporting uncertainty quantification of water quality predictions in heterogeneous mining waste rocks: a review

Polluted drainage from weathering of sulfide-rich waste rock deposits can cause long-term impairment to waterways and biodiversity near mining sites. Mechanistic models represent established tools to support the predictions of the quantity and quality of waste rock drainage, and their associated risks. Yet, model-based predictions in typical waste rock systems are ubiquitously uncertain because of the strongly heterogeneous nature of these waste deposits. Embedding heterogeneity within predictive modeling is complicated by the magnitude and level of knowledge of the waste rock heterogeneity, and the large number of scale-dependent parameters feeding the model equations. This review encompasses deterministic and stochastic modeling approaches that emphasize consolidated tools and emerging modeling solutions to deal with heterogeneity for the modeling of waste rocks. Physical (e.g., variability of texture, hydraulic and pneumatic properties), geochemical (e.g., variability of mineralogy and kinetic parameters), and thermal heterogeneities are evaluated. The review points out the importance of stochastic modeling as a fundamental approach to embed uncertainty in long-term model-based decisions. Regulators and decision makers must be convinced of the benefit of using stochastic modeling, which is still considered to belong mainly to the academic sphere.


Introduction
Mine waste management represents one of the biggest environmental and socio-economic concerns for mining operators, governments and citizens. A critical concern is related to the weathering of sulfide-rich mining byproducts such as waste rocks and tailings, which can generate a toxic leachate enriched in metals and metalloids. This drainage is very often called acid rock drainage (ARD). The formation, fate and management of ARD, including passive and active remediation approaches, has been addressed in multiple reviews (Blowes et al. 2003;Dold 2017;Muniruzzaman et al. 2018;Nordstrom et al. 2015;Nordstrom and Alpers 1999;Plumlee 1999;Wolkersdorfer et al. 2020;Wolkersdorfer and Bowell 2004), guidelines (EPA 1994;European Commission 2009;INAP-GARD 2014;MEND 1998;Price 2009;Sobek et al. 1978) and books (Blowes et al. 2014;Lottermoser 2010).
The deleterious impact of ARD on water resources, and in turn on public health and aquatic ecosystems, is well known and studied. From a health risk perspective, ARDcontaminated water supplies, such as pumping wells (Christenson 1995), can expose the population to elevated concentrations of metal(loid)s, such as lead (Pb), arsenic (As), manganese (Mn), and cadmium (Cd). Excess of metal(loids) has been linked to multiple physiological alterations, e.g. neurological alteration, kidney disease, chronic obstructive pulmonary disease, hypertension and heart disease (Neuberger et al. 2009;Wright et al. 2006). The transport of ARD to water resources affects aquatic ecosystems by a number of direct and indirect pathways, which include habitat modification, niche loss, bioaccumulation within the food chain, loss of food sources, elimination of sensitive species, reduction in primary productivity, and food chain modification (Betrie et al. 2016;Gray 1997).
Waste rocks are among the most common forms of mine waste. Amos et al. (2015) extensively reviewed the physical and mineralogical characteristics of waste rock piles, providing a list of the principal processes related to sulfide oxidation and solute loading. They concluded with a discussion on acid mine drainage prediction and prevention techniques, which included mechanistic model based analyses. Amos et al. (2015) indicated in particular that ''numerical models will provide a valuable tool to incorporating greater complexity and mechanistic process description into leach water quality predictions''.
Model-based predictions of ARD from waste rocks are severely complicated by the presence of spatially variable physical and biogeochemical properties. A conceptualization of the type of heterogeneities that this paper tackles is shown in Fig. 1: (a) physical heterogeneity controls the movement of water and associated solute dynamics in ARD-generating environments; (b) biological and geochemical heterogeneity affect the reactions involving gas, water and minerals, with biota as the catalyzers of geochemical processes, during the formation and evolution of ARD plumes; (c) thermal and pneumatic heterogeneity have strong influence on more specific aspects, such as mechanisms sustaining the supply of oxygen towards oxidizing sites. Note that the ARD-receiving environment could also be a heterogeneous system, such as an alluvial aquifer (Fig. 1) or a fractured aquifer.
Heterogeneity generates uncertainty in the parameters feeding mechanistic models used to make predictions of ARD from waste rock. Consequently, these model-based predictions are bound to be uncertain in a system without adequately resolving the magnitude of the involved heterogeneity. While stochastic modeling is an appropriate framework to deal with heterogeneity-driven uncertainty in hydrogeological problems (Freeze 2004;Guadagnini and Fig. 1 Schematic diagram illustrating the generation of acid rock drainage and the subsequent propagation of leachate from waste rocks (as piles) to a receiving environment (groundwater aquifer and surface water bodies). In the overburden pile, darker colour represents the compacted fine-grained materials whereas the light colour represents coarse waste rock material. The inset illustrates heterogeneity associated with the water-gas-solid phases as well as bacterial species at pore-scale. The drainage from the pile percolates to the aquifer, which can be also a heterogeneous system, with spatially variable hydraulic properties, such as the hydraulic conductivity (K) Tartakovsky 2010;Lichtner and Tartakovsky 2003;Rubin et al. 1994;Serre et al. 2003), the use of stochastic models for waste rocks modeling is still limited. From the perspective of mine waste simulations, even the deterministic models embedding heterogeneity has been also rarely developed. A main reason is that ARD from waste rocks is controlled by multidimensional multiphase coupled processes described by strongly nonlinear partial differential equations. Intense computational efforts must be considered when resolving heterogeneous waste rock simulations. The development of multicore workstations, parallel computing and access to supercomputers is rapidly changing this scenario, particularly regarding stochastic modeling.
The purpose of this work is to fill a gap in the current literature by presenting a review of the evolution of the approaches that can be used to improve waste-rock modeling considering heterogeneity as a complicating factor for the model predictions. In the first part of the paper, we address an overview of the mechanistic approaches that have been adopted so far for long-term predictions of ARD from heterogeneous waste rocks. This includes both deterministic and stochastic models. Secondly, we present our perspective regarding how the discipline could evolve to increasingly include heterogeneity in future ARD studies. We emphasize the importance of embedding stochastic heterogeneity for risk assessment and other type of studies that require an explicit simulation of heterogeneity as a driver of model uncertainty.
With this review, we aim at providing the reader with a better understanding of model-based tools that can be used to quantify, control or reduce the heterogeneity-driven uncertainty regarding the predictions of the deleterious impact of ARD on water resources. The review targets researchers and technicians, particularly hydrogeological and geochemical modelers, as well as mining practitioners and operators. A goal of this work is to illustrate that the use of heterogeneous modeling is no longer prohibitive, as it could have been just one decade ago: nowadays, modeling is facilitated by easier-than-ever-before accessibility to powerful workstations, super-computers and other solutions for parallel computing, which require simple notions of programming skills to be pre and post-processed. In this sense, heterogeneity-driven uncertainty can be explicitly quantified by modelers and the results presented to decision makers in the form of probabilistic results.
The review assumes that the reader is already familiar with basic elements and processes influencing ARD in waste rocks as well as basic corollary notions associated to mine waste management, such as depositional methods to create waste-rock deposits. Specific information on ARD formation and attenuation mechanisms can be found in a wealth of scientific articles, books and reviews-e.g. about pyrite oxidation mechanisms (Battistel et al. 2019;Evangelou and Zhang 1995;Rimstidt and Vaughan 2003), about the formation of iron and aluminum hydroxysulfates (Bigham and Nordstrom 2000), about mineralogical assemblage (Parbhakar-Fox and Lottermoser 2015), about methods to remediate ARD (Akcil and Koldas 2006;Johnson and Hallberg 2005;Sheoran and Sheoran 2006;Tripathy 2014) or about biologically-related processes (Baker and Banfield 2003;Nordstrom and Southam 1997;Papirio et al. 2013). Ultimately, it is assumed that the reader also possesses notions of flow and reactive transport modeling in porous and fractured media. Steefel et al. (2014) provides a comprehensive review of reactive transport modeling in subsurface environments.

Evolution of mechanistic modeling for heterogeneous waste rocks
This section revises some of the key mechanistic (i.e. process-based) models that have deterministically incorporated heterogeneity in their formulations or been used to compute ARD from waste rocks using heterogeneous parametrizations. While not all the revised models have been used for quantifying the heterogeneity-driven uncertainty in ARD predictions, these models could be used for this purpose. Indeed, most of them provide the basis for rigorous uncertainty quantification analysis, as addressed in the following Sect. 3.

First models
Influential initial contributions were the breakthrough studies by ; Ritchie (1987, 1986). They presented a set of model formulations of pyritic oxidation in mine wastes that incorporated the assumptions of a previous simpler model (Ritchie 1977) and coupled the effect of transport of oxygen into the particles comprising the waste rock deposit with transport through the pore space of the waste rocks. While their study was based on diffusive gas transport only, they first presented ) the general formulations of the model and a simplifying assumption which made it possible to resolve the problem analytically. The approximated analytical model was then tested numerically ). As the previous models were based on single-size particle radii, in a follow-up study (Davis and Ritchie 1987) they extended the model to evaluate a range of particle sizes in the waste rocks, which is a more realistic description of true waste rock piles. A central equation is the 1-D diffusion equation for gaseous oxygen in the waste rocks, with the form: where x is the spatial coordinate, t is time, / is the porosity of the waste rock, C pores is the oxygen concentration within the pore space, D g is the diffusion coefficient of gaseous oxygen in pore space, and q is a source-sink term that account for the mass transfer between the pore spaces and the reactive minerals, with the form: where a is the size of the particles, with range ½a 0 a n ,f a ð Þ is the density of particles between a and da, C part is the oxygen concentration within the reactive particle, and D w is the diffusion coefficient of aqueous oxygen in water. These expressions are then included into the general reactive transport model  to compute sulfide leaching rate and heat generation from the exothermic oxidation of sulfides. The model, called ''distributed particle size model (DPM)'', is one of the first documented heterogeneous reactive transport models for waste rocks.
In Davis and Ritchie (1987), the DPM was used to simulate weathering of an idealized waste rock system with parameterization inspired on the White's waste rock piles at Rum Jungle in the Northern Territory of Australia ). The pyritic waste material was characterized by a pyrite content of 3.3%, a porosity of 0.4, an oxygen diffusion coefficient through the pore space of 6:72 Â 10 À6 m 2 s -1 , and depth of 18 m. A distribution of particles spanning between 1 and 200 mm was adopted, with a larger number of particles corresponding to the larger sizes and with half the mass composed of particles with radii less than approximately 32 mm. The DPM results were compared against the results of two equivalent homogeneous models embedding respectively particle sized of 10 and 50 mm. Davis and Ritchie (1987) found that the heterogeneous DPM generated a drastic difference in the model results compared to the homogeneous models. While the latter could still be accurate for predictive purposes of sulfate production rates at specific time scales and for the tested conditions, predicted heat source distribution at 26 years was sensibly impacted by the selected particles distributions. The authors concluded that it was ''essential'' to include the particle size distribution (i.e. heterogeneity) in the model.
Since Davis and Ritchie (1987), several other mechanistic modeling formulations and analyses have followed. The vast majority of research on ARD modeling from waste rocks has focused on numerical modeling, with a few notable exceptions where analytical solutions were obtained (Binning et al. 2007). The rapid development of numerical models was promoted on the one hand by the increasingly high power and low costs of computational machines, and on the other hand by the increasing awareness of the importance of embedding heterogeneities for model predictions (heterogeneity is hardly embedded into analytical solutions). Nowadays, multiphase multicomponent reactive transport models are able to simulate a wide range of processes, including both advective and diffusive gas transport, the effects of oxidizing agents in contact with sulfidic materials and the interaction between minerals and pore-fluids (da Silva et al. 2009;Ma et al. 2019;Mayer et al. 2015;Molson et al. 2005;Muniruzzaman et al. 2020;Sracek et al. 2004;Wang et al. 2019).

Models embedding physical heterogeneity
Flow and transport heterogeneity in waste rock deposits are primarily controlled by the inhomogeneous spatial distribution of textural properties controlling water infiltration and seepage through the subsurface, and the associated advective and dispersive mechanism of solutes. The degree of heterogeneity depends on natural heterogeneity of the background of the deposits from which the mine waste has been generated and the additional physical heterogeneity created by the artificial alternation of the ore deposits, for instance through blasting, milling, dumping and layering (Bailey et al. 2013;Day 1994;Smith et al. 1995;Smith and Beckie 2003). For instance, end-dumping deposition (Blackmore et al. 2014;Peterson 2014;St-Arnault et al. 2020;Vriens et al. 2019) results in material segregation with coarser material falling for larger distances than the fine particles, which remain on top of the dumping sequence or pile. The sequencing of dumping stages determine a layering structure with compacted traffic surface embedded among the dumps lifts (Lefebvre et al. 2001a, b).
The selection of the most-suitable modeling approach to reproduce and predict flow and transport heterogeneity in waste rocks is not a simple choice and it is still a debated issue. A critical initial element to consider is the importance of capillary and non-capillary flow. So far, mechanistic models for flow in waste rocks have mostly relied on the concept of equivalent porous media, in which continuum approaches based on the Darcy's equation are valid and characterized by properties such as the hydraulic conductivity (K) and the porosity (/). Since waste rocks are unsaturated systems, the majority of waste rock models presented so far have been based on single-porosity capillary-based Richards (1931) equation, in which the key parameterization is based on soil-water characteristic curves (SWCCs) and relative permeability (HF) functions. In general, this approach is valid when the waste rocks have a ''soil-like'' regime, where flow is generally homogeneously rinsing the waste rocks, resulting in a matrixcontrolled flow pattern. An increasing number of studies have found that flow in coarse-textured waste rocks may be still dominated by capillarity and matrix dominated water flow and transport (Appels et al. 2018;Neuner et al. 2013).
The importance of combining textural heterogeneity and hydrological regimes to properly understand flow in waste rocks is well known. Newman et al. (1997), for instance, developed a 2-D deterministic two-column model embedding two distinct materials, a fine-grained soil and a coarse layer soil. Adopting a finite element code based on a single-porosity version of the Richards equation, Newman et al. (1997) showed that the type of material where infiltrating water preferentially travels through is strongly controlled not only by the soil texture but also by the rainfall regime and the selected SWCCs and HFs for the two materials. Under heavy rainfall, the influx rate at the model top boundary exceed the infiltration capacity of the fine-grained material and water enters preferentially the coarse-grained materials. Under light rainfall, the material enters preferentially the fine-grained material which has (under unsaturated conditions) a higher relative permeability than the coarse-grained material.
Strong contrasts in textures can result in more pronounced preferential channels embedded into a less permeable matrix. Multiple studies (e.g. Blackmore et al. 2018Blackmore et al. , 2014 have suggested that the presence of non-capillary-driven flow can significantly complicate the processes conducting recharging through the pore space similar to what occurs for instance in structured soils (Beven and Germann 2013;Simunek et al. 2003). When preferential flow dominates the waste rock's seepage regime, only a few zones can collect majority of the infiltrating water within the deposits (Nichol et al. 2005). Morin and Hutt (2008) suggested that as little as 5% of the entire mine deposit may be actually rinsed by the infiltrating water. In the Diavik experimental waste rocks, Neuner et al. (2013) reported that the porosity (/) of the bulk waste rock in the test piles is estimated to be / = 0.24, the void or macropore porosity is estimated to be / = 0.18, while the pore space affected by capillarity is approximately the remaining of the bulk volume of the test piles.
The heterogeneity in the flow distribution affects directly advective and dispersive transport mechanisms of solutes (as well as reactions, as described in detail in the next sections) in waste rocks. Internal and external tracer tests have been widely used to elucidate the implication of preferential flow in heterogeneous waste rocks, and the type of models to be selected. The presence of ''dead pores'' largely conditions the fate of a conservative tracer within saturated waste rock (Murr 1979). Wetting front velocities in coarse heterogeneous material travel several orders of magnitude faster than the mean water velocity in granular waste rocks (Nichol et al. 2005;Webb et al. 2008). Internal tracers, in particular blasting residuals, have been also adopted to characterize heterogeneous mine waste deposits (Bailey et al. 2013;Blackmore et al. 2018;Hendry et al. 2018;Mahmood et al. 2017). At the Antamina waste rock experimental site, Blackmore et al. (2018) found that, while a waste rock pile characterized by coarse-grained materials shows rapid flow and long tailing in the estimated breakthrough curves (BTCs) for the externally injected and internal tracers, the finer grained pile also shows clear tailing on both tracers' BTCs. At the Aitik mine in Sweden,  estimated that flow and transport could have accounted preferentially for about 55-70% of the total water content in the local waste rocks.
Modeling solute transport in heterogeneous waste rocks can be done in multiple forms. One option is through multidimensional models that explicitly simulate the variability of textural properties of the waste rocks. Molson et al. (2005) presented one of the first examples in this sense by coupling the models presented by Fala et al. (2005) to POLYMIN (Molson et al. 2004), a reactive mass transport and sulfide oxidation model based on the wellknown multicomponent geochemical model (Allison et al. 1991) (this code is also described in the next section for what concerns its reactive modeling capabilities). They showed that solute transport tends to channel along the preferential flow created by flow instabilities. In real site, these instabilities, correspond to flow channels, would occur naturally due to intrinsic heterogeneities in the materials and along the interface and also due to pore-scale variations in the moisture content. Another possible modeling approach is based on an effective (or lumped) parameterization of the 1-D advection-dispersion equation, which can account for multiple porosities, each one characterized by a specific water mobility (Simunek et al. 2003). Blackmore et al. (2018) used the dual-porosity version of HYDRUS-1D (Simunek et al. 2005) to fit multiple tracer tests in two waste rock experimental piles in Antamina concluding that a possible conceptual model for transport in waste rocks would consist of at least three porosity components: (1) a porosity associated with fast-flow channels, activated under strong recharge events; (2) a porosity associated with a slower-flow matrix material, activated or controlled by natural rainfall events; (3) a no-flow porosity, i.e. the immobile domain, in which transport occurs through diffusive processes only.

Models embedding geochemical heterogeneity
Polluted drainage from mining waste deposits involves acidic, circumneutral or basic aqueous conditions. Likely, the most common of the geochemical process causing ARD is the weathering of potentially acid generating minerals (PAGs) and in particular sulfides. PAGs oxidation triggers several other processes, such as metal leaching or carbonate and silicate dissolution to compensate water acidity. The weathering rate of PAGs depends on multiple aspects such as textural properties of the minerals, the presence of water and oxidizing agents (mainly oxygen) existing or supplied to the reactive sites and the presence of acidophilic bacteria. For instance, pyrite oxidation can be described through the following reactions (Harries and Ritchie 1981): where DH is the enthalpy change.
Heterogeneity with respect to the spatial variability of primary mineral contents in the mine waste occurring during the mine waste construction or deposition (i.e., the pre-weathered mineralogical content) is a first element controlling the extension of ARD (Morin and Hutt 2000). Such variability can be a signature of the randomness of the nature, and it is usually associated with the formation of the soils and rocks owing to long-term geological processes. Mineralogical heterogeneity can occur at any scale of interest, from impurities (Savage et al. 2008) and overgrowth textures (Weisener and Weber 2010) on sulfide grains to the variability in average mineralogical content within kilograms of rocks, as evaluated during static or kinetic tests (Parbhakar-Fox and Lottermoser 2015). Heterogeneity in primary minerals is also observed at larger scales, from experimental to operational-scale (fullscale) mine waste sites (St-Arnault et al. 2020;Vriens et al. 2019). Since fine-grained materials are generally much more reactive than coarse grained materials (Strömberg and Banwart 1994), there is a strong connection between physical and geochemical heterogeneities that can have a key implication for the release and mobility of ARD in heterogeneous waste rocks. Heterogeneity associated with the distribution of secondary mineral contents (i.e., the post-weathered mineralogical content) is another important element controlling the extension of ARD (Servida et al. 2013;Wisotzky 1994). Al et al. (2000) concluded that ''reactions with the secondary minerals, rather than the primary mineral substrate, probably represent the principal controls on trace-element distributions in the pore water''. This conclusion seems supported by several other studies (e.g. Laurenzi et al. 2015; van der Sloot and van Zomeren 2012).
From a modeling perspective, embedding geochemical reactions in predictive heterogeneous waste rocks is computationally challenging. Early models embedded largely simplified geochemical systems compared to the multicomponent and strongly nonlinear reactions characterizing true geochemical systems. The TOUGH-AMD model (Lefebvre et al. 2001a, b) was developed and used to evaluate the role of fine-grained zones as capillary barriers impeding oxygen penetration across waste rocks dominated by diffusion and convective transport regimes. Their results indicated that the presence of such physical heterogeneities strongly limits the distribution of convective cells, reducing in turn the generation of ARD. However, the geochemical formulation in TOUGH-AMD is also limited to a few selected processes (Lefebvre et al. 2001a, b), particularly those generating sulfide leachate, without the possibility to consider, for instance the effects of mineral buffering or the formation of secondary minerals.
One of the first deterministic multicomponent reactive transport model on waste rocks was presented by Molson et al. (2005), who coupled HYDRUS-2D and POLYMIN model. Their study provided important evidences of the role of heterogeneity, included for instance the impact of preferential flow for the assessment of ARD. Molson et al. (2005) coupled unsaturated flow models with a reactive model that accounted for the most important elements controlling the formation of ARD from the weathering of sulfide-rich waste rocks, including oxygen supply and the presence of buffering minerals, such as calcite. Heterogeneity was imposed to these models by layering the textural properties of the waste rocks, and in addition by inclining these layers. The results highlight that heterogeneity has a key influence on the overall emergence of ARD. Molson et al. (2005) concluded that ''it is necessary to understand the moisture distribution and flow system together with the grain size distribution, oxygen concentrations, geochemistry and waste rock mineralogy''. For example, low-pH water which accumulates within stagnant flow zones will not contribute to ARD loading.

Quantifying the uncertainty: stochastic modeling
When dealing with geological media, an individual model outcome cannot be considered as the valid result for decision making purposes. Indeed, the lack of exhaustive knowledge of the subsurface often renders the mapping of properties and/or the parameterization adopted to describe a variable uncertain and prone to error, leading to both epistemic and aleatory uncertainties (Christakos 2000;Tartakovsky 2013). For predictive purposes, stochastic modeling provides a more versatile manner to deal with such uncertainty. In stochastic modeling, a probabilistic range of possible model results is provided instead of a single model output resulting from deterministic modeling. This enables the modeler to estimate not only the expected value of the distribution but also, and more importantly, the range of uncertainty surrounding the expected value. Among the multiple approaches for stochastic modeling in subsurface hydrology, the most popular is likely Monte-Carlo (MC) based modeling embedding geostatistical descriptions of one or more random spatial properties of the aquifers. In MC modeling, instead of recreating only one possible version of the reality, a set of equally plausible ''worlds'' is reproduced. Real-life sparsely distributed measurements and observations can be fully incorporated in (or ''honored by'') the models. Conditional or unconditional simulations are described through probability density functions and other relationships, such as spatial correlations, that allow for a statistical mapping of the model outputs. Chief among these properties is the variability of K used in stochastic modeling since the 1970s (Freeze 1975). However, any variable can be virtually chosen as random spatial functions. For instance, waste rocks are partially water-saturated deposits, hence, in the stochastic analyses of such unsaturated systems, the full range of parameters characterizing the characteristic unsaturated curves, such as the van Genuchten (1980) model parameters, can be treated as random spatial functions (e.g. Russo 2012).
In the context of ARD modeling in waste rocks, limited stochastic models have been presented in the literature. One reason is that the variability of the model parameters is not only limited to physical properties of the waste rocks (e.g. the hydraulic parameters) but also on geochemical parameters (e.g. Gerke et al. 1998). This issue poses additional challenges in the numerical solution of the fully coupled partial differential equations within MC context. However, stochastic reactive transport modeling remains particularly useful for predictive purposes under uncertainty. This type of modeling is appealing for risk-based analysis, given that a result can be expressed using a range of possible values, with a quantitative description of the likelihood of ARD (Pedretti et al. 2017a).
The following sections provide a summary of the stochastic analyses that have been presented tackling ARD from mining waste deposits. A summary of the main stochastic reactive transport models and related parameters can be found in Table 1. Note that not all the analyses presented below were developed with the explicit goal of quantifying the uncertainty in ARD. Yet, these analyses considered the randomness of model inputs and outputs, and therefore are more appropriate to uncertainty analysis than the modeling studies presented in Sect. 2 (even though in some cases the codes or algorithms were the same).

Approaches based on Eulerian modeling
Most stochastic analyses on heterogeneous multidimensional and multicomponent reactive transport models for the assessment of ARD from waste rocks have relied on Eulerian-based modeling. The main advantage of Eulerian models is that they can incorporate the full range of physical and geochemical processes occurring in waste rocks. The main drawback is that these models require high computational costs that may prevent their use for stochastic applications (particularly those involving many model solutions, such as in MC frameworks). Gerke et al. (1998) presented one of the first-documented multidimensional multicomponent stochastic-based reactive transport models for the assessment of ARD from waste rocks. Their focus was on physical and chemical heterogeneity caused by mixing of soil materials that may have already been oxidized to different degrees during the deposition of the spoil pile. To this end, they generated a two-dimensional layered flow and reactive transport model where the generation of ARD occurred under diffusive gas transport conditions. According to Gerke et al. (1998), this assumption ignores geochemical reaction kinetics and the role of bacteria by assuming the oxidation rate at the mineral grain surface is fast compared to the oxygen diffusion rate. In order to account for the depositionally induced structure of the spoil pile, they considered the relative (unsaturated) permeability K r , the unoxidized grain radius in the shrinking-core model (r c ), and the fraction of sulfide mineral (f s ) as spatially variable random fields. An unconditional geostatistical field (i.e. an ''individual realization'') was constructed using the Sequential Gaussian Simulations (SGSIM) module of GSLIB (Deutsch and Journel 1998), using an exponential autocorrelation approach. While SGSIM can create anisotropic correlated random spatial fields, Gerke et al. (1998) opted for spatial distribution of the fraction of sulfur to be more dependent on the properties of the geologic material than the deposition technique. Thus, they chose autocorrelation lengths for f s to be uniform throughout the section, while those for K r and r c to differ for the two layers in accordance with the different methods of overburden deposition.
The results by Gerke et al. (1998) suggested that the main effects of heterogeneity compared to a 1-D homogeneous model embedding equivalent properties appeared to be a strong ''smearing'' of the oxidation front and pronounced variations in concentrations throughout the respective buffering zones. The minimum pH of 2.0 arrived approximately 35 years later than in the homogeneous case, and the return to neutral conditions at the 50 m depth occurs 50 years later. Note that the simulation time of the model by Gerke et al. (1998) was reported as about 2.1 days for one individual realization, making it prohibitively at that time to perform a true stochastic analysis based for instance on multiple realizations, as in a MC framework. Fala et al. (2013) developed a stochastic study based on their 2D unsaturated models presented above Molson et al. 2005). They extended the homogeneous models to incorporate a geostatistical variability of K, suction (w), and volumetric water content (h). As in Gerke et al. (1998), Fala et al. (2013) used an anisotropic exponential autocorrelation function to generate individual sequential Gaussian simulations. Four different combinations of correlation lengths (k) of the three selected random properties were used to generate layered, vertically oriented and mixed distribution of parameters. Each field was then implemented in an individual reactive transport model and the results were compared. Table 1 Summary of stochastic reactive transport modeling analyses presented since 1997. SGS = Sequential Gaussian Simulations; SIS = Sequential Indicator Simulation; MPS = Multiple Point Geostatistics. n.c. = data not readily available from references. For the specific name of the variables, acronyms of models and specific information, we refer the reader to the text. Parameters: l = mean; r 0 = nugget, r = standard deviation, k = correlation length; v = proportion of matrix flow

References
Approach / adopted codes Geostatistical model ( The results by Fala et al. (2013) suggested that ARD is strongly sensitive to the selected material in the piles and the orientation and correlation lengths of the anisotropic autocorrelation model. For instance, in the layered case, the flow field is predominantly subvertical (due to the gravitydriven conditions imposed to the model), however k and the variance of the properties (r 2 ) generate local deviations which affect, for instance, the impact of preferential flow. The rates of sulfide mineral oxidation are influenced by the grain size, sulfide fraction and moisture content distribution. An important aspect discussed by Fala et al. (2013) is that the multiple realizations should be compared to find the general effects of heterogeneous properties on the expected behavior of the waste rocks. They finally indicated that the simulated models ''can take many days'', suggesting that ''it may not be practical to conduct fully stochastic analyses (for most situations)''. As discussed later, we contend that this may be the case today, however the increasing capacity of computational machines can generate a breakthrough also in the extended use of stochastic methods for ARD studies, bypassing the difficulties to computationally resolve MC simulations of nonisothermal multicomponent, multiphase, and multidimensional reactive transport of systems with complex geometries. Lahmira et al. (2017) extended the model setup presented by Lefebvre et al. (2001a, b) to study the effects of a random variability in physical properties within the waste rocks on the resulting ARD. They kept an identical conceptual model and base model setup, which consists of a fine-grained layer at the top of the pile with a lateral inclined batter. In Lahmira et al. (2017), the authors compared three random spatial distributions of the heterogeneous materials. The realizations represent a mixed distribution without a specific correlation pattern, as no geostatistical model has been used to create the random spatial fields.
The results by Lahmira et al. (2017) confirmed that lowpermeability compacted layers strongly limit convection, as in the deterministic models. The presence of fine-grained material near the boundary of a pile can limit air entry while beneath the pile surface favors the internal condensation of water vapor and thus minimize water loss. Coarse materials on the other hand seem to promote preferential flow of gas and water vapor. Lahmira et al. (2017) indicated that secondary gas convection cells can exist in the pile because of heterogeneity. The initiation of such secondary cells requires a minimal degree of heterogeneity, which is closely related to the ratio of permeabilities between coarse and fine materials. Lahmira et al. (2017) pointed out that the presence of coarse-grained material in the pile does not necessarily lead to more convection and higher ARD production, as the magnitude of convection seems to be controlled more by the amount of fine-grained material and its distribution in the pile. Specifically, the authors conclude that the presence of fine and very fine materials that remain at high water saturations constitute barriers to gas flow.

Approaches based on Lagrangian modeling
Lagrangian, or particle-based, models are alternative options to Eulerian based models which are particularly appealing for stochastic reactive transport modeling. Compared to Eulerian based mechanistic models, in Lagrangian modeling the complexity of a fully-developed heterogeneous flow and reactive transport model is usually simplified, resulting in efficient computational time. Lagrangian models do not yet include the full range of nonlinear processes controlling the formation and fate of ARD in waste rocks, although research is very active and promising regarding the possible use of Lagrangian models for multicomponent reactive transport (Engdahl et al. 2019;Henri and Fernàndez-Garcia 2015;Lu et al. 2018;Schmidt et al. 2019). Engdahl et al. (2017) provides an example of application of a particle-based model analysis including most of ARD generating processes in a sulfide-rich mining environment.
Lagrangian models embed statistical distributions of selected parameters that enable reproducing the salient random nature of a variable of interest. As such, they become very useful for uncertainty analyses that require multiple realizations to be resolved in computationallyefficient times. A classic example is for instance the amount of mass or the concentration of solutes passing a control section over time, i.e. the BTC. This BTC can be a random output of a model that embeds a random distribution of key parameters controlling the uncertainty (e.g. K), while leaving the other parameters constant Fiori et al. 2017;Pedretti and Bianchi 2018).
One of the first effective stochastic models for waste rocks was presented by , who developed a stochastic dynamic modeling framework focused on the implication of flow heterogeneity on copper leaching through the Aitik waste rock heaps. The approach by  was based on a probabilistic Lagrangian framework for reactive subsurface transport. Although the conceptual model is simplified compared to a real-case scenario, the Eriksson and Destouni (1997) model emphasizes the need to account for a distribution of possible model outcomes with a given range of uncertainty rather than a deterministic result. Their primary variable of interest was the water travel time, t, which was expressed by a bimodal probability density function (pdf), f t ð Þ, of form: where v is the fraction of slow flow path in the system, i.e. the matrix flow (1 À v being the proportion of preferential flow paths), while t G and r 2 t are the mean and variance of the population of residences time, respectively for slow flow paths (i = 1) and the population of preferential flow paths (i = 2). The ensemble of values t is then embedded into a 1-D analytical solution of the reactive transport, in which, t becomes a scaling variable controlling the firstorder kinetics of copper dissolution rates. The model also considers the equilibrium precipitation of secondary copper minerals.
The results by  showed that a stochastic selection of input parameters controlling preferential flow generates substantially different results than a deterministic selection of these parameters. The model by  has inspired a number of subsequent model-based analyses of ARD in mining sites (e.g. Blackmore et al. 2018;Malmström et al. 2008;Pedretti et al. 2017aPedretti et al. , 2020 and is considered a pioneering work in this field. The approach by  has the main advantage of being computationally efficient, given that most of the key equations were linearized by the physically-based approximations and assumptions made by the authors. However, such simplification reduced the number of processes compared to a real site. Indeed, the model does not ''predict'' long-term copper leaching from the Aitik site. As acknowledged by the authors, this would have required a detailed geochemical model that had accounted for multiple processes, which cannot be treated analytically. Rather, the model showed that probabilistic Lagrangian models can serve to obtain a sensitivity to selected model variables, which help focusing field studies in waste rocks, such as tracer tests in the case of the assessment of preferential flow.

''Hybrid'' approaches
Lagrangian approaches account for oversimplified geochemical systems to be predictive for waste rock systems, while Eulerian models are still too computationally demanding to perform a true stochastic analysis based on MC simulations. Although the capacity of computational machines is always evolving, as discussed later in the paper, nowadays alternative approaches must be sought if MC models are performed for routine applications, such as quick uncertainty analysis performed by practitioners, which may not have access to special computation platforms (e.g. supercomputers) or have more limited programming and modeling skills than advanced researchers.
In a series of recent papers, Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020 presented and applied an efficient stochastic modeling framework able to account for the majority of mechanisms leading to ARD from waste rocks within a probabilistic context. The method is promising to obtain a quick estimation of the likelihood of ARD occurrence in mineralogically and physically heterogeneous context. The method includes ARD-and waste-rockspecific processes, such as rate-limited gas transport or multicomponent geochemical systems. This approach portrays a multidimensional system as discretized into a set of flow-paths, loosely described as ''streamtubes''. After splitting the system into N ''streamtubes'', stochastic modeling was conducted considering each streamtube as a 1-D reactive transport problem, along which ARD forms from the interaction between gas, infiltrating and pre-existing pore water and primary and secondary waste rock minerals. The amount of streamtubes depends on the degree of flow heterogeneity, which is described via travel time pdfs, i.e. in the spirit of the Lagrangian approach by . However, the solution of transport in each streamtube is achieved via Eulerian modeling. As such, the stochastic modeling framework by Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020 can be seen as an ''hybrid'' modeling approach. In Pedretti et al. (2020Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020) the multicomponent multiphase reactive transport code MIN3P (Mayer et al. 2002) was utilized with simple MATLAB-based scripting for input and output processing. Alternative solutions for the implementation of the geochemical solvers in the multipurpose reactive transport simulations and parallelization can be used for similar purposes. For instance, when adopting codes from the family of geochemical code PHREEQC (Parkhurst and Appelo 2013), two new modules, IPhreeqc (Charlton and Parkhurst 2011) and PhreeqcRM (Parkhurst and Wissmeier 2015) are specifically designed to access PHREEQC's reaction capabilities from any scripting languages. This facilitates a faster communication between the geochemical code and transport simulators without the need of reading/writing external files (Muniruzzaman et al. 2020;Rolle 2016, 2019;Rolle et al. 2018;Sprocati et al. 2019).
The results from the framework proposed by Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020 are presented in the form of pdfs or cdfs of the ensemble of model outputs of target variables, as in any other stochastic analysis. These model outputs are derived from the mixed concentration of components forming the geochemical models and arriving at a specific control section (e.g. the bottom of the waste rock piles), where they mix over time using a flux-weighted formulation. So far, the approach has been successfully applied to estimate the effective neutralizing capacity of heterogeneous waste rock piles. Results are expressed as pdfs of the concentration of a solute component (e.g. SO 4 ) or the pH. For instance, the expected pH after 100 years from the starting of waste rocks weathering is obtained, along with the degree of uncertainty (e.g. the 95% confidence interval).
The key advantage of the approach by Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020 is that it is extremely efficient from a computational perspective, being easily parallelizable. A MC analyses composed by 100 realizations, each one composed by 100 streamtubes, correspond to 10,000 individual 1-D multicomponent reactive. For instance, tested simulations in Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020 ran overnight to complete a MC analysis on a workstation based on 4-cores IntelÒ XeonÒ E-2124G 3.40 GHz. The main drawback is that the discretization into 1-D models is a necessary simplification to avoid the problems linked to multidimensional Eulerian models. The limitations of this simplification are described in the cited literature.

A perspective on the evolution of the discipline
For several decades, research and applications in mining and hydrogeology have been evolved separately. Nowadays, thanks to multiple projects sponsored by mining companies looking at integrating research and mining practices, international conferences and workshops, scholars and practitioners seem to converge on the idea that a better understanding and more reliable prediction of acid rock drainage (ARD) requires multidisciplinary approaches able to bridge cutting-edge scientific advances with practical needs. Mathematical modeling has been part of this evolution and mechanistic models can now embed almost all known partial differential equations controlling the main physical and geochemical processes generating ARD in waste rocks (Amos et al. 2015). The increasing capability of workstations and supercomputers is however the key aspect to define the perspectives regarding the evolution of the discipline in the incoming future. The following is the authors' perspective regarding possible topics that will be progressively addressed by researchers and practitioners in the following years.

Communicate the need of stochastic modeling beyond the academic boundaries
When dealing with waste rocks modeling, it must be understood by practitioners that the result from any model will be uncertain. If deterministic analysis is used, that results will be certainly wrong and possibly too unreliable to make a decision out of it. The first reason is that waste rocks are complex environments in which physical, chemical, and biological processes are interrelated (Fig. 1) and the spatial and temporal variability of the properties controlling these processes are often insufficiently characterized to the scales which are relevant to make use of deterministic modeling approaches. The second reason is that, even when the full range of virtually-available spatially-distributed heterogeneous properties of a waste rock were known, the modeler could not incorporate this knowledge into a code, as it would result in too high computational burden. Stochastic modeling has been proved to assist modeler when dealing with these difficulties, such as through moment-based analysis to cope with data limitations or via parameter upscaling. Stochastic models generate probabilistic results, such that decision makers can decide, for instance, whether a waste rock deposit will generate some polluted drainage, how much it will cost to remediate to such contamination or how to ensure that any drainage is properly controlled. All these decisions will be made with a known degree of uncertainty.
Unfortunately, stochastic hydrogeological modeling struggles to be fully accepted beyond the academic boundaries (Bode et al. 2018;Freeze 2004;Renard 2007;Sanchez-Vila and Fernàndez-Garcia 2016), and convincing mining practitioners is not an easy task. In the words of Sanchez-Vila and Fernàndez-Garcia (2016), who rephrased Renard (2007), modelers must understand that ''deterministic models do not represent reality at all. The reason is the combination of unsampled natural heterogeneity and scenario uncertainty. Only stochastic models have a chance of providing the answers needed for proper groundwater management efforts.'' While recent works (Fala et al. 2013;Lahmira et al. 2017;Pedretti et al. 2017a) suggest that it is already possible to use stochastic reactive transport modeling for predictive purposes in waste rocks, the true ''revolution'' is expected to occur when modelers in consultant agencies will be convinced by researchers and, more importantly, will have easy access to parallel computing for reactive transport modeling. Parallel computing has bloomed years ago (Steefel et al. 2014), and several codes have been adapted since then (e.g. Parkhurst and Wissmeier 2015;Su et al. 2017). Monte Carlo modeling on multidimensional multicomponent reactive transport models for waste rocks will be possible to generate stochastic scenarios of multiple realizations, which will provide expected behavior of the waste rocks along with quantified degree of uncertainty. For instance, this would help by-passing the simplifications of the approach proposed by Pedretti et al. (2016Pedretti et al. ( , 2017aPedretti et al. ( , b, 2020.

Pore-scale analyses
ARD processes are controlled by microscopic phenomena (Pantelis and Ritchie 1992), while models are usually performed at much larger scales that lumps together multiple reactive sites. This is one of the reasons leading to the well-known scale dependence in reaction rates and a major debated issue in waste rocks studies (Amos et al. 2015).
We contend that in the next years the emerging characterization of modeling of flow and reactive transport heterogeneities at the pore scale (Blunt et al. 2013) will be a key focus to reduce uncertainty in waste rock modeling (and ARD modeling in general). One reason is the development and diffusion of modern machines for pico-and nano-scale imaging methods, in particular X-ray micro tomography (Sayab et al. 2015) and similar methods. These techniques have now allowed acquisition of three-dimensional reconstructions from a series of two-dimensional projections taken at different angles. Another major breakthrough in this discipline has been the increased availability of high resolution experimental investigations at the pore-scale, such as microfluidic experiments (de Anna et al. 2014;Zhang et al. 2010), or at the Darcy scale including spatially resolved measurements (Muniruzzaman et al. 2014;Rolle 2015, 2017;Rolle et al. 2013).
Enhanced capability to resolve such small-scale physicochemical processes will lead to significant development not only in establishing process understanding but also in effective upscaling the effects of ARD mechanisms during large-scale transport. At the same time, the increasing capability of computational methods to resolve partial differential equations, such as the Stokes equation, on high-resolution grids has allowed a prompt blooming of studies focusing, for instance, on the role of local dispersion, reactive transport, or mineral-solution interactions at the scale of the pore (Yoon et al. 2015), thus circumventing the problems of scales in hydrogeological modeling.

Connectivity and spatial order
Transport in heterogeneous media is very sensitive to the degree of connectivity of the clusters (Bianchi et al. 2011;Pedretti 2017, 2018;Knudby et al. 2006;Molinari et al. 2015;Renard and Allard 2013) and to the spatial ordering of the structures forming a geological medium. Connectivity, clustering and spatial ordering can be measured in multiple manners, e.g. through the analysis of geological entropy Pedretti 2018, 2017;Pedretti 2020;Pedretti and Bianchi 2019;Ye et al. 2018). However, connectivity is not yet an established aspect for waste rock modeling. It is expected that this will become an important issue for future development in the area, particularly related to the impact of preferential flow and channeled transport in the waste rocks. One of the most recent works in this sense is Appels et al. (2018). Although they did not implement a reactive transport model analysis as in the above-mentioned contributions, their work was one of the first to adopt multiple point geostatistics (MPS), which allow to incorporate connectivity of textural properties and associated hydraulic and transport properties (e.g. K clusters). In particular, the connectivity of extreme patterns is honored by MPS compared to, for instance, Sequential Gaussian Simulations, where extremes tend to be smoothed out and underestimate the impact of connected features on transport (e.g. Renard and Allard 2013).

Biologically-controlled kinetics
Microbially-controlled reactions have a dominant role to control the scales of ARD in waste rocks. The spatial variability of biologically-related parameters can be as important as the variability of physical and geochemical parameters. Yet, little evidences of waste rock models embedding heterogeneous distribution of parameters directly related to microbial activity and their impact on ARD were found while reviewing the literature for this paper. The most adopted approach so far seems to be the use of effective kinetic rates that lump together both the reaction rates associated to geochemical properties of waste rocks and the effects of microbially-induced reactions.
Models nonetheless exist to include explicitly the spatial effects of microbially-induced reactions, for instance as shown by Casas et al. (1998). Modeling the spatio-temporal variability of microbial colonies using stochastic approaches could enable reproducing the uncertainty associated to the scarce knowledge of the biologically controlled kinetics of ARD generation. A key starting point could be the definition of model zonations based on geostatistical distribution of ''acidic microenvironments'' (Dockrey et al. 2014) at sulfide mineral surfaces in order to catalyze the sulfide oxidation processes even under circumneutral pH conditions (Southam and Beveridge 1993). Such acidic microenvironments are thought to be sustained by the simultaneous occurrence of ferric (oxy)hydroxide and/or (oxy)hydroxysulfate precipitation, which allows an isolation of these microenvironments from the bulk pH-neutral pore water solution within the close proximity of the sulfide surface Pace et al. 2018).
A second approach could be based on the simulation of the spatial distribution of the habitats for these colonies. For instance, the porous layer of the secondary Fe 3? phases (e.g., K-jarosite, schwartmannite) essentially provides important habitats for developing microbial colonies and biofilms of acidophilic bacteria under bulk non-acidic conditions (Dockrey et al. 2014). Consequently, in these so-called acidic microenvironments, Fe 3? solubility is increased, sulfide weathering is locally enhanced, and the pH limiting growth condition of acidophilic species is overcome (Dockrey et al. 2014). Additionally, these bacterial species can deposit extracellular polymeric substances (EPS), which are also known to complex with Fe 3? and subsequently enhance sulfide dissolution kinetics by making it available for sulfide oxidation reactions (Sand et al. 1995). the structure and content of this paper. MM has been supported by internal funding for research and development (TUMMELI project) at the Geological Survey of Finland, European Regional Development Fund (KaiHaMe and Kaivasu projects). The results here presented have been developed in the frame of the MIUR Project ''Dipartimenti di Eccellenza 2017-Le Geoscienze per la società: risorse e loro evoluzione''.
Funding Open access funding provided by Università degli Studi di Milano within the CRUI-CARE Agreement.

Compliance with ethical standards
Conflicts of interest The authors declare that there is no conflict of interest regarding the publication of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.