In situ leach mining (ISLM) is an advanced selective mining technique used for the exploitation of copper oxide deposits. ISLM involves the injection of a solution into a deposit, the dissolution and mobilization of target metals, and the recovery of this solution for processing. The vast supply of raffinate, a sulfuric acid by-product resulting from conventional copper ore processing, makes it possible to maintain continuous, high-frequency, applications during ISLM. Since solution injection and recovery in a copper oxide deposit is governed primarily by the hydraulic properties, the latter are the major factors affecting the ISLM process.

To date, the ISLM of copper oxide in a fractured crystalline ore deposit has enjoyed limited success, as compared with ISLM of uranium. The primary setback is attributed to the lack of attention given to the scale-dependence of hydraulic properties. To be able to develop a rational and sound hydrologic design, engineers must first identify the underlying hydraulic process at the scale of leaching in the deposit. This paper describes the application of a stochastic continuum approach for describing heterogeneous hydraulic properties at the Cyprus Casa Grande in situ leach mine, near Casa Grande, AZ. The use of simulated annealing (SA) in conjunction with deterministic flow modeling is presented as a means for overcoming scale dependence; ultimately, it extends the engineers ability to more objectively develop ISLM designs in a heterogeneous domain.

Scale Dependence

Representative Elemental Volume

Traditionally, the representative elemental volume (REV) is defined as that volume over which an ensemble of grains (or micropores) statistically describes the medium as a representative continuum. Specifically, this is the volume over which physical property measurements are statistically homogeneous and can be represented by a single value. By virtue of its definition, the REV concept can and is commonly extended to hydraulic properties, e.g. permeability, hydraulic conductivity, or transmissivity. Furthermore, the REV implicitly acknowledges the possibility of scale dependent properties. That is, there exists some characteristic volumes associated with at least a lower and an upper limit beyond which heterogeneity exists. In the ISLM design of uranium operations, the use of a single-valued (average) field hydraulic property proved satisfactory, because the scale of the well pattern (5-spot) is small relative to the scale of heterogeneity that exists in a typical roll front deposit.

Representative Elemental Domain

In figure 1, the original REV concept is generalized to include various scale dependent hydraulic heterogeneities encountered in a fractured copper oxide deposit. In this figure, the line around which measured hydraulic properties fluctuate represents, statistically, the scale-dependent mean. The domain bounded by a minimum and maximum REV, e.g. REVlto REVm is defined as a representative elemental domain (RED). Those deposit volumes associated with a RED can be expected to behave in a statistically homogenous manner, i.e. as an equivalent porous media. Conversely, departures from the hydraulic mean represent increasing statistical heterogeneity.

In a fractured copper oxide deposit, the scale dependence is attributed to the introduction of hierarchial hydraulic heterogeneities. For example, at some relatively small-scale the hydraulic properties may be ascribed to variations in capillary pore size distribution. At this scale, the RED1 is analogous to those volumes over which ISLM design occurred for uranium operations. At an intermediate- scale, the hydraulic behavior may be related predominantly to the introduction of microfractures. These two types of hydraulic heterogeneities were identified in a recent laboratory analysis. At the large-scale, typical hydraulic field measurements in fractured crystalline rock tend to be erratic and subject to the volume of rock influenced by the test. This behavior is attributed to test volumes intersecting large-scale fractures (joints, partings, and other macro-fractures), many of which are hydraulically inactive. While a REV may exist in such a fractured deposit, it would be much larger than the ISLM scale at which detailed hydraulic information is required and can be practically obtained. Hence, the typical scale of ISLM design in a crystalline deposit falls somewhere between the RED2 and RED3. These observations suggest that the critical problem is to extend the engineers ability to evaluate ISLM designs over scales not characteristic of a REV. One promising approach is to treat the hydraulic properties as a realization of a random process. With this in mind, the remaining analysis focuses on a stochastic treatment of hydraulic properties measured at the Cyprus Casa Grande ISLM facility.

Hydrologic Field Measurements

In Situ Leach Mining Facility

The Cyprus ISLM test site, near Casa Grande, AZ, is developed in a porphyry stock consisting of an intensely fractured granodiorite copper oxide cap overlying a sulfide base. The ore body is bounded on its lower side by the lakeshore fault, while its upper side is overlain by roughly 700 to 1,500 ft of cemented, unmineralized gravels (fig. 2). Chrysocolla and brochantite are the predominant copper minerals in the oxide zone.

