Comparison of unsaturated flow and solute transport through waste rock at two experimental scales using temporal moments and numerical modeling

Journal of Contaminant Hydrology 171 (2014) 49–65

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology journal homepage: www.elsevier.com/locate/jconhyd

Comparison of unsaturated flow and solute transport through waste rock at two experimental scales using temporal moments and numerical modeling Sharon Blackmore ⁎, Leslie Smith, K. Ulrich Mayer, Roger D. Beckie The University of British Columbia, Vancouver, BC, Department of Earth, Ocean, and Atmospheric Sciences, 2020-2207 Main Mall, Vancouver, BC, Canada, V6T 1Z4

a r t i c l e

i n f o

Article history: Received 14 April 2014 Received in revised form 8 October 2014 Accepted 11 October 2014 Available online 23 October 2014 Keywords: Waste rock Scale Tracer Dual-porosity Mobile–immobile (MIM) Temporal moments Mass-transfer

a b s t r a c t This study analyzed and compared unsaturated flow response and tracer breakthrough curves from a 10-m high constructed pile experiment (CPE) in the field (Antamina, Peru) and two 0.8 m high laboratory-based columns. Similar materials were used at both experimental scales, with the exception of a narrower grain size distribution range for the smaller column tests. Observed results indicate that flow and solute transport regimes between experimental scales were comparable and dominated by flow and solute migration through granular matrix materials. These results are supported by analogous breakthrough curves (normalized to cross-sectional area and flow path length) that suggest observation- or smaller-scale heterogeneities within the porous media have been homogenized or smoothed at the transport-scale, long breakthrough tails, and similar recovered tracer mass fractions (i.e., 0.72–0.80) at the end of the experiment. CPE breakthrough curves do indicate a portion of the fluid flow follows rapid flow paths (open void or film flow); however, this portion accounts for a minor (i.e., ~0.1%) component of the overall flow and transport regime. Flow-corrected temporal moment analysis was used to estimate flow and transport parameter values; however large temporal variations in flow indicate that this method is better suited for conceptualizing transport regimes. In addition, a dual-porosity mobile–immobile (MIM), rate-limited mass-transfer approach was able to simulate tracer breakthrough and the dominant transport regimes from both scales. Dispersivity values used in model simulations reflect a scaledependency, whereby column values were approximately 2× smaller than those values applied in CPE simulations. The mass-transfer coefficient, for solute transport between mobile and immobile regions, was considered as a model calibration factor. Column experiments are characterized by a larger “mobile to immobile” porosity ratio and a shorter experimental duration and flow path, which supports larger mass-transfer coefficient values (relative to the CPE). These results demonstrate that laboratory-based experiments may be able to mimic flow regimes observed in the field; however, the requirement of scale-dependent dispersivities and mass-transfer coefficients indicates that these tests may be more limited in understanding larger-scale solute transport between regions of different mobility. Nevertheless, the results of this study suggest that the reasonably simplistic modeling approaches utilized in this study may be applied at other field sites to estimate parameters and conceptualize dominant transport processes through highly heterogeneous, unsaturated material. © 2014 Elsevier B.V. All rights reserved.

⁎ Corresponding author at: BGC Engineering Inc., 800-1045 Howe St., Vancouver, BC, Canada, V6Z 2A9. Tel.: +1 604 684 5900x41275. E-mail addresses: [email protected] (S. Blackmore), [email protected] (L. Smith), [email protected] (K. Ulrich Mayer), [email protected] (R.D. Beckie).

http://dx.doi.org/10.1016/j.jconhyd.2014.10.009 0169-7722/© 2014 Elsevier B.V. All rights reserved.

1. Introduction The construction of waste rock piles is part of the mining process at many mine sites. These piles are generally large, unsaturated structures containing material with grain diameters

50

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

spanning at least six orders of magnitude (i.e., 10−6 m–1 m; Nichol et al., 2005; Fala et al., 2005). Predictions of solute loads released at the base of a waste rock pile require the characterization of fluid flow and solute transport. Estimating solute transport parameters for use in predictive models becomes increasingly difficult at larger scales due to the effects of spatial variability in the many factors that control the infiltration of water and the mobility of solutes (e.g., initial/boundary conditions, material properties, preferential flow paths; see Jury and Flühler, 1992; Li and Ghodrati, 1997; Vanderborght et al., 2000; Vanderborght and Vereecken, 2007). These complexities observed at the large scale are generally viewed as difficult, or impossible, to reproduce in smaller experiments (e.g., field barrels, laboratory columns, or humidity cells). For example, a number of studies have compared mineral weathering rates across scales with the intention of determining correction (or scaling) factors (Eriksson et al., 1997; Frostad et al., 2005; Malmström et al., 2000; Otwinowski, 1995). However, the results reveal large discrepancies in behavior among scales, which suggest large field experiments may be necessary to adequately characterize system response (Malmström et al., 2000). The material characteristics of waste rock can have a substantial impact on its hydrological response (Smith and Beckie, 2003). This observation is supported by the study of Butcher et al. (1995), who state that in stony soil the flow path of water is an intrinsic property of the soil medium (at a given water content). Flow through waste rock can be categorized into two types: preferential flow through either open voids and/or the coarse-sized fraction of the matrix, and matrix flow through finer-grained (or matrix) materials. Preferential flow paths can be expressed as channelized flow whereby the water flux is concentrated into spatially distinct areas smaller than the total crosssectional flow area (Nichol et al., 2005). In this paper, the term matrix flow describes the movement of water through the granular material present between the cobbles and boulders and is dictated by moisture contents, capillary tension and gravity. The occurrence of some degree of preferential flow in unsaturated porous media has long thought to be the rule rather than the exception (e.g., Brusseau and Rao, 1990; Flury et al., 1994). In general, three factors exert a strong control on the degree of preferential flow in waste rock: the infiltration rate, the spatial arrangement of matrixsupported and matrix-free zones, and the particle size distribution (PSD). The first factor is dependent upon the local climate and is not pertinent to this study as infiltration conditions are similar for all experiments reported here. For the second factor, Herasymuik (1996) observed that the method of construction controls the spatial orientation of material within the pile. For example, waste rock piles that are created by enddumping may result in a fining-upwards gradation, as coarser material falls to the base of the lift and finer materials generally remain near the surface (Fala et al., 2005). An alternative is pushdumping, whereby waste rock is deposited in a series of lifts starting at the top of a pile and “pushing” the material to a level surface. This method generally results in a coarse lower zone and a non-uniform upper zone with horizontal traffic surfaces between lifts (Corazao Gallegos, 2007). The PSD of waste rock is a function of its mineralogical composition, material friability, mine-specific blasting

techniques, among other factors. PSD analyses generally encompass size fractions that range over several orders of magnitude and, regardless of the coarse proportions, typically exhibit a long tail at the finer size fractions (Smith and Beckie, 2003). Structured porous media exhibiting flow heterogeneities, or a combination of preferential and matrix flow paths, are frequently described by multi-domain models (Gerke and van Genuchten, 1993a). In the simplest multi-domain model (i.e., dual-porosity), water/ solute are partitioned between a mobile and an immobile domain. Using this approach, the mobile domain constitutes preferential flow and matrix flow paths, and the immobile domain includes stagnant (or non-flowing) waters. The immobile domain is a dynamic region whereby water and/or solute may transfer in and out of this domain, as a result of pressure or concentration gradients (Šimůnek and van Genuchten, 2008). Model formulations using a fixed water content in the immobile domain and solute only transfer between the immobile and mobile domains represent the most basic dual-porosity approach and are termed mobile–immobile (MIM) models (e.g., van Genuchten and Wierenga, 1976). A dual-permeability model also assumes two overlapping pore domains, like the dual-porosity approach, but replaces the immobile domain with a low-permeability domain, which allows active, albeit slow water movement (Gerke and van Genuchten, 1993a, 1993b; Šimůnek and van Genuchten, 2008). Water and solute transport in a dual-permeability model can occur within and between mobile and immobile domains. In this study, flow and tracer response are examined at two experimental scales, in 0.8 m tall laboratory columns and in a 10 m high experimental waste rock pile constructed using the end-dumping method. This study is akin to previous work in that it compares flow and conservative solute transport at two scales (Eriksson et al., 1997; Strömberg and Banwart, 1994, 1999a, 1999b), but differs in its effort to maintain similar compositions and proportions of the different waste rock types at both scales. The objective of this research is two-fold: 1) to determine the main flow regimes controlling solute transport in these waste rock assemblages under site-specific and laboratory conditions, and 2) to determine if the principal influences on solute transport linked to waste rock flow regimes in larger test piles can be adequately reproduced at more practicable scales by managing “controllable” variables (i.e., material compositions and precipitation inputs). Both experiments share the same rock type and the rainfall applied to the laboratory columns mimics the daily rainfall recorded at the experimental pile. Flow is recorded from basal lysimeters and a tracer test is used to characterize the partitioning of infiltration among different flow paths. Tracer results are analyzed using a flowcorrected time approach (Eriksson et al., 1997) and numerical modeling. The flow-corrected time approach has been employed in other studies with transient flow behavior (e.g., Destouni, 1991; Jury, 1982; Small and Mular, 1987; Wierenga, 1977). Numerical modeling aimed to apply the simplest multi-domain modeling approach (i.e., the MIM model) to permit the diagnosis of dominant features of solute transport processes, while minimizing the number of model parameters.

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

2. Materials and methods 2.1. Site description The experimental field site is at the Antamina Mine, located in north-central Peru (9° 32 S and 77° 03 W), approximately 270 km north of Lima, Peru. The average annual temperature and precipitation range between 5.5 and 6.0 °C, and 1200 and 1300 mm, respectively. Climate at Antamina is governed by the tropical rain belt, which creates distinct “wet” and “dry” seasons. Wet seasons occur between October and April, with the remainder of the year described as the dry season. Approximately 80–90% of the total annual precipitation at Antamina occurs during wet season months. Antamina divides its waste rock into three main classes based on geochemical parameters (Class A, B and C; Table 1). Class A and C material represent end members, with the former considered reactive and the latter less-reactive; Class B material is intermediate between the two. The majority of Antamina waste rock is carbonate-rich and contains variable sulfide mineral contents (i.e., N 3% Class A; b 3% Class B and C) and minor ore minerals (e.g., chalcopyrite, molybdenite, sphalerite) (Antamina, 2007). 2.2. Constructed pile experiments (CPE) From 2005 to 2009, five waste rock piles (or CPEs) were built at Antamina for the purpose of understanding controls on waste rock leachate quality and quantity (i.e., Pile 1 through 5). Each CPE is 36 m (l) × 36 m (w) × 10 m (h) and contains approximately 25,000 tonnes of material (Fig. 1-A). Piles were constructed through multiple tips (or tipping phases) of enddumped material producing approximately 37° slopes. The focus of this study is on Pile 5 (herein referred to as 'CPE'), which incorporates Class A and C waste rock. Details of the construction of all five CPEs can be found in Corazao Gallegos (2007) and Bay et al. (2009). This study's CPE (Pile 5) consists of six unique tipping phases, each containing 2600–4500 tonnes of Class A intrusive or Class C hornfels/marble material waste rock. The ratio of Class A to C material in this CPE is ~9:11 (as wt.%:wt.%) (Fig. 1A). Particle size distribution (ASTM D 5519–94) and bulk density (ASTM D 5030–89) were analyzed at the mine site (Golder, 2010). Tests were completed on samples consisting of 4–8 tonnes of waste rock. Averaged PSD curves of samples relevant to this study are shown in Fig. 2. Generally, a typical PSD curve of waste rock shows a well-graded distribution with a coefficient of uniformity Cu (=D60/D10) of 20 or higher (McKeowen et al., 2000; Morin et al., 1991). Calculated Cu values present very different physical properties for Class A and C materials. Class A material shows high Cu values (i.e., 257) and a well-graded distribution whereas Class C materials contain low Cu values (i.e., 12) and a poorly graded distribution.