The ISLM operation, located approximately 700 ft below the surface, consists of 58 wells drilled, in a series of 15 divergent fan patterns, along the 700 level drift (fig. 3). Wells penetrate the cemented gravel deposit, leached cap rock, and oxide ore zone, and terminate at the top of a sulfide zone. Each well is cemented in place in the interval between the top of the oxide zone and well collar, and perforated at one-foot intervals throughout the oxide ore. Pressure and flow sensors were installed on the wells for continuous monitoring of conditions during two years of leaching operations at the mine.

Transmissivity Estimates

At the Casa Grande mine, the flow of leach solution is assumed to be confined to macrofractures in the copper oxide ore. This conjecture is supported by post leach coring which showed evidence of solution mineral contact only in macroscale fractures. The confinement of leach solution to macrofractures is attributed largely to the unsaturated conditions occurring at the site.

Recognizing that the contact-between cemented gravels and copper oxide undulates in cross section, a two-dimensional approximation is constructed for modeling purposes. In this approximation, both wells and drifts are projected onto what is called the “hanging wall” cross section. Thus, this cross section shows the locations of the 58 wells and associated drifts, as they intercept a plane passing through the ore zone and along the top of the perforated interval in each well (fig. 3).

Since pressure and flow measurements reflect conditions occurring over a borehole interval, this investigation focuses on the hydraulic property called transmissivity. The transmissivity, as defined here, represents the product of the borehole length and effective fracture conductivity. Recognizing that the pressure diminishes exponentially with distance from a well, the sample support effectively represents a cylinder on the average 100 ft in length. Transmissivity provides a convenient mechanism for assessing the hydraulic behavior occurring normal to the orientation of injection and recovery wells.

Direct inversion of a numerical flow simulator was used to obtain local estimate of transmissivity, at the 58 well locations using the analytic element technique. From the recorded hydrologic data, an 11-data period was identified during which approximately steady-state pressure and the conditions existed. An appropriate data se was assembled including relevant head and flow specified boundary conditions to initiate the inversion process. Both well and drift locations were head specified and described analytically. In using this technique, a reference point, described in terms of head, was placed outside the influence of the well field, but coincident with an observation well where hydraulic head was zero.

Univariate Statistics

Since the transmissivity estimates were determined to be approximately log-normally distributed, all of the data was transformed by taking the logarithm. To examine the distribution of log-transmissivity, univariate statistics including the mean values and standard deviations, were determined for the hydraulic field. Table gives these summary statistics.

From the summary table, the transmissivity distribution is noted to be heterogeneous, spanning roughly 4 orders of magnitude. The distribution is for the most part symmetric, but has a minor positive skew. While summary statistics give information as to central tendencies, they do not provide a means to quantify the spatially dependent hydraulic process that may exist at the site. By computing the variogram, the spatial structure associated with the hydraulic properties can be investigated.

Structural Analysis

Omni-directional Variogram

Both the experimental and modeled omni-directional (90 degree tolerance) log-transmissivity variograms are shown in figure 4. For log-transmissivity, the spatial function displays a nugget effect, 0.3 ft 4/min², increasing with distance until it stabilizes with a limiting value (sill) of 0.8 ft 4/min² at a range of about 135 ft. The isotropic experimental variogram was fit to a theoretical variogram using a spherical model. The crossvalidation technique was used to test qualitatively the variogram fit.

The fact that it is possible to derive a variogram for the Casa Grande ISLM transmissivity data implies that the hydraulic process is spatially correlated. In such a case, the deposit volume over which the transmissivities are derived is not characteristic of a REV. In fact, the range may be considered as the lower bound to the RED3, and therefore a continuum approach cannot be applied to ISLM design at this site.

Directional Variogram

The classical variogram function models the spatial continuity of a random function. Change in the variogram model as a function of direction is referred to as geometric anisotropy. In this study, geometric anisotropy is reflected by a changing range with azimuth (constant variance). The directional dependence of log-transmissivity is interpreted as varying degrees of fracture connectivity throughout the well pattern; ultimately, it manifests itself as directional transmissivity. Figure 5 shows the corresponding range ellipse; the major axis occurring at 135 degrees (from the vertical) has a range of 175 ft, while the minor axis occurring at 45 degrees has a range of 70 ft (ratio = 0.4).

Estimation

In practice the total number of hydraulic property measurements available will be limited due to concerns such as time cost, and the environment. Thus, it becomes desirable to be able to estimate this phenomenon at unknown locations. Estimation of values at these unknown locations demands, however, that additional information or assumptions be made. While a deterministic model could be used, it requires that there is sufficient knowledge about the processes that drive the hydrologic response. At the field scale, however, the complexity of in situ processes force the engineer to admit that uncertainty exists between sample locations.