51

The footprint of the CPE is underlain by a 75-mil geomembrane liner serving as a basal lysimeter. Within this large basal lysimeter are 3 smaller (4 m × 4 m) geomembrane liners located along the centerline of the CPE. These smaller catchments act as sub-lysimeters to characterize flow or infiltration through a smaller cross-sectional area within the CPE (Fig. 1-A). These smaller sub-lysimeters and the large basal lysimeter are referred to as Lysimeter A–C and Lysimeter D, respectively (Fig. 1-A) Drainage is channeled from the lysimeters to an instrumentation hut that records flow rate/volume, electrical conductivity, and temperature. The experiment was initiated in the dry season and corresponds to the first day of recorded data (i.e., June 1, 2009). The instrumentation hut includes sampling ports for analysis of drainage water quality. Each CPE contains six instrumentation lines with multiple sensors to record in situ temperature, electrical conductivity, gas concentrations and water content. Sensors recording in situ moisture contents (θ) [L3∙L−3] are used to determine initial conditions of waste rock material. Data from the remaining instruments (i.e., conductivity, temperature and gas concentrations) are not discussed in this study, but further details may be found in Bay et al. (2009) and Peterson et al. (2012). This study focuses on flow and transport to Lysimeter B of the constructed pile. Lysimeter B is located along the centerline of the pile and encompasses material from the first 3 enddumping events used to construct the CPE. In total, approximately 430 tonnes of waste rock overlies the 4 m (w) × 4 m (l) footprint of Lysimeter B in a 1:4 ratio of Class C material underlying Class A materials. Class C material was used in the first end-dumping event, whereas the second and third events contained Class A material. The placement of the two Class A materials was separated temporally by the installation of instrumentation line 1 (Fig. 1-A). Both Class A materials were excavated from the same location within the Antamina open pit but at 9 days apart. As such, these materials are assumed to be similar, and are not differentiated in this study. A lithium bromide tracer with 3050 mg L−1 of Br− was applied on the CPE on January 24th, 2010 or 237 days following experiment initiation. Bromide was selected as the tracer due to its low background concentrations at Antamina (i.e., 0.2– 0.8 mg L−1) and its conservative nature. Bromide was applied to the crown of the pile through a single applied rainfall event of 24 mm and 5 hours duration (4.8 mm∙hour−1). The application rate corresponds to a heavy precipitation event with a 4-year return period at the Antamina mine. The Christiansen uniformity coefficient value (CUC; Christiansen, 1942) of the tracer application was determined to be 77, which confirms a moderate to even tracer application across the pile crown. Observed ponding was minimal and no surficial runoff from the pile occurred during tracer application. Bromide concentrations from water samples were measured using a spectrophotometric method (modified from Presley

Table 1 Antamina waste rock classification (modified from Antamina, 2007). Class

Description

Rock types

Zn (ppm)

As (ppm)

Sulfides (%)

Visual oxides (%)

A B C

Potentially acid-generating Potentially acid-generating Non acid-generating

Skarn, intrusive, marble, limestone Marble, limestone, hornfels Marble, limestone, hornfels

N1500 700–1500 b700

N400 N/A b400

N3.0 b3.0 b3.0

N10.0 b10.0 b10.0

52

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

A A

C

C

Construcon of CPE with a enddumped p phase

A

10 m

Instrumentaon lines and drainage

Line 5

Line 6 Sub-lys B

Lysimeter D

Sub-lys A

Sub-lys C

36 m

0.45 m

B

Eaves

A-1

0.8m

A-1

A-2

A-2

C Silica sand Tipping Bucket

C Silica sand

Fig. 1. Large-scale constructed pile experiment (CPE) (A) and small-scale laboratory column (B), including complementary photographs of experiments for the purposes of scale reference. Note, box schematic in (A) represents design template for smaller scale experiment shown in (B).

(1971)). One of every 25 sample analyses was repeated for quality assurance. For a brief time period when tracer samples were inadvertently not collected, bromide concentrations were estimated by correlation with dissolved lithium concentrations, which is a weakly sorbing cation that was co-applied with bromide. Lithium concentrations were obtained from routine weekly water-quality samples collected and analyzed by the mine. This study's approach assumes vertical flow from the CPE surface to the free-draining sub-lysimeter B. Despite significant material property differences between Class A and Class C material, Blackmore et al. (2012) showed this CPE (i.e., Pile 5) is dominated by vertical flow regimes with a negligible horizontal flow component. Specifically, measured tracer concentrations from sub-lysimeter C, outside of the tracer application area, were at or below detection for the majority of the experiment and released an insignificant (i.e., 0.03%) amount of applied tracer mass over the duration of the experiment. 2.3. Laboratory columns The smaller scale column experiments were designed to closely resemble the material assemblage and structure of the CPE. Class A and C materials for the column study were transported to the University of British Columbia (UBC, Vancouver, BC, Canada). Class A material was comprised of two subsets, A-1 and A-2. Class A-1 is comprised of Class A waste rock from three different regions of Antamina's open pit, which includes material used in Pile 5. Class A-2 is composed of material from a single location and is more physically and geochemically

homogeneous than Class A-1. Although these materials are not an exact replica of waste rock overlying Lysimeter B in Pile 5, they exhibit the range of physical heterogeneity within the Class A waste rock classification, analogous to the materials used in the CPE. The Class C material used in the columns was collected from 2 tipping phases during construction of the CPE and is representative of the majority of the Class C material found in the Pile 5 CPE. Waste rock encompassing each of the three subsets (Class A-1, Class A-2 and Class C) was combined and cone-andquartered (Kellagher and Flanagan, 1956) prior to placement in two 0.8 m (h) × 0.45 m (d) columns (Fig. 1-B). One part Class C underlies two parts Class A, which is equally divided between the two Class A samples, similar in design to the layering in the CPE (i.e., Pile 5). A 5 cm layer of silica sand was placed at the base of the waste rock mixture to promote drainage to the basal outlet. The two columns are referred to as C1 and C2 and each contain ~170 kg of waste rock. Waste rock was placed in each column and the construction and placement of material required caution and care to reproduce the gross internal structure in terms of layering Class A over Class C materials, ensure a random representation of grain sizes from each material type, avoid disruption of surrounding materials and prevent damage to instrument installations. Eight sensors (Decagon ECH2O probes) are located within each column to record moisture content (θ) [L3∙L−3]. These sensors require an adequate hydraulic contact with the surrounding soil and are sensitive to air gaps and soil disturbance. Accordingly, the moisture content sensors were placed in a 1:1 ratio of silica

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

75

Cobbles

100

Gravel

4.75

Sand

Percent passing (%)

80

53

0.075 Fines (silt/clay)

CPE: Class A-avg Column: Class A-1 Column: Class A-2

60

40

20

A) Class A

0 1000

100

10

1

0.1

0.01

Parcle size diameter (mm) 75

Cobbles

100

Gravel

4.75

Sand

0.075 Fines (silt/clay)

80 Percent passing (%)

CPE: Class C-avg Column: Class C

60

40

20

B) Class C

0 1000

100

10 1 Parcle size diameter (mm)

0.1

0.01

Fig. 2. Average particle size distribution curves for Class A (A) and C (B) waste rock material used in CPE and column experiments. Note: n = 5 (CPE: Class A-avg); n = 4 (CPE: Class C-avg); n = 4 (Column: Class A-1); n = 1 (Column: Class A-2); n = 4 (Column: Class-C). Dashed lines represent min-max PSDs for averaged curves. The Unified Soil Classification System is shown on upper y-axes (USAEWES, 1960).

sand to sieved waste rock to ensure infiltrating waters contacted the sensors. A grain size limitation is imparted by the column dimensions and material is restricted to cobbles with diameters less than 10 cm. PSD curves for Class A and C material used in column experiments are truncated and normalized (i.e., 10 cm = 100% passing) in order to reflect grain sizes used at this scale and enable a better comparison with field scale PSD results (Fig. 2). Class A and C materials showed comparable Cu values to those calculated from CPE curves (i.e., Class A: 210 and 133 (columns — Class A1 and A2, respectively) versus 257 (CPE); Class C: 11 (columns) versus 12 (CPE)). These results indicate that the physical properties of Class A and C waste rock are comparable at both experimental scales, despite a grain diameter constraint imposed for the columns. Porosity and residual water content measurements were made in the laboratory by completely saturating oven-dried material to obtain water content at saturation θsat [L3·L−3],

which was assumed to equal porosity, and allowing this water to drain for residual water content θr [L3·L−3] (Table 2). Distinct PSD curves were generated for Class A-1 (n = 4) and A-2 (n = 1) materials and averaged PSD curves of Class C (n = 4) are applicable to both the CPE and column experiment (Fig. 2). The percentage of material less than 2 mm in Class A and C material is 23% and 4%, respectively, implying that Class A waste rock is considerably finer-grained. Each column contains a basal lysimeter, which channels all drainage to a tipping bucket for recording flow rate/volume and allows for collection of aqueous samples (Fig. 1-B). Precipitation events were applied using an air-atomizing spray nozzle (Lechler Airmist® Series 136.2), combining compressed air (at 100 psi) and flow rates of 100 mL min−1 (377 mm∙hour− 1). Rainfall was applied once daily and five times weekly to columns for durations of 40 seconds to 35 minutes. This range mimics dry season and wet season precipitation amounts (respectively) observed at Antamina

54

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

Table 2 Initial model simulation input parameters. CPE

Laboratory Column

Parameter

Class A

Class C

Class A-1

Class A-2

Class C

Equivalent soil type

Clay loam – loam

Sandy clay loam

Clay loam – loam

Clay loam

Sandy clay loam

1.8 1.3 0.06 1950 0.225 0.15 0.24

5.9 1.5 8.0 1950 0.17 0.06 0.29

1.9 1.3 0.06 1950 0.165 0.16 0.25

3.6 1.4 0.21 1960 0.16 0.16 0.26

5.9 1.5 3.0 1950 0.10 0.06 0.29

[m−1] [−] [m·day−1] [kg·m−3] [m3·m−3] [m3·m−3] [m3·m−3]

α n Ksat ρb a θ θrb θsatb a b

Golder (2010). Obtained from laboratory analyses.

over the June 2009–May 2010 water year. Applied rainfall yielded a moderate Christiansen uniformity coefficient (CUC value = 67). To avoid infiltration of precipitation along the column sides, a small eave was placed along the inner diameter of both columns (see Fig. 1-B). Eave catchment water was funneled to a collection bottle for volumetric measurement and corrected applied precipitation volumes reflect approximately 70–80% of initial applied volumes (i.e., 250–300 mm∙hour−1). Similar to the CPE, a tracer was applied in a single rain event of 10 mm over 0.2 hours on both C1 and C2 at the height of the wet season (i.e., t = 241 days). Bromide, as LiBr, was also used in the column experiments (Co = 3100 mg L− 1 Br−) and leachate concentrations were analyzed using the same spectrophotometric method as described for the CPE tracer test.

The expression for the first temporal moment ðτ Þ [T], or the arrival time of the center of tracer solute mass at depth z [L], is: Z

τSðτ; zÞ∂τ : τðzÞ ¼ Z 0∞ Sðτ; zÞ∂τ

e where In this equation (S) [M·T−1] is S(τ,z) = cðτ; zÞQ c(τ, z) [M·L−3] is the measured solute concentration at time τ e [L3·T−1] = Vtot/ϒ is the and a vertical travel distance z and Q steady-state water flow that corresponds to flow-corrected time, τ. In this analysis, the first moment is evaluated at the base of the pile or column (z = 10 or 0.8, respectively). Assuming that the conservative tracer is unaffected by reactions and mass diffusion into immobile domains (Eriksson et al., 1997), the τ value can be used to calculate the apparent mobile water content θm⁎ [L3·L−3], 

τ≡ϒ

V ðt Þ V tot

ð1Þ

where ϒ [T] is the experimental tracer-test period, V(t) [L3] is the cumulative flow up to time t and Vtot [L3] is the total flow volume for the experiment duration. Breakthrough curves represented as functions of τ are therefore stretched out during periods of high flow and compressed during low-flow periods (relative to real-time curves).

ð2Þ

0

2.4. Parameter estimation Eriksson et al. (1997) investigated the existence of preferential flow in waste rock at the Aitik mine site, Sweden using a flow-corrected time (τ) approach. This method accounts for temporal variability in flow and allows for the comparison of breakthrough curves through their temporal moments, from multiple experiments at different scales. The following is a brief discussion of the flow-corrected methodology. Flow-corrected time (τ) reflects time as a relative proportion of the total flow exiting a boundary relative to the duration of the tracer experiment and related to the time lapsed following tracer application (Eq. (1)):

∞

θm ¼

e τQ Az

ð3Þ

where A [L2] is the discharge area of drainage and z [L] is the travel distance. The degree of preferential flow (v; [−]) can be assessed by comparing the apparent mobile water content to the actual mobile content values θactual [L3·L−3], v¼

θm : θactual

ð4Þ

The actual mobile water content is assumed to equal the measured water content θmeas (from sensors) and measured water contents may be comprised of both immobile and mobile water domains (i.e., θmeas = θm + θim). However, the Eriksson et al. (1997) approach assumes no transfer between the immobile and mobile domains or the absence of an immobile domain (i.e., θactual ≈ θmeas = θm + θim; if θim ≈ 0 or negligible transfer) and is conceptualized as a uniform flow model. Indeed, in the presence of immobile porosity, the assumption that measured water contents equal actual water contents (to be used in Eq. (4)) results in an underestimation of v values. Therefore, the ratio of “apparent to actual mobile water contents” (i.e., v) can be used as an indicator of preferentially flowing water and/or an

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

immobile domain. A low ratio suggests that a system may be better represented by a dual-domain approach. The spreading of the breakthrough curves around the mean arrival time ðτ Þ is quantified by the second temporal moment (στ) [T] (Eriksson et al., 1997), as: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uZ ∞ u ðτ−τ Þ2 Sðτ; zÞ∂τ Z ∞ : στ ¼ u t 0 Sðτ; zÞ∂τ

ð5Þ

0

Relative spreading around the mean arrival time can be expressed by the normalized second moment (CVτ) [T], which allows spreading to be compared between experiments of different scales (e.g., columns to CPE). στ ð6Þ τ A moderate to high amount of spreading coupled with a low ratio of apparent to measured mobile water contents (or v bb 1; Eq. (4)), is suggestive of the existence of preferential flow paths and/or diffusive mass-transfer between mobile and immobile domains (assuming measured moisture contents significantly over-estimate actual mobile water contents). Finally, Eriksson et al. (1997) estimate the local longitudinal dispersivity λ [L] by fitting an advection-dispersion equation (ADE) (adapted from Kreft and Zuber, 1978) with a timeaveraged velocity to the measured flow-corrected breakthrough curve: CV τ ¼

0

1

2  2 3 τ C B 1− 7 6 Sðτ; zÞ 1 τ C B exp6 −  7 ¼ B  1 C 5 4   A @ λ τ Mo   3 2 4 τ 4π λz ττ z τ

ð7Þ

where Mo [M] is the applied tracer mass; the dispersion coefficient D [L2·T−1] (from Kreft and Zuber, 1978) is replaced by D ¼ λv, and v is the time-averaged pore water velocity [L·T−1] (Destouni et al., 1994). Three assumptions are associated with the use of Eq. (7): the observed tracer spreading is due only to local hydrodynamic dispersion or a Fickian model of transport (Roberts et al., 1986) and not diffusional solute exchange to immobile domains, the model fits assume a constant velocity and dispersion, and the use of first temporal moment values implicitly accounts for potential early tracer arrival due to preferential flow (Eriksson et al., 1997). The flow-corrected time approach and application of Eqs. (1) to (7) define the first and second temporal moment (Eqs. (2) and (5) or (6), respectively) and assist in identifying the existence of preferential flow paths (Eqs. (3) and (4)) and characterizing the dispersivity of the system, assuming no solute transfer to the immobile domain (Eq. (7)). 2.5. Numerical simulations Flow and solute (i.e., bromide) transport are modeled using HYDRUS1D, which includes several formulations for variably saturated flow in one dimension including single and dual domain approaches (Šimůnek and van Genuchten, 2008). For the former, variably saturated flow through a uniform, or single-

55

porosity, system is described by the Richards' equation (Eq. (8) (for one-dimensional systems); Richards, 1931; Šimůnek et al., 1998).

∂θðhÞ ∂ ∂h ¼ K ðhÞ þ1 ∂t ∂z ∂z

ð8Þ

where h is the pressure head [L], θ is the volumetric water content [L3·L−3], t is time [T], z is the spatial coordinate [L], and K(h) is the unsaturated hydraulic conductivity as a function of h [L·T−1]. Transport of a conservative solute is described by the classical advection-dispersion equation (Eq. (9)):

∂θc ∂ ∂c ∂qc − ¼ θD ð9Þ ∂t ∂z ∂z ∂z where c is the solution concentration [M·L−3], D is the dispersion coefficient accounting for hydrodynamic dispersion [L2·T−1], and q is the specific discharge [L·T−1]. Dual-porosity (including the MIM formulation) and dual permeability are two examples of a dual domain approach that may be modeled using HYDRUS1D (Šimůnek and van Genuchten, 2008). The selection of an appropriate model to describe solute transport through a porous medium is difficult to establish a priori (Roth et al., 1991). In an effort to maintain simplicity, the dualporosity MIM model is selected as the applied dual-domain approach for model simulations in this study. The dual-porosity approach assumes that the total water content θ of a system can be partitioned into its mobile θmo and immobile θim components (Eq. (10)). θ ¼ θmo þ θim

ð10Þ

In a dual-porosity model, either water and solute, or solute only, may transfer between the mobile and immobile domain. The former assumes flow and transport is restricted to the mobile domain, and both may transfer between the mobile and immobile porosities (Šimůnek and van Genuchten, 2008; van Genuchten and Wierenga, 1976). The latter, solute only transfer approach, is better described as a mobile–immobile (MIM) model in HYDRUS1D and assumes uniform water flow in the mobile domain, a constant immobile water content, and solute exchange between the two domains controlled by a masstransfer coefficient (ωmim; [T−1]). Two equations are required to describe water flow in the mobile domain (Eq. (11)) and moisture dynamics in the immobile domain (Eq. (12)) (Šimůnek et al., 2003):

∂θmo ðhmo Þ ∂ ∂hmo ¼ K ðhmo Þ þ 1 − Гw ∂t ∂z ∂z

ð11Þ

∂θim ðhim Þ ¼ Гw ∂t

ð12Þ

where Гw [T−1] is the mass-transfer rate for water between the mobile and immobile domains. Refer to Šimůnek et al. (2003) and Köhne et al. (2004) for alternative formulations of the masstransfer rate.

56

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

The dual-porosity conservative solute transport model adapts the standard advection-dispersion equation (Eq. (9)) to represent solute exchange between the two domains as a sum of an apparent first-order diffusion and dispersion process and advective transport (where applicable). Dual-porosity solute transport for a conservative solute is described below, with Eq. (13a) describing solute transport in the mobile domain, 13b is the mass balance for the immobile domain, and 13c is the transfer of solute mass between the mobile and immobile domains.



∂θmo ðcmo Þ ∂ ∂c ∂q c ¼ θmo Dmo mo − mo mo −Г s ∂t ∂z ∂z ∂z

∂θim ðcim Þ ¼ Гs ∂t  Г s ¼ ωmim ðcmo −cim Þ þ Г w c

ð13aÞ

ð13bÞ ð13cÞ