Kriging

Kriging, e.g. ordinary or indicator, is typically used to estimate field parameters, e.g. ore grades and transmissivity that show a definable spatial structure. Kriging, like other estimation procedures, e.g. polygonal, triangulation, and inverse distance squared, avoid large errors by being close on average. By being close on average, however, these techniques result in a tendency to smooth, or filter, extreme values. Since the extreme values exert the primary hydraulic control during ISLM in a fractured ore deposit, it is advantageous to use a stochastic approach which honors both the statistical structure (expressed through the variogram), as well as extreme values.

Another disadvantage with kriging (and other estimation procedures) is that it involves statistics which are deterministic. That is, repeated application to the same problem yields the same estimate. Thus, kriging does not have a mechanism to account for uncertainty in the transmissivity estimates. This is particularly important, since a project engineer must typically meet an objective set within technical, economic and environmental constraints. One approach that accounts for the spatial structure, extreme values, and provides for analysis of alternatives is a stochastic simulation technique called simulated annealing (SA).

Stochastic Simulation

Simulated Annealing

Generating alternate conditional stochastic images of field variables using SA is a rather new approach. In this approach, the initial transmissivity image is sequentially modified by swapping values in pairs of grid nodes not involving a conditioning datum, i.e. well data. A swap is accepted if the objective function (defined essentially by the square difference between the simulated and model variogram) is lowered. Since not all swaps which raise the objective function are rejected, this approach allows for getting out of local minima.

The ultimate success of the SA method depends on a slow “cooling” of the image controlled by a temperature function that decreases with time. The higher the control parameter, the greater probability that an unfavorable swap will be accepted. The criteria whereby a swap is accepted or rejected lies in the acceptance probability distribution given by the Boltzman distribution. When additional swaps do not lower the objective function, the simulation is complete.

The specification of how to lower a temperature is known as the “annealing schedule.” While there are mathematically based annealing schedules that guarantee convergence, they are too slow for practical application. The annealing schedule used in this study is based on an empirical approach outlined by Deutsch. Here the maximum number of perturbations at any one temperature is set equal to 100 times the number of nodes, while the acceptance target is set equal to 10 times the number of nodes. The image is considered to have cooled sufficiently, i.e. annealing process stops, once the convergence criterion, i.e. a low objective function equal to 0.001, is reached.

Simulated Transmissivity Field

In this section, SA is used to estimate spatial distributions of log-transmissivity from the Cyprus field measurements. The final transmissivity images are subject to the condition that the results match known values at their respective well location and have desired statistical properties imposed. Since SA is conditioned at well locations, the transmissivity values in each stochastic image reflect the full range observed in the field, e.g. from -6.5 to -2.5. From one realization to another, however, local variations are present between wells.

The required size of the simulation field depends on the type of random function and range of its heterogeneities. Providing that the random function is stationary, and the simulation field is large enough, the statistics of a realization should exactly equal the model statistics. This implies that accurate statistical reproduction is not based on discretization of the domain, but on the relative dimensions of the simulated field. A rule-of-thumb commonly used when conducting stochastic simulations is that the simulated field be roughly 10 times the variogram range.

In the present study, it was determined that the model statistics could be matched using a field of roughly 2½ times the variogram range. Figure 6 displays both the model and simulated variogram derived from the first realization. The excellent agreement between these two variograms is expected, since the simulation algorithm aims at drawing realizations that reflect those statistics modeled from the field data. It should be pointed out, however, that the results are only as good as the hydraulic information used as input to the stochastic simulator. Also, since the model statistics are inferred from sample statistics that are uncertain, i.e. limited field data, particularly at the tails of the cumulative density function, the exact reproduction of the model statistics may not be desirable.

In all of the simulations, the discretization of the field was roughly 1/25th of the variogram range; this resulted in 5 ft by 5 ft grid elements. Since the variogram range changed depending on modeled structure, the field size and total number of grid nodes also changed. For example, the smaller range associated with the isotropic structure required only 4500 grid nodes, while that for the anisotropic structure required 7600 grid nodes. Despite the use of these relatively modest number of nodes, the execution times for simulations were considerable. For example, a single isotropic simulation took about 115 min and 470 min using an HP-760 (small work station), and 486-33 MHz (personal computer), respectively. With the addition of nodes and/or anisotropic structure, the simulation time increases.