where cmo and cim are solute concentrations of the mobile and immobile regions [M·L−3], respectively, Dmo is the dispersion coefficient in the mobile region [L2·T−1], qmo is the specific discharge in the mobile region [L·T−1], ωmim is the masstransfer coefficient [T−1], and Гs is the mass-transfer term for solutes between mobile and immobile regions [M·L−3·T−1]. In the dual-porosity model, c* is equal to cmo or cim for Гw N 0 or Гw b 0, respectively. The dual-porosity MIM approach used in this study applies the uniform flow equation (Eq. (8)), where θ refers to the θm (or the mobile water content), and solute transport shown in Eq. (13). In the latter, Гw is zero as solute transfers between immobile and mobile porosities are restricted to diffusive mass transfer controlled by a mass-transfer coefficient (ωmim) and the concentration gradient between the two domains. 2.5.1. Model discretization and parameterization Model simulations in this study assume vertical flow (Section 2.2); therefore CPE simulations are conceptualized as a 10 m column consisting of 2 material types: 8 m of Class A (emplaced during two separate discharge events) overlying 2 m of Class C. The one-dimensional profile is discretized into 100 elements from the top surface of the pile to a lower seepage-face boundary (i.e., pressure head = atmospheric), with a uniform grid block size of 10 cm. Evaporation and subsequent infiltration at the upper surface of the CPE was computed from observed site precipitation (see Supporting data, Fig. S1) and the Hargreaves method was used for estimating evaporation (Jensen et al., 1997). The Hargreaves method was selected for two reasons: a water balance approach would add uncertainty by assuming the absence of non-vertical flow beyond the areal footprint (4 m × 4 m) of Lysimeter B and the Hargreaves evaporation estimation approach involves fewer input parameters (i.e., minimum/ maximum temperatures and latitude) than the Penman– Monteith method (i.e., net radiation, soil heat flux, wind speed, temperature, atmospheric pressure, relative humidity). Temperatures are measured hourly at a weather station, near the CPE site. Mean night-time and day-time temperatures are assumed to represent robust estimates of the minimum and maximum

daily temperatures, respectively, required for the Hargreaves evaporation equation. In the model, precipitation and evaporation are applied daily and a single instantaneous pulse of bromide is applied, representing the timing of the tracer test at the height of its first wet season. The model simulates flow and solute transport for the first 804 days following completion of construction, with minimum and maximum time steps of 0.005 days and 1 day, respectively. Daily input values for CPE simulations can be found in Supporting data (Table S1). The smaller-scale columns are described as 1D-vertical domains, 0.8 m in length, consisting of 3 materials; 2 × Class A (as Class A-1 and Class A-2) and Class C, all in equal proportions. The domain was discretized into 100 elements, with grid blocks of 8 mm, and a flux condition is assigned to the upper boundary. A seepage-face (i.e., h = 0) is applied to the lower boundary, which assumes the boundary flux will remain zero with a negative pressure head and outflow is triggered by the saturation of the lower soil profile. Daily precipitation and evaporation values were applied to column simulations. The former reflect area-normalized applied precipitation volumes. Evaporation was estimated using a water balance method (i.e., E = P − Q), by applying two assumptions. Firstly, all recorded outflow is assumed to be vertical and constrained by the spatial boundaries of the column geometry. Secondly, daily precipitation or rate of infiltration is month-dependent (to reflect seasonal changes at Antamina) and assumed to be constant within the days of that month. Therefore, the largest uncertainties with evaporation estimates occur at the first days of a month or when there is significant change in the daily applied precipitation rates. A single pulse of bromide was applied on day 241, to represent the start of the tracer test, and flow and solute transport was simulated for 372 days with similar minimum and maximum time steps as used in the CPE simulations. Daily input values for column simulations can be found in Supporting data (Tables S2 and S3, for C1 and C2 respectively). Six physical parameters are required to describe the matrix properties of Class A and C material used in the CPE and the laboratory columns; residual water content θr [L3·L−3], saturated water content θs [L3·L−3], dry bulk density ρb [M·L−3], saturated hydraulic conductivity Ksat [L·T−1], and two van Genuchten parameters α [L] and n [−]. The CPE was constructed over a period of 6–8 months and the material had a significant “wet-up” time. Estimates of water content θm [L3·L−3] for CPE simulations were based on measurement ranges from in situ volumetric water-content sensors (i.e., 0.18–0.36 m3∙m− 3) following the 2009 wet season (i.e., May–June 2009). Materials used for the laboratory columns were flushed and drained prior to construction to measure material-specific soil parameter input values (i.e., θr, θs) described in previous section (Table 2). Therefore, initial water contents were assumed to reflect residual water-content values. Initial estimates for dispersivity λ [L] values were obtained by fitting flow-corrected data to an ADE (Eq. (7)), using the method described previously by Eriksson et al. (1997), and are to be further refined in model simulations. The van Genuchten parameters were obtained using a fourstep procedure: 1) renormalizing PSD results to material passing the #4 sieve (i.e., b 4.75 mm); 2) apportioning the renormalized volumes to sand–silt–clay fractions (i.e., 4.75 mm N sand N

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65 2.0 Cumulave oulow (m 3·m-2)

0.42 mm N silt N 0.074 mm N clay); 3) applying these percentages to the USDA soil classification triangle (Brown, 2003; refer to Fig. S2 in Supporting data); and 4) referencing a table developed by Carsel and Parrish (1988) describing typical van Genuchten parameters for the 12 USDA-classified soil types. Reference to this table provided estimates for the shape parameters α and n (Eq. (14)), and saturated hydraulic conductivity Ksat values, which are used as initial input values for the simulations at both scales.

A

1.5

1.0 C1 C2 CPE

0.5

0.0 0

ð14Þ

3. Results 3.1. Measured flow and bromide breakthrough Flow monitoring for the CPE and the laboratory columns began in the dry season. In Fig. 3-A, the wet season (October– April) is reflected in rapid increases to CPE cumulative outflow (i.e., t = 150–300 days and 500–700 days), with significantly less flow and low-flow plateaus in the dry seasons (i.e., May–

100

200

300

400 500 Time (days)

600

700

800

B

300 [Bromide] (mg·L-1)

The #4 sieve size, described in step 1, was selected based on studies by Yazdani et al. (2000), which found porous media containing median grain diameters greater than 4.75 mm do not have significant capillarity in unsaturated conditions. It was inferred that material passing this sieve size actively contributes to the mobile and immobile domains hosting matrix flow paths and some preferential flow paths. Preferential flow related to non-capillary or film flow may be related to pore spaces generated by particles with greater than 4.75 mm. Final van Genuchten parameters (Table 2) used in the flow simulations were manually calibrated by matching simulated to measured outflow. Initially estimated Class A parameters (from laboratory measured and tables by Carsel and Parrish (1988)) are very similar to final model-calibrated values (shown in Supporting data, Table S4). Final Ksat values (3.0 m·day−1) for Class C were more closely related to the upper end of its soil classification (i.e., sandy loam (1 m·day− 1). This result is likely attributed to the significantly lower Class C component passing the 4.75 mm sieve. Note that saturated porosity estimates used in model simulations were obtained from values measured in the laboratory, not from typical soil property tables (i.e., θsat = 0.39–0.45). Dual-porosity models require estimates of immobile porosity θim [L3·L−3] and mass-transfer rates ωmim [T−1]. Immobile porosity values are defined as the portion of matrix material porosity where water is effectively stagnant. Initial estimates of immobile porosity were assumed to be comparable to laboratory-measured residual water content (shown in Table 2), but were refined with model calibration. It is important that residual moisture content values reflect the lowest possible immobile water content values. The mass-transfer coefficient ωmim was estimated by fitting simulated to measured results. Uniform flow and solute transport simulations were also carried out using the Richards' equation (Richards, 1931), for the purposes of comparing results from dual-porosity to single porosity model approaches. Input values for single porosity simulations were identical to those used in dual-porosity simulations (i.e., θsat, θr, λ, Kd).

400

200 CPE 100

0 200

400

[Bromide] (mg·L-1)

ðθs −θr Þ θðhÞ ¼ θr þ ð1 þ jαhjn Þm

57

300

400

500 Time (days)

600

700

800

C

300

200 C1 C2 100

0 200

250

300 Time (days)

350

400

Fig. 3. Cumulative area-normalized outflow (A) and bromide tracer concentration for CPE (B) and laboratory columns (C1 and C2; C) versus total experimental time. Arrows represent tracer application at 237 days (B) and 241 days (C).

September or t = 0–150 days, 300–500 days, and N700 days). Successful efforts to mimic the Antamina recharge cycle in the column experiments are supported by the similarity between column and CPE-recorded outflows over the total experimental time (Fig. 3-A). However, unlike the CPE that drains year-round, the column drainage occurs almost exclusively in wet season months. Initial column outflows begin 8 days into the wet season (i.e., t = 129 days) and no drainage reports to the base of the column over the majority of the dry season. Discrepancies between cumulative outflows from C1 and C2 are a result of difficulties encountered during experiment construction. Eaves were installed on the inner diameter of both columns to recover and remove rainfall along the column walls. These eaves were located at a lower level in C2 (i.e., closer to the upper boundary), thereby capturing ~10% more side-flow and contributing to lower cumulative outflow than C1. Bromide breakthrough from the CPE was bi-modal, with a high single-point peak shortly following application and a long release over the first 1.5 years at a much lower maximum concentration. Bromide was first measured 72 hours after tracer application at concentrations ranging from 0.2 mg∙L−1

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

to 7 mg∙L− 1. A single maximum concentration peak of 289 mg·L−1 Br− (C/C0 = 0.09) was observed at 9 days following tracer application (i.e., 245 days of total experimental time). A non-conservative tracer (uranine) was also applied concurrently with bromide for the purpose of aiding visual observation of rapid tracer breakthrough (Blackmore et al., 2012). The high, single-point bromide peak (i.e., 289 mg·L−1 Br−; t = 245 days or 9 days following tracer application) is coincident with the highest uranine concentration measured in the outflow from the waste rock pile (Blackmore et al., 2012), confirming that the early high concentration arrival of the tracer is not a measurement artifact. The rapid arrival of tracer at high concentrations provides strong evidence for preferential flow in this CPE; however, the bromide mass released at this peak concentration represents only a small proportion (i.e., 0.1%) of the total applied bromide mass (Fig. 3-B). The majority of bromide is released from the test pile over the first year, with peak concentration values of 130 mg·L−1. The subsequent breakthrough follows a Gaussian distribution from day 250 to 650 of the experimental duration (see Supporting data, Fig. S3). Measured tracer concentrations decrease (or tail) to detection limit values (i.e., ~1–2 mg·L−1) to the end of the experiment. Tracer tailing is generally attributed to mass exchange between regions of different mobility (Kissel et al., 1973). For the laboratory columns, bromide arrival occurred at 17 hours (C1) and 21 hours (C2) following tracer application at initially low concentrations (i.e., 4 mg∙L−1 Br−; Fig. 3-C). Maximum concentrations in C1 were 333 mg·L−1 Br− (C/C0 = 0.11) at 6 days and 307 mg·L−1 Br− (C/C0 = 0.10) in C2 at 9 days following tracer application (i.e., day-246 (C1) and day249 (C2) of total experimental time). In both C1 and C2, column breakthrough curves show trends of rapid release and long tails, similar to tailing observed in the CPE experiment.

Cumulative tracer recoveries for the CPE and laboratory columns, as a function of area-normalized flow and flow path length, are presented in Fig. 4. This method, as opposed to the commonly applied pore volume approach (i.e., cumulative area-normalized out flow multiplied by water content), was used due to uncertainties associated with measured absolute water content values (due to the requirement of placing sensors in relatively fine-grained material; Section 2.3) and to avoid the need of selecting a single effective water content value for highly heterogeneous materials. Results shown in Fig. 4 show the CPE and two columns released a similar amount of the applied bromide mass at the end of each experiment, with the CPE at ~80% and C1 and C2 at 72% and 75% (respectively). This result, in addition to similar area-normalized flow between scales (Fig. 3-A) and the conservative nature of bromide, suggests that solute transport is retarded or influenced in all experiments by similar processes and characteristics. Although both experiment scales present similar cumulative tracer recoveries, it is recognized that in terms of flow-path normalized time the CPE experiment released half of the applied tracer solute mass (i.e., M/Mo = 0.5) quicker than the column experiments (Fig. 4). This slightly faster response may be a result of the larger grain size distribution of materials used

0.08 Normalized mass flux rate (d-1)

58

A CPE

0.06

0.04

0.02

1.0

0 0

100

200

300

400

500

600

700

Flow corrected me (days)

0.8

Normalized mass flux rate (d-1)

0.08

M/Mo

0.6

0.4 C1 0.2

C2 CPE

0.0 0.0

0.2 0.4 0.6 Cumulave area and flow path normalized oulow (-)

B

C1 C2

0.06

0.04

0.02

0.8

Fig. 4. Cumulative tracer mass recoveries from CPE and column experiments versus normalized outflow. Tracer mass is normalized to the total mass of applied tracer (i.e., 5.4 g (C1); 4.0 g (C2), 1309 g (CPE)), whereas flow (m3) is normalized to the cross-sectional area (m2) (0.16 m2 C1/C2 versus 16 m2 CPE) and flow path length (m) (0.8 m C1/C2 versus 10 m CPE). Note x-axis reflects outflow flux following tracer application or flow recorded at t ≥ 241 days (C1/C2) and t ≥ 237 days (CPE).

0 0

20

40

60

80

100

120

Flow corrected me (days) Fig. 5. Bromide breakthrough curves for CPE (A) and column (B) experiments, in terms of normalized mass flux rate as a function of flow-corrected time (τ; Eq. (1)). Normalized mass flux rate is calculated as [S(τ,z)/Mo], where S(τ,z) e Note x-axis reflects time lapsed following tracer application, [M·T−1] = cðτ; zÞQ. (i.e., experimental time minus 237 days (CPE) or 241 days (columns)).

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

59

0.08 Normalized mass flux rate (d-1)

in the CPE experiment, relative to the column scales, which contributes to increased material heterogeneities and therefore a wider variety of transport regimes (Coppola et al., 2011). However, this difference did not significantly impact the dominant flow characteristics observed, which were similar in both experiments.

3.2. Temporal moment analysis

A) CPE

Measured λ = 0.095 m λ = 0.30 m

0.06

0.04

0.02

0 0

Normalized mass flux rate (d-1)

0.08

100

200 300 400 500 Flow-corrected me (d)

B) C1

600

700

Measured λ = 0.12 m λ = 0.15 m

0.06

0.04

0.02

0.00 0

20

40 60 80 Flow-corrected me (d)

100

120

0.08 Normalized mass flux rate (d-1)

Fig. 5 presents bromide breakthrough curves from all tracer experiments in terms of normalized mass flow rate S/Mo [T−1] and flow corrected time (τ) [T]. For the former, Mo [M] is the total tracer mass applied to the system and S(τ) was defined previously in Eq. (2). Flow-corrected breakthrough curves stretch and compress results during periods of high and low flow, respectively, and are calculated for time periods following tracer application. The CPE flow-corrected breakthrough curve (Fig. 5-A) shows that the majority of bromide release has occurred by 300 (flow-corrected) days (or day-545 in experimental time), followed by a significantly lower bromide release. The jagged line observed beyond 300 (flow-corrected) days in Fig. 5-A is reflective of decreasing bromide concentrations and increasing flow rates during the second wet season (2010–2011). Column flow-corrected breakthrough results (Fig. 5-B) are similar between C1 and C2. At the column-scale, tracer concentrations are only measured in the wet season as drainage does not report to the column base during the dry season. As such, column flow-corrected curves are relevant for the wet season only and the absence of periods of low flow and shorter flow paths (i.e., residence time) results in smoother curves relative to those observed at the CPE-scale. Data shown in Fig. 5 is used to calculate temporal moments (Table 3). For conservative tracer experiments, only the first and second temporal moments are used to estimate parameters such as the pore water velocity and dispersivity coefficient (Aris, 1958). The CPE experiment has a first temporal moment or mean residence time ðτ Þ of 295 days and a normalized second temporal moment or spreading (CVτ) of 0.20. Both column experiments (C1 and C2) have similar mean residence times and spreading (i.e., τ = 14 days and 20 days, CVτ = 0.53 and 0.49; respectively). Calculated apparent mobile water contents are similar for all experiments and v, the ratio of these values to in situ moisture content (Eq. (4)), is lower in the CPE. Lower calculated v values suggest a higher likelihood of preferential flow at the CPE scale, relative to the column scale. Low v values suggest a dual-domain approach, or immobile– mobile model conceptualization, may be applicable to all experiments.

C) C2

Measured λ = 0.08 m λ = 0.12 m

0.06

0.04

0.02

0.00 0

20

40 60 80 Flow-corrected me (d)

100

120

Fig. 6. Fitted ADE model results with measured normalized mass flux rate (as shown and described in Fig. 5). Model fits are calculated from Eq. (7).

3.3. Numerical modeling 3.3.1. Estimation of dispersivities using the flow-corrected ADE Dispersivity values estimated by fitting Eq. (7) to normalized mass flow rates are shown in Fig. 6. Dispersivity could not be estimated at the larger CPE-scale due to a poor agreement between theoretical (ADE-modeled) and measured breakthrough curves (Fig. 6-A). For column experiments, a moderate fit is observed between flow-corrected-time

Table 3 Calculation of effective flow and transport parameters.

CPE C1 C2

Water flow

First moment

Second moment

Normalized second moment

Apparent mobile water content

e ðτ Þ (L·day−1) Q

τ (days)

στ (days)

CVτ (−)

 m

38 0.72 0.70

295 14 20

59 7.6 9.9

0.20 0.53 0.49

0.07 0.08 0.11

(m3·m−3)

Ratio of apparent to actual mobile water content va (%) 32–41 40–44 53–58

a Minimum–maximum range of measured moisture content values (θ) used in calculating v values for CPE and columns; i.e., CPE [0.17 m3·m−3–0.22 m3·m−3], C1 [0.18 m3·m−3–0.20 m3·m−3] and C2 [0.19 m3·m−3–0.21 m3·m−3].

60

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

breakthrough curves and modeled ADE curves and dispersivity values are estimated to be 0.12–0.15 m (C1) and 0.08–0.12 m (C2) (Fig. 6-B and C). However, modeled ADE curves could not simulate peak concentrations (at τ = 10 days) for C1 or breakthrough tails (between τ = 19 days–40 days) for C2 (Fig. 6-B and C, respectively). This result suggests that a single porosity ADE may be able to partially simulate solute transport at the smaller scale, which is not possible in an increasingly complex and larger scale system. The ADE model assumes that the flow regime (advection and dispersion) is temporally invariant, which is applicable to the column scale due to its smaller size and one-half wet season was required to flush the majority of the tracer mass. At the CPE-scale, several wet-dry seasons were required to flush a similar tracer mass and the change in flow regimes (between seasons) cannot be accommodated well by the ADE model.

3.3.2. Flow and solute transport modeling using the MIM approach Numerical simulations implemented the mobile–immobile (MIM) dual-porosity approach, following the inability to simulate transport at the CPE-scale and inadequacy to simulate breakthrough peaks (C1) and tails (C2) using a uniform ADE. Simulations held the immobile water contents constant and only solute transfer occurs between the mobile and immobile domains (Figs. 7 and 8). A good agreement was observed between measured and simulated solute breakthrough at the CPE- and column-scale (Figs. 7 and 8). Specifically, measured breakthrough peaks and tails from column experiments are captured with the MIM approach. In regards to the CPE, the dominant flow and solute transport are captured with the MIM approach; however simulations did not capture the single-point peak concentration in CPE breakthrough (at 9 days after tracer release corresponding to day 245 since experiment initiation). Although CPE transport

400

A) CPE

2.0

A) CPE 300

1.6

[Bromide] (mg·L-1)

Cumulave oulow (m)

2.4

1.2 0.8 Measured Modeled

0.4 0.0

100

0

0

200

400 Time (Days)

600

800

200

1.2 1.0

400

500 Time (days)

600

700

800

B) C1

0.8 0.6 0.4

300

200

100

0.2 0

0.0 0

100

200 Time (days)

300

200

400

250

300 Time (days)

350

400

400

1.0

C) C2

C) C2 0.8

[Bromide] (mg·L-1)

Cumulave Oulow (m)

300

400

B) C1 [Bromide] (mg·L-1)

Cumulave Oulow (m)

200

0.6 0.4

300 Measured Dual-Porosity Model Uniform Model

200

100

0.2 0

0.0

200

0

100

200

300

400

250

300 Time (days)

350

400

Time (days) Fig. 7. Measured (solid lines) versus simulated (dashed lines) flow using a mobile–immobile model. Time represents total experimental time.

Fig. 8. Measured (symbols) versus simulated (lines) solute transport using a dual-MIM model (dashed lines) and a best-fit uniform (dotted lines) modeling approach. Time represents experimental time.

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