Figure 7 gives the contour reference (inverse distance squared), ordinary kriging and simulated REV imaging results. For comparative purposes, the well location and corresponding magnitude of transmissivities are superimposed onto the contoured reference diagram (fig. 7a). Both of these deterministic images indicate two regions of low transmissivity; one that transects the image from top to middle right, and one that protrudes from the bottom right corner inward. As expected, the kriging estimator resulted in a characteristically smooth and somewhat filtered distribution (fig. 7b). This smoothing effect is typically most pronounced in areas of sparser sampling. Despite its poor representation of extreme values, the kriged map remains useful for identifying and interpreting global heterogeneity at the ISLM design scale.

In contrast to the smoothly varying kriged map, the stochastic simulations (figs. 7c, 8, and 10) appear more ragged. Despite the fact that the stochastic image shown in figure 7c is conditioned to the field data, it bears no resemblance to any of the related figures (7a-b, and 8). The reason is that this image reflects a statistically homogeneous transmissivity distribution, or an REV. This image was produced by neglecting the scale-dependent structure reflected statistically by the variogram, i.e. range is set equal to zero resulting in a “pure nugget” effect.

Figure 8 gives five typical stochastic log-transmissivity realizations for the same period, but characterized by statistical anisotropy. As expected, the these images bear no resemblance to the simulated REV. The four higher transmissivity regions show in the reference figure (7a) are also represented here. This fact reinforces the ability of SA to reproduce global structure Locally, the degree of randomness that appears from one stochastic image to another represents the degree of uncertainty present.

In the previous example, the success associated with simulating the log-transmissivity field is attributed largely to the incorporation of anisotropic structure. In practice, however, the limited number of hydrologic measurements will likely preclude the identification of statistical anisotropy. Hence, the question remains as to what impact this may have on the spatial imaging of hydraulic features. To address this issue, and the question of time invariance, a second interval of hydrologic measurements were analyzed. Variogram analysis revealed an isotropic structure with the following characteristics: nugget – 0.3 ft4/min², range – 85 ft, and sill – 1.80 ft4/min². The contoured reference map is shown in figure 9.

From one interval to another, there is a ubiquitous change in log-tranmissivity. This emphasizes the temporal nature of scale-dependent ISLM hydraulic properties. While many factors may contribute to this type of behavior, this phenomenon is attributed largely to a change in well conditions, i.e. change in operational mode and injection/recovery pressures. The primary transmissivity features include three high and three low capacity regions. Two of the three high transmissivity regions occur adjacent to the drifts, while the third appears at the 1050 ft elevation extending roughly 20 ft inward. The three low transmissivity regions appear along a diagonal of negative slope at the 1050 ft, 975 ft, and 925 ft elevations.

Figure 10 gives six log-transmissivity realizations simulated based on the interval 2 statistics and those variogram characteristics outlined above. In comparing these isotropic simulations to the reference figure, the only globally consistent feature appearing in all of the images is the high transmissivity region adjacent to the lower drift. The remaining high/low capacity regions appear in three or four of these six images. These observations suggest that the identification and use of directional dependence in stochastic imaging is important for accurate imaging of the hydraulic process. This assertion is consistent with the results obtained from interval 1 simulations conducted using isotropic structure.

Identifying Uncertainty-Risk Assessment

The total number of realizations to be computed depends on the degree of uncertainty present in the analysis. While several realizations should be computed to identify uncertainty, it is tempting to retain only one for subsequent conditioning of a deterministic flow simulator. The selection of a single image based on subjective appreciation effectively disregards any uncertainty. A more objective approach is to quantify the degree of uncertainty by constructing probability distributions for any/all critical field locations.

In these probability plots, the mean represents the “expected log- transmissivity” that might be achieved on the average, if this same set of statistics were used. For this study, the standard deviation embodies the degree of risk involved in assigning a log-transmissivity to a given node. The more random and uncertain the hydraulic process is, the higher the standard deviation (or risk) will be by identifying locations where the risk is too great (as established by some cutoff probability), it is possible to target additional drilling and hydraulic test locations. Unfortunately, the imposition of various economic and environmental constraints renders this approach unlikely. Hence, improvements to the simulation technique must be considered to further constrain the nonuniqueness associated with stochastic imaging, while using an isotropic structure.

To further reduce this nonuniqueness (and thus uncertainty), additional constraints may be incorporated into the objective function. The ability to modify the objective function is one distinct advantage that SA has over other stochastic simulation approaches. For example, the spatial statistics associated with permeable zones identified from seismic or radiowave difference tomograms might be incorporated into the objective function.

ISLM Design Application

The ISLM design problem focuses on optimizing well conditions to produce a uniform distribution and recovery of leach solutions. In a heterogeneous deposit, the prediction of hydraulic variables at recovery and observation wells will be subject to uncertainties, particularly since the hydraulic properties represent the input parameter of greatest uncertainty when flow modeling. The following approach is suggested as a general guide to hydraulic ISLM design:

1. Assemble site specific hydraulic field property set. Hydraulic properties, e.g. hydraulic conductivity, or transmissivity, may be derived from either direct, e.g. pump or injectivity tests, or indirect, e.g. inversion based on pressure methods.
2. Simulate hydraulic property field. Using the simulated annealing technique, compute many alternate realizations. In areas of greatest uncertainty, additional drilling and/or hydraulic field testing may be considered.
3. Conduct flow simulations. Since a hydraulic property distribution (realization) is used to condition the deterministic flow simulator, a compatible element grid must be generated. If the hydraulic properties were nonlinearly transformed for geostatistical analysis, they should be back transformed prior to conditioning. Since well spacing is generally fixed, the engineer must decide as to which, and how many, wells will be used for injection. Following the assignment of a constant hydraulic condition at the injection well(s), a series of steady-state flow simulations could proceed. During each simulation, recovery well pressures would be iteratively updated until some objective criteria is matched, e.g. uniform flow rate. This process is then repeated using other hydraulic property realizations.
4. Risk Assessment. By constructing a pressure head probability density function for each of the recovery wells, the amount of uncertainty with regards to the desired flow race can be quantified, e.g. defined in terms of one standard deviation. Similar statistics could be calculated for other hydraulic concerns, such as traveltimes and localized flow path deviations. These hydrologic concerns are particularly important in assessing the operators ability to capture pregnant leach solution.

Heap Leach Design

In addition to ISLM mining, the stochastic continuum approach can be extended to address the general problem of heap leaching. While the objective of uniform leaching remains the same, the basic design strategy would focus on optimizing the spacing and discharge characteristics of the sprinklers. Instead of measuring discrete values of a hydraulic property, a nonlinear hydraulic property function would be determined at each location, due to its dependence on moisture content (or saturation). By combining basic hydrologic principles, the spacing and sprinkling rates can be determined for optimum exploitation of the neap.

Other Considerations

In using the stochastic ISLM design approach, the problem associated with a scale-dependent hydraulic process is overcome. This approach, however, does not account for coupled processes, i.e. flow-deformation-chemistry. In most cases, the coupled flow-deformation process can be avoided by maintaining an injection pressure below that required for dilation. Other processes such as flow- chemistry cannot be neglected.

Assuming that the leach solution has a low pH and is not saturated with respect to gangue cations, the dissolution and mobilization of metals will increase the effective porosity and thus permeability. Conversely, if the fluid chemistry is not effectively maintained, localized precipitation may occur. Recognizing the dynamic influence that ISLM exerts on hydraulic properties with time, it is recommended that the stochastic ISLM hydraulic design approach be conducted at intervals throughout the evolution of the leaching process.

Conclusions

Knowledge of the spatial variability of a hydraulic process is important for management and control of leach solutions in both a copper oxide deposit or heap. In many cases, the variogram can be used to evaluate whether heterogeneous conditions exist at the scale of an existing or proposed ISLM operation. At the Casa Grande ISLM underground test facility, near Casa Grande, AZ, the hydraulic process is spatially correlated. This implies that hydraulic process associated with ISLM is characteristically heterogeneous. Under these conditions, the hydraulic property distribution must be described using a stochastic approach.

Since the extreme values exert the primary hydraulic control during ISLM, it is advantageous to use a conditional stochastic simulation technique, such as SA. Contrary to kriging, SA, retains the statistical structure, full range of hydraulic values, and it allows for the derivation of many alternate stochastic realizations. The latter is important, since the uncertainty at locations between conditioned well points can be reviewed objectively.

By using the hydraulic property realizations as a conditioning datum for flow modeling, an engineer can evaluate ISLV designs regardless of deposit conditions. If the limited number of hydraulic measurements precludes the identification of statistical anisotropy, however, the uncertainty (risk) increases. To reduce the degree of uncertainty, it is suggested that spatial statistics derived from low cost geophysical information, e.g. permeable zones identified from seismic or radiowave difference tomograms, be incorporated into the SA objective function. Also, recognizing that ISLM exerts a dynamic influence on hydraulic properties with time it is recommended that the stochastic ISLM hydrologic design be reevaluated at intervals throughout the evolution of the leaching process.