may be better described by a model that includes a small fraction of open void or film flow, more complex simulations were not carried out as this peak represents a small flow component in relation to total flow (i.e., 0.1%). Secondary simulations that assume uniform flow and solute transport (i.e., uniform model) were performed to validate the efficacy and level of simplicity required to simulate observed breakthrough curves from all experiments. At the CPE-scale, the uniform model predicted tracer arrival only after a near complete flushing of the tracer from the pile. This result confirms that a uniform flow/solute model is inadequate at the CPE scale and supports the previous inabilities to model CPE breakthrough with a single-porosity flow-corrected ADE. Results for the column simulations show that the uniform model is able to reproduce tracer arrival in general terms (Fig. 6B and -C; fine dashed lines); however, this model predicted late arrival and reduced peak concentrations. Differences between simulated breakthrough curves obtained from uniform ADE versus MIM approaches are likely attributed to the use of a timeaveraged flow velocity versus transient flow (respectively). Nevertheless, these results suggest that the column experiments can be simulated, in a general sense, with a single domain ADE or uniform model approach; however a dualporosity MIM model better simulates the observed peak concentrations and breakthrough tails. For the CPE, only the dual porosity MIM approach yields a fit for measured breakthrough curves. Model simulations for column experiments used identical flow parameters to simulate columns C1 and C2, which further supports the previous observations (from Fig. 4 and Fig. 5-B) that both columns contain material with comparable characteristics and transport regimes. Final dispersivity estimates, immobile porosity and mass-transfer coefficient values used in the CPE and column simulations are given in Table 4. Masstransfer coefficients of Class A material differed by a factor of 3 to 30 between experimental scales, whereas immobile porosity values were similar. Higher dispersivity values were applied in the MIM model simulations, in comparison to those estimated using a singleporosity flow-corrected ADE (shown in Fig. 6). A similar observation was noted by Sassner et al. (1994), who suggested that estimated dispersivity values from ADE-modeled curves may indeed be larger. Weighting of dispersivity values shown in Table 4 to flow path length (e.g., CPE: Class A 75%: Class C 25%) produces a CPE dispersivity values that is approximately two times higher than those used in column simulations (i.e., CPE: 0.25 m versus C1/C2: 0.14 m), which agrees with well-established concepts of scale-dependency of dispersivity (Coppola et al., 2011; Vanderborght and Vereecken, 2007). Changes to final dispersivity values were applied in several sensitivity analyses and did not substantially impact the results.

61

Model simulations of solute transport in the CPE could not capture the initial tracer peak (i.e., 289 mg·L−1 at day-245 or 9 days following tracer application); however, this single-point peak represents only ~0.1% of the captured breakthrough mass. This minor component likely reflects fast preferential or film flow as a result of more heterogeneous, coarser-grained material present at the larger scale. Although the dual-porosity MIM approach was not able to reproduce the fast preferential flow component on the CPE scale, the model results adequately captured fast matrix flow and solute mass-transfer to slower matrix flow and immobile waters at both scales and should be considered non-unique.

4. Discussion This study compares flow and tracer results at two scales, differing by an order of magnitude in height (i.e., 10 m versus 0.8 m). The objective of the experimental design/construction was to ensure that similar materials and conditions were used in all experiments and reduce the likelihood of discrepancies between scales. PSD curve similarities (refer to Fig. 2) support the comparability of study materials used, as per a ≤10% difference between the percent passing the 2 mm sieve diameter for CPE and column experiments. The conclusion that these efforts were adequate in replicating one system at two scales was supported by 3 observations: similar cumulative outflow (Fig. 3) and measured tracer mass release (Fig. 4) at both scales, and the application of identical or similar hydrologic parameters (i.e., α, n, Ksat, θs, θr, θim) in model simulations, with the exception of dispersivities and, in particular, mass transfer coefficients. CPE breakthrough curves show that the majority of bromide is released over the first year (following tracer application), at a mean residence time of 295 days, and preferential flow was evidenced as a minor (~0.1%) flow component due to a single maximum concentration bromide peak within days of application. Breakthrough curves from C1 and C2 were similar, with both showing a rapid rise in drainage concentrations to maximum values, a decrease in measured concentrations with increasing time, but no indications of preferential flow. Results depicted in Fig. 4 show that cumulative tracer recoveries (normalized to cross sectional area and flow path length) from all experiments are comparable. This comparability suggests both experiments are representative of a similar transport scale, or an averaging of observation (or smaller) scale heterogeneities within the porous media (Roth et al., 1991). This result indicates that this study was successful in mimicking processes active at the larger scale in the laboratory and demonstrates the applicability of estimated solute parameters from column experiments to field scale experiments

Table 4 Final flow and solute transport parameters used in model simulations. CPE Parameter λmo θim ωmim

[m] [m3·m−3] [day−1]

Laboratory column

Class A

Class C

Class A-1

Class A-2

Class C

0.30 0.175 4.5 × 10−5

0.10 0.043 5.0 × 10−8

0.12 0.16 1.4 × 10−4

0.17 0.22 1.4 × 10−3

0.12 0.06 5.0 × 10−8

62

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

(Coppola et al., 2011) — with the exception of dispersivities and mass transfer coefficients. The temporal scale of tracer breakthrough differed between the two experiments. Two water years (or wet-dry seasons) were required for the CPE to release approximately 81% of the total applied bromide mass. At the column scale, one-half wet season resulted in the release of 74% (C1) and 72% (C2) of the solute mass. Although temporal differences are noted, bromide continued to be detected in all drainage waters at concentrations near detection limit values (i.e., 1–2 mg·L−1) until the end of the experiment duration. These long tails suggest extremely slow matrix flow in zones of lower conductivity and/or masstransfer from stagnant or immobile domains is present in both experiments. The validity of the flow-corrected time method to calculate temporal moments was evaluated by comparing flow-corrected breakthrough curves to modeled curves using a single-domain advection dispersion equation with homogeneous parameters (i.e., Eriksson et al., 1997; modified from Kreft and Zuber (1978)). A moderate fit between column-scale flow-corrected data (C1 and C2) and the ADE-modeled curves indicates that temporal moments may be estimated as single effective values and implicitly account for the existence of variable flow paths and macrodispersive solute spreading (Eriksson et al., 1997). This agreement suggests a single-domain uniform model may be applicable at the column scale. However, only transient flow models that applied the dual porosity MIM approach were able to capture the breakthrough curve and tailing observed from the column experiments. These results indicate at least one of the assumptions associated with the flow-corrected time method (i.e., effective first temporal moments or macrodispersion) is not valid. The same single-domain ADE could not simulate flowcorrected breakthrough at the CPE scale. This result was not surprising as aquifer heterogeneity is typically proportional to experimental scale and the ability to calculate parameters as single or effective values is increasingly difficult in highly heterogeneous systems. Additionally, the CPE breakthrough curve extends over two water years (Fig. 3) that alternates between dominantly drain-down matrix flow (dry seasons) and event flow including preferential and fast matrix flow (wet seasons). This large temporal variation in flow regimes may not be adequately captured by the flow-corrected time compression and add to the difficulty in resolving effective parameters using this method. Breakthrough curves that present long tails and tracer mass recoveries less than unity (i.e., CPE = 0.80; C1 = 0.75; C2 = 0.72) suggest that a rate-limited mass-transfer model may better explain solute transport at both experimental scales and support the application of a MIM approach. The application of this approach does not infer the absence of macrodispersive processes at either scale, but points to its suitability in representing observed bromide breakthroughs. In this approach, waste rock can be conceptualized as containing measurementscale mobile and immobile domains with variable conductivities that are associated with pore-scale mass-transfer processes (Harvey and Gorelick, 2000). Solute mass is partitioned between either the moving (or mobile) fluid and a less permeable (or immobile) domain (Carrera et al., 1998; Haggerty and Gorelick, 1995; Haggerty et al., 2000) and solute movement between domains is defined by a mass-transfer coefficient (ωmim; day−1

(Eq. (13c))) and occurs by physical or chemical processes (Haggerty et al., 2004). In this study, mass-transfer coefficient values were estimated by fitting the simulated dual-porosity MIM curves to observed breakthrough tails (Table 4) and final values present the most significant identifiable difference between column and CPE experiments. Specifically, Class A material mass-transfer coefficients were approximately an order of magnitude smaller at the CPE scale, relative to values estimated for the columns (i.e., ωcolumn = 0.14 -1.4 × 10−3 day−1; ωCPE = 4.5 × 10−5 day−1). The application of the same conservative, non-sorbing tracer in all experiments indicates that chemical processes are ruled out as a contributing factor to the difference in mass-transfer coefficient values. Additionally, the similarity of immobile and mobile porosities, (area-normalized) outflow, and (flow-path normalized) tracer recovery suggests that physical processes do not differ considerably between the two scales. Temporal and/or spatial dependencies are two possible explanations for the mass-transfer coefficient difference. For the former, temporal dependencies may be related to differences in experiment duration. This concept is described by Haggerty et al. (2004), which showed that mass-transfer coefficients, estimated by fitting the breakthrough curve tail, are representative of mass-transfer processes contributing to the observed tail and operating over the experimental duration. In this scenario, the timescale of mass transfer estimated from a single-rate model will increase with the experiment duration and produce lower mass-transfer coefficients (Haggerty et al., 2004). Moreover, breakthrough curve tails from highly heterogeneous systems, such as materials used in this study, may not be captured adequately by a single-rate mass-transfer model. Therefore single or effective rates likely reflect an averaging of multiple mass-transfer rates as well as the timescale and, given the larger flow path/duration of the CPE experiment, may result in lower mass-transfer coefficient estimates relative to column experiments. Mass-transfer spatial dependences in a dual-porosity approach refer to pore-scale solute transfer process between immobile and mobile domains. Two lines of evidence suggest that the smaller-scale column experiment presents a greater likelihood of contact between immobile and mobile domains; larger flow-corrected (normalized) second moments (CVτ) and v values, relative to CPE values (refer to Eqs. (6) and (4) (respectively) and Table 3). Larger second moments are suggestive of possible preferential flow and diffusive mass-transfer between mobile and immobile domains (Eriksson et al., 1997). Measured column breakthrough curves demonstrate that preferential flow, or the concentration of flow into spatially localized areas (Nichol et al., 2005), was a negligible component of the overall flow regimes at this scale and the majority of tracer solute infiltrated the matrix pore-space. However, the larger mass-transfer coefficient from column simulations does complement these larger second moments relative to the CPE scale. Calculated v values are ratios of estimated apparent mobile water contents (θm⁎; m3∙m−3) and total moisture content held in the granular matrix. Ratios represent the proportion of water associated with mobile pore water and were calculated to be 40–44% (C1), 53–58% (C2), and 30–40% (CPE). The remaining proportion (1 – v; [−]) is associated with slowly moving or immobile pore water and, as such, results indicate that columns have a larger "mobile-to-immobile" porosity ratio and a greater

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

likelihood for contact between advecting solutes and immobile domains. 4.1. Parameter — estimation uncertainties and future prediction All input parameters included in dual-porosity simulations are considered to be independently constrained and within reasonable boundaries of each experiment. The calculation of temporal moments using the flow-corrected method is shown to be reasonably robust at the smaller scale, whereas larger scales present a higher likelihood for spatially heterogeneous flow behavior. A dual-porosity and rate-limited mass-transfer approach is shown to be better suited to simulating measured breakthrough response. This result suggests temporal moment calculations using the Eriksson et al. (1997) method may be more suited to conceptualize, as opposed to define, likely flow regimes and estimate parameter values. The assumption that measured water contents (θmeas) are equal to actual mobile water contents (refer to Section 2.4) may introduce some uncertainty to the calculation of v values. Measured water contents from in situ sensors may overestimate mobile water contents for two possible reasons: sensors used in this study were purposely placed in finer-grained material to ensure adequate contact between infiltrating waters and sensor surfaces and may not be representative of the true grain size distribution of the surrounding porous medium; and, water may be immobilized or cling to waste rock surfaces. These factors imply that sensor-measured water contents may be greater than the true average water content, which was also highlighted in Eriksson et al. (1997). However, as was mentioned previously, the comparability of materials used in both experiments suggests overestimations of water content will be present at both scales to a similar degree. Therefore, calculations of the proportion of preferential flow (i.e., v values) are better regarded as qualitative estimates and assist with conceptualization of flow domains. Simulation results were insensitive to grid refinement, as finer grid blocks (i.e., 2.5 cm (CPE) and 2 mm (C1 and C2)) produced nearly identical results as those shown in Figs. 7 and 8. To gain insight into the length of time required to flush tracer from the system, C1 simulations were extended to 2× and 4× the duration of this experiment (i.e., t = 744 days and 1488 days

1000 Measured Base Case

100 [Bromide] (mg·L-1)

t = 2X t = 4X 10

1

0.1

0.01 200

400

600

800 1000 Time (days)

1200

1400

Fig. 9. Multi-year simulations (dashed lines) of C1 bromide breakthrough. Note, measured and base-case results reflect the first year of analyses and are shown previously in Fig. 8-B. Time refers to total experimental time, or time lapsed following experiment initiation.

63

(respectively)) (Fig. 9). These simulations assume identical yearly recharge and input parameters (as the base-case scenario (i.e., t = 372 days)). Results indicate that tracer discharge will occur for time periods well beyond this experiment and at concentrations at or below minimum measureable limits. Similar results (not shown) were observed from simulations of C2 and the CPE. As such, tracer mass not shown in Fig. 4 will likely be released for many years to decades for 100% recovery. 5. Conclusions Outflow data and breakthrough curves from tracer tests were collected from two laboratory columns (C1 and C2) and an experimental waste rock pile (CPE) at the Antamina mine. The experiments contain similar waste rock assemblages and grain size distributions for fractions with diameters less than 10 cm, but differ in their vertical scale (i.e., h = 0.8 m versus h = 10 m (respectively)). Rainfall rates applied to the laboratory columns mimic the seasonal rainfall pattern at the mine site. Tracer recovery in the laboratory columns, excluding a long, lowconcentration tail, was complete in approximately 40 days; a similar tracer recovery in the CPE required approximately 400 days, approximately proportional to the length of the flow path. The following conclusions are highlighted: 1. Cumulative tracer recovery, normalized to outflow area and flow path, from the CPE and columns is similar. This result indicates that transport regimes are comparable between experiments and the larger scale system was reasonably replicated at a smaller, more manageable scale. As well, this comparability suggests that the presence of smaller-scale heterogeneities (that may influence solute transport) at both experimental scales have been smoothed and solute transport is governed by similar processes. 2. While the CPE indicates the presence of rapid flow paths (open void or film flow), this flow mechanism was a minor component of the overall mass-transfer in this experiment (~ 0.1%). The majority of the tracer mass, in all experiments, migrated through the granular matrix materials located between the coarser fractions (e.g., gravel, cobbles, and boulders). 3. A rate-limited mass-transfer model is suggested to explain observed solute transport and is supported by three observations: long-tail breakthrough curves, tracer mass fractions less than 1 (i.e., 0.72–0.80), and the application of a dualporosity MIM model to simulate measured results. The masstransfer model can better simulate the large physical differences in waste rock materials used in both experimental scales. 4. Dual-porosity MIM models applied similar or identical input values at both scales, with the exception of the scaledependent solute mass-transfer coefficient and dispersivity values. Column experiments yielded a mass-transfer coefficient upwards of an order of magnitude higher than the value applied in CPE simulations. Spatial and temporal reasons exist for a higher mass-transfer coefficient values at the columnscale. On a temporal scale, mass-transfer coefficients are estimated by fitting a breakthrough curve tail and are a function of experimental duration (i.e., longer temporal scales may yield lower mass-transfer coefficients). At a spatial scale, materials with strongly contrasting physical properties, such

64

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65

as those used in this study, are generally better represented using a multiple-rate mass-transfer approach. Therefore, a single- or effective-rate mass-transfer coefficient may reflect an averaging of the range of mass-transfer values encountered along the experimental flow path. Column experiments also yielded larger secondary temporal moments (CVτ) values and larger ratios of "mobile to immobile" porosities (v), relative to those values from the CPE. This supports the likelihood of greater transfer and/or contact between the two domains in the column experiments and a higher mass transfer coefficient, relative to the CPE. Although results from model simulations should be considered non-unique, differences in input and model calibrated parameters assist with conceptualizing unsaturated flow and transport at both scales. 5. In this study, the application of a flow-corrected method was best utilized as a qualitative estimation or conceptualization of the flow regimes at both scales. Specifically, the heterogeneity of materials and the scale discrepancy between experiments make it difficult to calculate robust effective or singleparameter values. A single porosity flow-corrected ADE model was able to moderately simulate breakthrough curves at the column scale only, but solute transport was better simulated using a MIM approach. This result suggests that large temporal variation in flow regimes, particularly at the larger scale, may not be adequately captured by the flowcorrected time compression and necessitate a more complex model approach. As well, several flow-corrected terms (e.g., θm⁎, v) are associated with values that are obtained in situ, which may introduce parameter uncertainty. This study employed a systematic method to simulate flow regimes at laboratory and field environments by stressing the importance of maintaining material homogeneity between experiments. Characterization of waste rock material, in regards to its PSD, is critical to evaluating the representativeness of scaling up results from manageable experiments. Numerical models incorporated identical material-specific physical property values (e.g., n, Ksat) for all experiment simulations as per similar physical properties. The observed similarity of column and CPE breakthrough curves from this study does not assume a precise replication of flow regimes between experimental scales, as these results are a function or summation of all possible complexities encountered within the unsaturated porous media. Instead, the comparability of these results signifies that it is possible to mimic field flow conditions in a laboratory-based, or more manageable, environment and enables a better understanding of the governing controls on solute transport regimes in larger systems. The results of this study may be applied to other locations that use waste rock piles as a method of mine waste deposition. As was mentioned previously, the prediction of flow/solute transport from waste rock piles is a requirement for most mine sites and the determination of a suitable modeling approach is largely dependent on the homogeneity of the porous media (Coppola et al., 2011). Typically, waste rock consists of highly heterogeneous materials both from a mineralogical and physical perspective. Therefore the ability for this study to apply a reasonably simplistic modeling approach to simulate the dominant flow and solute transport characteristics from heterogeneous waste rock suggests that the methods described therein may be equally applicable to other mine sites, and possibly also other

unsaturated sedimentary deposits composed of gravelly and stony media. Acknowledgements Partial funding of this project was provided by the Natural Sciences and Engineering Research Council (NSERC) of Canada and a Canadian Foundation for Innovation (CFI) grant awarded to Dr. K.U. Mayer. The authors would like to thank the personnel at the Antamina mine (Antamina Peru) that helped make this project a success and to H. Peterson and M. Gupton for their assistance in waste rock tracer tests. As well, the authors would like to thank the Spectroscopy and Kinetics Hub, in the Laboratory of Molecular Biophysics (University of British Columbia; Vancouver, BC, Canada), for access and training in tracer solute analyses. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx.doi.org/10.1016/j.jconhyd.2014.10.009. References Antamina, 2007. Ore Control Procedure Manual. Companía Minera Antamina S.A., Geology Department, Yanacancha, Peru. Aris, R., 1958. On the dispersion of linear kinematic waves. Proc. R. Soc. Lond. A Math. Phys. Sci. 245 (1241), 268–277. http://dx.doi.org/10.1098/rspa.1958. 0082. Bay, D.S., Peterson, H.E., Singurindy, O., Aranda, C.A., Dockrey, J.W., Sifuentes Vargas, F., Mayer, K.U., Smith, J.L., Klein, B., Beckie, R.D., 2009. Assessment of Neutral pH Drainage from Three Experimental Waste-rock Piles at the Antamina Mine, Peru. Securing the Future, 8th International Conference on Acid Rock Drainage. Skellefteå, Sweden. Blackmore, S.R., Speidel, B., Critchley, A., Beckie, R.D., Mayer, K.U., Smith, L., 2012. The Influence of Spatial Heterogeneity and Material Characteristics on Fluid Flow Patterns through Two Unsaturated, Mixed Waste-rock Piles. 9th International Conference on Acid Rock Drainage. Ottawa, Ontario. Brown, R.B., 2003. "Soil Texture." Fact Sheet SL-29. Soil and Water Science Department, Florida Cooperative Extension Service. Institute of Food and Agricultural Sciences, University of Florida, September 2003. 8 pages. Brusseau, M.L., Rao, P.S.C., 1990. Modeling solute transport in structured soils: a review. Geoderma 46 (1–3), 169–192. http://dx.doi.org/10.1016/00167061(90)90014-Z. Butcher, B., Hinz, C., Flury, M., Fluhler, H., 1995. Heterogeneous flow and solute transport in an unsaturated stony soil monolith. Soil Sci. Soc. Am. J. 59 (1), 14–21. Carrera, Jesús, Sánchez-Vila, Xavier, Benet, Inmaculada, Medina, Agustín, Galarza, Germán, Guimerà, Jordi, 1998. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeol. J. 6 (1), 178–190. http://dx.doi.org/10.1007/s100400050143. Carsel, Robert F., Parrish, Rudolph S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24 (5), 755–769. http://dx.doi.org/10.1029/WR024i005p00755. Christiansen, J.E., 1942. Irrigation by Sprinkling, (Bulletin No. 670, Berkley). Coppola, Antonio, Comegna, Alessandro, Dragonetti, Giovanna, Dyck, Miles, Basile, Angelo, Lamaddalena, Nicola, Kassab, Mohamed, Comegna, Vincenzo, 2011. Solute transport scales in an unsaturated stony soil. Adv. Water Resour. 34, 747–759. http://dx.doi.org/10.1016/j.advwatres.2011.03.006. Corazao Gallegos, J.C., 2007. The Design, Construction, Instrumentation and Initial Response of a Field-scale Waste Rock Test Pile(Masters' thesis) University of British Columbia. Destouni, Georgia, 1991. Applicability of the steady state flow assumption for solute advection in field soils. Water Resour. Res. 27 (8), 2129–2140. http:// dx.doi.org/10.1029/91WR01115. Destouni, Geogia, Sassner, Mona, Jensen, Karsten Høgh, Destouni, Geogia, Sassner, Mona, Jensen, Karsten Høgh, 1994. Chloride migration in heterogeneous soil — 2. Stohastic modeling. Water Resour. Res. 30 (3), 747–758. Eriksson, Nils, Gupta, Archana, Destouni, Georgia, 1997. Comparative analysis of laboratory and field tracer tests for investigating preferential flow and transport in mining waste rock. J. Hydrol. 194 (1–4), 143–163. http://dx. doi.org/10.1016/S0022-1694(96)03209-X.

S. Blackmore et al. / Journal of Contaminant Hydrology 171 (2014) 49–65 Fala, O., Molson, J., Aubertin, M., Bussière, B., 2005. Numerical modelling of flow and capillary barrier effects in unsaturated waste rock piles. Mine Water Environ. 24 (4), 172–185. http://dx.doi.org/10.1007/s10230-005-0093-z. Flury, Markus, Flühler, Hannes, Jury, William A., Leuenberger, Jörg, 1994. Susceptibility of soils to preferential flow of water: a field study. Water Resour. Res. 30 (7), 1945–1954. http://dx.doi.org/10.1029/94WR00871. Frostad, S., Klein, B., Lawrence, R.W., 2005. Determining the weathering characteristics of a waste dump with field tests. Int. J. Surf. Min. Reclam. Environ. 19 (2), 132–143. http://dx.doi.org/10.1080/13895260500157988. Gerke, H.H., van Genuchten, M.T., 1993a. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res. 29 (2), 305–319. http://dx.doi.org/10.1029/ 92WR02339. Gerke, H.H., van Genuchten, M.T., 1993b. Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resour. Res. 29 (4), 1225–1238. http://dx.doi.org/10.1029/92WR02467. Golder, 2010. Waste-rock and Tailings Geochemistry — Field Cell Monitoring, Dec 2002 to May 2009”. Antamina, Peru. Haggerty, Roy, Gorelick, Steven M., 1995. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour. Res. 31 (10), 2383–2400. http://dx.doi.org/10.1029/ 95WR10583. Haggerty, Roy, McKenna, Sean A., Meigs, Lucy C., 2000. On the late-time behavior of tracer test breakthrough curves. Water Resour. Res. 36 (12), 3467–3479. http://dx.doi.org/10.1029/2000WR900214. Haggerty, Roy, Harvey, Charles F., von Schwerin, Claudius Freiherr, Meigs, Lucy C., 2004. What controls the apparent timescale of solute mass transfer in aquifers and soils? a comparison of experimental results. Water Resour. Res. 40 (1), W01510. http://dx.doi.org/10.1029/2002WR001716. Harvey, Charles, Gorelick, Steven M., 2000. Rate-limited mass transfer or macrodispersion: which dominates plume evolution at the Macrodispersion Experiment (MADE) site? Water Resour. Res. 36 (3), 637–650. http://dx.doi.org/10.1029/1999WR900247. Herasymuik, G.M., 1996. Hydrogeology of a Sulfide Waste Rock Dump(Masters' thesis) University of Saskatchewan. Jensen, D.T., Hargreaves, G.H., Temesgen, B., Allen, R.G., 1997. Computation of ETo under nonideal conditions. J. Irrig. Drain. Eng. ASCE 123 (5), 394–400. Jury, William A., 1982. Simulation of solute transport using a transfer function model. Water Resour. Res. 18 (2), 363–368. http://dx.doi.org/10.1029/ WR018i002p00363. Jury, William A., Flühler, Hannes, 1992. Transport of chemicals through soil: mechanisms, models, and field applications. In: Sparks, Donald L. (Ed.), Advances in Agronomy vol. 47. Academic Press, pp. 141–201. Kellagher, Richard C., Flanagan, Francis James, 1956. The multiple-cone sample splitter. J. Sediment. Res. 26 (3), 213–221. http://dx.doi.org/10.1306/ 74D705E2-2B21-11D7-8648000102C1865D. Kissel, D.E., Richie, J.T., Burnett, E., 1973. Chloride movement in undisturbed clay soil. Soil Sci. Soc. Am. Proc. 37, 21–24. Köhne, J. Maximilian, Köhne, Sigrid, Mohanty, Binayak P., Šimůnek, Jirka, 2004. Inverse mobile–immobile modeling of transport during transient flow effects of between-domain transfer and initial water content. Vadose Zone J. 3 (4), 1309–1321. http://dx.doi.org/10.2113/3.4.1309. Kreft, A., Zuber, A., 1978. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chem. Eng. Sci. 33 (11), 1471–1480. http://dx.doi.org/10.1016/0009-2509(78)85196-3. Li, Yimin, Ghodrati, Masoud, 1997. Preferential transport of solute through soil columns containing constructed macropores. Soil Sci. Soc. Am. J. 61, 1308–1317. http://dx.doi.org/10.2136/sssaj1997.03615995006100050004x. Malmström, Maria E., Destouni, Georgia, Banwart, Steven A., Strömberg, Bo H.E., 2000. Resolving the scale-dependence of mineral weathering rates. Environ. Sci. Technol. 34 (7), 1375–1378. http://dx.doi.org/10.1021/es990682u. McKeowen, R., Barbour, S.L., Rowlett, D., Herasymuik, G., 2000. Characterization of the Grain-size Distribution for Waste Rock from Metal Mines — A Review of Existing Grain Size Data and an Evaluation of the Implications for Hydrogeologic Behaviour”. Proceedings of the 28th Annual Conference of the Canadian Society for Civil Engineerings. London, Ontario. Morin, K.A., Gerencher, E., Jones, C.E., Konasewich, D.E., 1991. Critical Literature Review of Acid Drainage from Waste Rock. Report 1.11.1. MEND. Nichol, Craig, Smith, Leslie, Beckie, Roger, 2005. Field-scale experiments of unsaturated flow and solute transport in a heterogeneous porous medium. Water Resour. Res. 41 (5). http://dx.doi.org/10.1029/2004WR003035 n/a–n/a.

65

Otwinowski, M., 1995. Scaling Analysis of Acid Rock Drainage. Mine Environment Neutral Drainage Program. Peterson, H.E., Beckie, R.D., Harrison, B., Mayer, K.U., Smith, L., 2012. Rapid Seasonal Transition from Neutral to Acidic Drainage in a Waste Rock Test Pile at the Antamina Mine, Peru. 9th International Conference on Acid Rock Drainage. Ottawa, Ontario. Presley, B.J., 1971. Techniques for analyzing interstitial water samples, part I: determination of selected minor and major inorganic constituents. Initial Reports of the Deep Sea Drilling Project vol. VII. U.S. Government Printing Office, Washington, DC. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. J. Appl. Phys. 1 (5), 318–333. http://dx.doi.org/10.1063/1.1745010. Roberts, Paul V., Goltz, Mark N., Mackay, Douglas M., 1986. A natural gradient experiment on solute transport in a sand aquifer: 3. Retardation estimates and mass balances for organic solutes. Water Resour. Res. 22 (13), 2047–2058. http://dx.doi.org/10.1029/WR022i013p02047. Roth, K., Jury, W.A., Fluhlur, H., Attinger, W., 1991. Transport of chloride through an unsaturated field soil. Water Resour. Res. 27 (10), 2533–2541. Sassner, Mona, Jensen, Karsten Høgh, Destouni, Georgia, 1994. Chloride migration in heterogeneous soil — 1. Experimental methodology and results. Water Resour. Res. 30 (3), 735–745. Šimůnek, Jirka, van Genuchten, Martinus Th., 2008. Modeling nonequilibrium flow and transport processes using HYDRUS. Vadose Zone J. 7 (2), 782. http://dx.doi.org/10.2136/vzj2007.0074. Šimůnek, J., Šejna, M., van Genuchten, M.T., 1998. The HYDRUS-1D Software Package for Simulating One-dimensional Movement of Water, Heat, and Multiple Solutes in Variably Saturated Water. Version 1.0. Colorado School of Mines, Goldon: IGWMC-TPS-70, International Groundwater Model Center. Šimůnek, Jirka, Jarvis, Nick J., van Genuchten, M.Th., Gärdenäs, Annemieke, 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272 (1), 14–35. http://dx.doi.org/10.1016/S0022-1694(02)00252-4. Small, Mitchell J., Mular, John R., 1987. Long-term pollutant degradation in the unsaturated zone with stochastic rainfall infiltration. Water Resour. Res. 23 (12), 2246–2256. http://dx.doi.org/10.1029/WR023i012p02246. Smith, L., Beckie, R.D., 2003. Hydrologic and geochemical transport processes in mine waste rock. Environmental Aspects of Mine WastesShort Course SeriesMineralogical Association of Canada, Ottawa. Strömberg, Bo, Banwart, Steven, 1994. Kinetic modelling of geochemical processes at the Aitik mining waste rock site in northern Sweden. Appl. Geochem. 9 (5), 583–595. http://dx.doi.org/10.1016/0883-2927(94)90020-5. Strömberg, Bo, Banwart, Steven, 1999a. Weathering kinetics of waste rock from the Aitik copper mine, Sweden: scale dependent rate factors and pH controls in large column experiments. J. Contam. Hydrol. 39 (1–2), 59–89. http://dx.doi.org/10.1016/S0169-7722(99)00031-5. Strömberg, Bo, Banwart, Steven A., 1999b. Experimental study of acidityconsuming processes in mining waste rock: some influences of mineralogy and particle size. Appl. Geochem. 14 (1), 1–16. http://dx.doi.org/10.1016/ S0883-2927(98)00028-6. US Army Engineer Waterways Experiment Station (USAEWES), 1960. "The Unified Soil Classification System.". Technical Memorandum No. 3-357. Vicksburg, MS. April 1960. van Genuchten, M.T., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous media I. Anal. Solutions 40 (4), 8. Vanderborght, Jan, Vereecken, Harry, 2007. Review of dispersivities for transport modeling in soils. Vadose Zone J. 6 (1), 29. http://dx.doi.org/10.2136/ vzj2006.0096. Vanderborght, J., Timmerman, A., Feyen, J., 2000. Solute transport for steadystate and transient flow in soils with and without macropores. Soil Sci. Soc. Am. J. 64 (4), 1305. http://dx.doi.org/10.2136/sssaj2000.6441305x. Wierenga, P.J., 1977. Solute distribution profiles computed with steady-state and transient water movement models. Soil Sci. Soc. Am. J. 41 (6), 1050. http://dx.doi.org/10.2136/sssaj1977.03615995004100060006x. Yazdani, J., Barbour, L., Wilson, Ward, 2000. Soil Water Characteristic Curve for Mine Waste Rock containing Coarse Material. Proceedings of the 6th Environmental Engineering Specialty Conference of the Canadian Society of Civil Engineers. London, Ontario.