1 USING ALTERATION HALOS TO DETERMINE THE MASS OF VOLATILES EXPELLED AND THE RATE OF EXPULSION IN A MAGMATIC PORPHYRY SYSTEM WITH APPLICATION TO THE P...

Author:
Amos Parrish

0 downloads 8 Views 833KB Size

A Thesis Presented to the Faculty of the Graduate School of Cornell University In Partial Fulfillment of the Requirements for the Degree of Master of Science

by Rachel Shannon January 2007

© 2007 Rachel Shannon

ABSTRACT

The rate of volatile expulsion, the duration of venting, the total volume of volatiles expelled, and the size of intrusion required to supply them are key parameters in the formation of any magmatic-hydrothermal ore deposit that can be estimated from the morphology of the hydrothermally altered rock. The radius of pervasively altered rock around the fluid source and the taper of alteration halo with distance from the pervasively altered zone constrain the duration and rate of volatile expulsion. The volume of altered rock records the total mass of volatiles expelled. An analysis of the Pittsmont Dome at the porphyry copper system at Butte, Montana provides an example. The moles of hydrogen ion consumed during alteration are calculated from the mineralogy of unaltered and altered rock. Published composition estimates are used to determine the concentration of reactable hydrogen ion in the magmatic fluid. Diffusion from the vein to a reaction front at the edge of the halo surrounding it constrains the rate of hydrogen ion loss from the vein. Semi-analytic and finite difference simulations of alteration halo formation show how the radius of the pervasively altered zone and the steepness of the taper beyond that zone depend on the rate of volatile expulsion; the faster the fluid velocity, the less steep the taper. Data from the Pittsmont Dome suggest 23 to 30 billion tons of magmatic fluid was expelled over a period of less than 20 years. Assuming 5 wt% magmatic water in the porphyry intrusion, a spherical intrusion ~7 km in diameter is needed to supply the volatiles for just this part of the Butte mineralization system.

BIOGRAPHICAL SKETCH

Rachel Shannon grew up in Sheridan, Wyoming, where she gained an appreciation for the outdoors from an early age. After high school, she moved from New Mexico to Arizona to Colorado, working a variety of odd jobs. When she entered college in 1999 geology seemed an obvious choice of major, and she received her Bachelor's from the University of Colorado at Boulder in 2002. She worked for two years in the petroleum industry before entering the graduate program at Cornell to pursue her Master's. Upon completion of this degree, she is planning to move to Washington, D.C. to work for a hydrology and environmental consulting firm

iii

ACKNOWLEDGEMENTS

I would like to thank Steve Czehura, the chief geologist at Butte, who gave me access to the mine and allowed me to collect samples. His extensive knowledge of the porphyry system was invaluable. All the data presented here was collected Deep Diamond Drill core stored at the Montana Museum of Mines, so I am grateful that it has been so carefully stored and indexed. Thanks to everyone for their assistance during my time in Butte, especially Camela Carstarphen of the Montana State Geological Survey and Chris Gammonds and Colleen Eliot of Montana Tech. John Dilles of Oregon State University allowed me to join his field trip to Butte and helped arrange my field work there. He is one of the best field geologists I know, and his understanding of the geology at Butte was integral in my understanding of it. Mark Reed of Oregon State University provided a great deal of guidance, particularly on the geochemistry and mineralogy at Butte, and I would like to thank him for many useful discussions. I am grateful to the Penrose organization for sponsoring a meeting on the Butte System. Without this meeting my involvement in this kind of research would have been impossible. I would like to thank all the scientists I have met at Cornell, whose continuing enthusiasm has been inspiring. First and foremost, I am very grateful to my major advisor, Larry Cathles, for innumerable meetings, discussions, and careful edits of this work. This thesis would not have been written without his scientific knowledge, patience and support. Thanks to my other committee member, Andrew Pershing, for teaching me how to code in Matlab and for careful edits of this paper.

iv

Finally, I appreciate all the encouragement I have received from my parents, John and Hannah Shannon, my brother Michael Shannon, and my sister, Renee Goyeneche. I am grateful for the support and humor of my friends and fellow graduate students at Cornell, especially Danica Wyatt, Patrick Carr, Tiffany Tchakirides, Ursula Smith, Katie Tamulonis, and Louise McGarry.

v

TABLE OF CONTENTS

Introduction .................................................................................................................... 1 Theory............................................................................................................................. 5 Semi-analytic analysis ................................................................................................ 5 Finite difference confirmation .................................................................................... 9 Extending finite difference simulations to three dimensions and including pervasive alteration ................................................................................................................... 11 Parameters for the Butte, Montana Porphyry Copper System ..................................... 16 Mineralogy and pre-main stage veining ................................................................... 16 Vein spacing ............................................................................................................. 18 Maximum halo width zmax, taper, and radius of pervasive alteration zone, p........... 24 Consumption of reactable hydrogen, G .................................................................... 26 Results .......................................................................................................................... 33 Discussion..................................................................................................................... 37 The base case interpretation...................................................................................... 37 Parametric variations ................................................................................................ 38 Total intrusion volume required by the Butte porphyry system ............................... 39 Implications .............................................................................................................. 41 Summary and Conclusions ........................................................................................... 43 Appendix: Matlab code for finite difference simulations............................................. 44 References .................................................................................................................... 53

vi

LIST OF FIGURES

Figure 1: An idealized porphyry copper system showing concentric zones of hydrothermal alteration around an intrusive stock (after Guilbert and Park 1986). Alteration occurs as halos around veins which are evenly distributed throughout the dome. When veins are close together, halos overlap and alteration appears pervasive................................................................................................................. 3 Figure 2: The taper of a vein halo at one stage of development for a single vein unaffected by neighboring veins. Hydrothermal fluid from a magmatic source with initial reactable hydrogen concentration C ∑ H + enters the vein at the left (x = 0

0) and diffuses through the alteration halo to react with unaltered host rock in a thin reaction zone. The reactable hydrogen in the vein fluid decreases with distance along the vein as illustrated at the top. The vein taper is approximately

z max x max

linear, and equals . Zmax occurs at the source where concentration of reactable hydrogen is highest. Both xmax and zmax increase as venting continues. After venting is complete, the vein and halo extend from the source to the radius of the porphyry shell; i.e. xmax = R. Three cells of the finite difference simulation (greatly enlarged) are indicated at the left and labeled a, b, c. The first stages of alteration in these cells are shown in figure 3. ....................................................... 4 Figure 3: First stages of alteration in the first three finite difference cells located as shown in figure 2. The white area at the bottom of each cell represents the vein. The grey areas represent the distance perpendicular to the vein which the halo moves during a timestep. Three timesteps are shown. The stippled pattern at the top represents the unaltered host rock. At time = 0, the source fluid with composition C∑ H + is introduced into the first cell (a). During the timestep, 0

reactable hydrogen is allowed to diffuse out of the vein and into the host monzonite. Over one timestep, the halo grows by ∆z according to ∆z =

ρbG j ∑ H + ∆t

,

where j ∑ H + is given by equation 1. The partially spent fluid is then moved into cell (b) and source fluid moved into cell (a). The halo of cell (b) advances one step and the halo of cell (a) advances a second step. As this process is repeated the vein halos in the cells grow and cells further along the vein develop halos... 10 vii

Figure 4: Linear single vein finite difference simulation of the formation of an alteration halo around one vein for EDM (top) and sericitic (bottom) alteration, showing the relationship between fluid flux and halo taper. The solid black lines show how the halos develop over time. The topmost line in each graph shows the shape of the alteration halo after the total time of halo formation, 4.5 years for EDM halos and 2.7 years for sericitic halos. Earlier profiles are given every 0.45, and 0.27 years, respectively. Lower fluxes yield a steeper taper, but all fluxes result in a linear decrease in alteration halo width with distance along the vein. Parameters are: C ∑ H + = 3 moles/kg fluid, f = 0.6, φ = 1%, τ = 5, G = 4.68 0

moles reactable hydrogen/kg rock and temperature = 600 C for EDM alteration, G = 3.55 moles/kg rock and temperature = 400 C for sericitic alteration. ........... 12 Figure 5: The transformation from the linear model with flow at a constant rate through a single vein to a spherical model with radial flow out from the center of the porphyry. In the single vein model, each cell is the same width, contains the same volume of vein fluid, and is evenly spaced along the length of the vein (a). The halo tapers linearly from x = 0 to x = xmax (b). To convert to the spherical model, each two-dimensional cell becomes a three dimensional shell nested around the center of the porphyry. Each shell contains the same volume of pore fluid; therefore, the outer edge of the shells get closer together with distance from the source (c). In the spherical transformation, the halos merge in the zone of pervasive alteration and there is a steep taper beyond this zone The average maximum halo width determined as described in figure 5 is set as a limit in the spherical model; once the halo width has reached this limit, reaction no longer occurs in that shell (d). ......................................................................................... 15 Figure 6: Cross section modified from Reed et al 1997, showing the porphyry shells (defined by the magnetite vein zones) and the location of drill holes. The east block of the Continental Fault has been restored to approximate pre-faulted location by adding 1300 meters in depth. The fluid source is the solid black circle located roughly at the center of a circle defined by the Pittsmont dome. All drill core measurements in figures 6 and 8 are plotted as a function of radial distance from this location. The minimum radius of pervasive alteration, p, is shown as a solid black curve around the source. In calculating distance from the source, 1300 meters depth is added to measurements taken from drill hole 10 in order to approximate pre-faulted conditions. ....................................................... 17 Figure 7: Pre-main stage vein density in the porphyry shell measured from drill core. Observed pervasive alteration where no distinct halos can be observed in the core is shaded in grey. The grey circles show measured vein density from drill core vs viii

distance from the center of the porphyry shell The dashed line shows a best-fit line which represents the average vein density as a function of distance. Because the average density changes very little with distance, here it is considered to be constant at 22 veins/ meter. The maximum halo width is constant at one-half the 1 1 distance between two veins, ⋅ = 0.023 m. .............................. 19 2 22 veins / meter Figure 8: Measured and calculated alteration halo widths versus distance from the fluid source (figure 4) for expelled volatile fluid mass of 30 billion tons for EDM alteration and 23 billion tons for sericitic alteration. Circles represent halo widths measured from core. Maximum halo width determined from vein spacing is shown as horizontal dashed lines; the radius of pervasive alteration, p, is the distance to which the calculated halo width equals maximum halo width. Each curve represents a different duration of venting: in the EDM case, the curves shown are for venting duration of 6.5, 13.0, and 26.0 years. In the sericitic halo case, curves shown are for venting duration of 3.1, 6.2, and 12.4 years. The results which best match measured alteration halos are shown in bold. In this case, duration of venting is 13.0 years for EDM alteration and 6.2 years for sericitic alteration. ................................................................................................ 25 Figure 9: Results of the spherical finite difference simulation of formation of alteration halos around veins for EDM (top) and sericitic (bottom) alteration, in which the fraction of reactable hydrogen concentration which drives diffusion, f was varied. Maximum halo width determined from vein spacing is shown as horizontal lines with dashed extensions. The solid black lines show how the halos develop over time. The topmost line in each graph shows the shape of the alteration halo after the total time of volatile venting, 13.0 years for EDM and 6.2 years for sericitic alteration. Earlier curves are at one-year intervals. Parameters are: C ∑ H + = 3 moles/kg fluid, f = 0.6, φ = 1%, τ = 5, G = 4.68 moles reactable 0

hydrogen/kg rock and temperature = 600 C for EDM alteration, G = 3.55 moles/kg rock and temperature = 400 C for sericitic alteration. .......................... 40

ix

LIST OF TABLES

Table 1 : Definition of terms and values used in the model ......................................... 20 Table 2: Mineral composition of fresh Butte Quartz Monzonite and the alteration halo assemblages .......................................................................................................... 27 Table 3: Proxy mineral assemblages for fresh Butte Quartz Monzonite and alteration halo assemblages .................................................................................................. 28 Table 4: Calculated difference between alteration halo assemblages and Butte Quartz Monzonite (negative sign indicates minerals destroyed during alteration).......... 30 Table 5: Stoichiometric matrices for alteration (based on mineral/fluid reactions listed in Geiger et al. 2002) ............................................................................................ 31 Table 6: Basis species change during alteration compared to fresh Butte Quartz Monzonite (negative sign indicates basis species removed during alteration)..... 32 Table 7: Calculation of hydrogen yield from CHILLER calculations at 400 deg Celsius. From Mark Reed, personal communication .......................................... 34 Table 8: Results for the Butte, Montana porphyry copper system .............................. 36

x

Introduction

More than half the world's supplies of copper and of molybdenum have been mined from porphyry copper systems (Singer 1995). While porphyry systems are arguably the single most important source of these minerals, they are becoming increasingly difficult to find. Most exposed porphyry ores have likely already been discovered (Richards 2003). As efforts focus on deeper deposits, explorationists must begin to rely on geophysical methods constrained by less direct observations (Richards 2003). Exploration will be more effective if the large-scale geologic features associated with porphyry systems are well constrained by knowledge of the physical and chemical processes that formed them.

Fortunately, many years of description

and study have provided an extensive knowledge base, and the source, evolution, and final morphology of porphyries are fairly well understood (see reviews by Richards (2003), Henley and Berger(2000)). Porphyry copper systems are dome shaped volumes of hydrothermally altered and mineralized rock which form around an intrusive stock. The porphyry domes show concentric zones of alteration: potassic at the center, surrounded by phyllic, argillic, and propylitic from the center outward (Guilbert and Park 1986). The mineralization is contained within a large number of veins which are usually thin (~millimeters wide) and are evenly distributed at any given location, (figure 1). Alteration often occurs as distinct halos surrounding these veins (figure 2).

Near the

center of the system, vein halos are so wide that the alteration appears pervasive. Porphyry copper systems are thought to have formed when a hot water-rich acidic fluid is expelled from a small apophysis at the top of a molten intrusion

1

(Richards 2003). As the magmatic fluid moves outward through fractures and reacts with rock, minerals are deposited in the fractures and alteration halos develop around them. The source of the fluid has been debated in the literature. Most recent authors, however, conclude that the fluid which supplies early-stage veins consists of volatiles that are original components of the magma (see reviews by Richards 2003 and Henley and Berger 2000). The duration of venting and the size of the intrusion required to supply the needed volume of magmatic fluid remain open questions, although some estimates have been made (e.g., see Geiger et al., 2002 for estimates for the Butte porphyry system). Here we address these questions with a semi-analytic model and finite difference simulations. Both assume the formation of an alteration halo is controlled by diffusion. Beyond the radius of pervasive alteration, halo width decreases with distance along a vein because the arrival of acidic fluid is delayed in time. The taper of the halo with distance from the edge of the pervasively altered zone is controlled by the mass flux of the fluid in the fractures. A slower mass flux results in a steeper taper, and a faster mass flux results in a shallower taper. The duration of venting is measured by the radius of the pervasively altered zone and the maximum halo width (figures 1 and 2). The semi-analytic model illustrates this concept, but finite difference simulations are needed to apply them to an actual porphyry system. When this is done it can be seen that only a few direct observations are needed to estimate fluid volume (and therefore mineral potential), intrusion size, and duration of volatile venting.

2

Figure 1: An idealized porphyry copper system showing concentric zones of hydrothermal alteration around an intrusive stock (after Guilbert and Park 1986). Alteration occurs as halos around veins which are evenly distributed throughout the dome. When veins are close together, halos overlap and alteration appears pervasive.

3

Figure 2: The taper of a vein halo at one stage of development for a single vein unaffected by neighboring veins. Hydrothermal fluid from a magmatic source with initial reactable hydrogen concentration C ∑ H + enters the vein at the left (x = 0) and 0

diffuses through the alteration halo to react with unaltered host rock in a thin reaction zone. The reactable hydrogen in the vein fluid decreases with distance along the vein

z max x max

as illustrated at the top. The vein taper is approximately linear, and equals . Zmax occurs at the source where concentration of reactable hydrogen is highest. Both xmax and zmax increase as venting continues. After venting is complete, the vein and halo extend from the source to the radius of the porphyry shell; i.e. xmax = R. Three cells of the finite difference simulation (greatly enlarged) are indicated at the left and labeled a, b, c. The first stages of alteration in these cells are shown in figure 3.

4

Theory

Semi-analytic analysis

The semi-analytic model assumes that the rate of formation of an alteration halo is governed by one dimensional diffusion of reactable hydrogen ion. The halo grows as the square root of time as hydrogen ion diffuses from the fracture to the unaltered rock at the edge of the halo where the pore fluid, including pH, is in equilibrium with the host rock. Within the halo, fluid does not react with the rock (figure 2). Reaction occurs only in a narrow zone very near the edge of the halo (e.g. Levenspiel, 1967, Cathles and Apps, 1975, Cathles 1979). At any time t, the reactable hydrogen ion concentration gradient that controls hydrogen transport from the vein to the reaction zone is the concentration in the vein,

C∑ H + , times the fraction of the total concentration which drives diffusion, f, divided by the halo width, z. At any time t, the flux of reactable hydrogen ion, j∑ H + , is:

j∑ H + = −

Dφ ρf

τ

⋅

f C ∑ H+ z

,

(1)

where D is the diffusion constant of reactable hydrogen ion, φ is the vein halo porosity, τ is the vein halo pore tortuosity, ρ f is the density of the fluid filling the pores, and z is the width of the halo. Over the time t required for the halo to grow to width z, the cumulative flux of hydrogen ion from the vein, J is:

5

J = 2⋅

Dφ ρf

τ

⋅

f C∑ H + z

⋅t .

(2)

Here average flux is estimated by multiplying the flux at time t by two to account for the fact that diffusion is faster at early times when z is small. The average rate of diffusion between time 0 and time t is about twice the diffusion rate at time t. Vein halo width at time t is calculated by setting the cumulative loss of reactable hydrogen ion from the vein equal to that in the vein halo plus that consumed by reaction in the rock. If G is defined to be the moles of hydrogen ion per kilogram rock consumed in converting fresh host rock to halo and ρ r is defined as the density of the host rock, the hydrogen ion consumed is ρ r G z . The reactable 1 hydrogen in the halo is z φ ρ f C ∑ H + . Setting the flux from the vein over time t 2 equal to the sum of these terms yields:

2⋅

D φ ρ f f C∑ H + t zτ

= ρr G z +

1 z φ ρ f C∑ H + . 2

(3)

Rearranging for the width of the halo as a function of time yields the familiar diffusion equation modified to include reaction at the halo edge:

z (t ) = 2

f D t /τ 2 G ρr 1+ φ ρ f C + ∑H

6

(4)

Consider a simplified case of flow through a single vein where the source is at x = 0 and no neighboring veins interfere with halo growth. In this case, the duration of volatile venting is recorded by the maximum halo width, zmax in the porphyry system (figure 2). From eq 4, the duration of venting is:

t venting =

2 ρr G 2 ⋅ z max 1+ φ ρf C + 4D f ∑ H0

τ

.

(5)

Total volume of magmatic volatiles expelled through one meter of one vein can be determined in this simplified case. Maximum halo width occurs at x = 0. Halo width is zero at x = xmax.. We will show with the finite difference simulation that the width of a halo in this simple case decreases linearly with distance from the magmatic fluid source along the vein (figure 2). Thus the moles of hydrogen ion reacted per unit length of vein between x = 0 and x = xmax is ρ r z max G x max . This must equal the moles of reactable hydrogen delivered over the duration of venting, Qv C ∑ H + t venting , so 0

x max =

QV C ∑ H + t max 0

ρ r z max G

.

(6)

The halo taper is the maximum halo width divided by the maximum length:

taper =

2 z max z max ρr G = x max QV C ∑ H + t venting 0

7

(7)

Substituting for tmax from equation 4 into equation 7 assuming that

2 ρr G >> 1 φ ρ f C∑ H + 0

yields the fluid flux through one meter of one vein:

QV = 2 ⋅

D f φ ρf

τ taper

.

(8)

The total volatile flux in the system can be roughly estimated by multiplying the flux through one meter of one vein by the area of a representative surface surrounding the source and the vein density (meter length vein per m2 rock), as shown in figures 1 and 7:

QT = Qv A σ .

(9)

Total fluid mass expelled from the porphyry is total flux times the time duration of venting:

M f = QT t .

(10)

The total intrusion mass is the volatile mass divided by the mass fraction of volatiles in the intrusion: MI =

Mf fV

(11).

The above discussion shows that we can determine the duration of venting from the width of the halos at the source (equation 5), the rate of venting from the vein taper (equation 8), and the mass of the intrusion responsible for the mineralization from the density of veins on a surface enclosing the porphyry fluid source (equations 9-11).

Unfortunately, the situation in an actual porphyry is not

8

this simple because alteration is pervasive at the center of the system where the vein halos have merged, and the flow of fluid away from the source is roughly spherical which means as the fluid moves away from the source it must move through an increasingly larger number of fractures. A finite difference model of halo formation can accommodate these complications as well as confirm critical aspects of the semianalytic model just presented.

Finite difference confirmation

The semi-analytic model described above contains all the concepts needed to analyze a porphyry system. To verify it, especially the assumption of a linear halo taper, and to extend it to account for actual porphyry geometry, a finite difference simulation was built which describes the spacial and temporal evolution of a halo around one vein. The numerical simulation was written in Matlab. The code is attached in the appendix. The operation of the finite difference model is illustrated in figure 3. Away from the vein, the alteration halo grows by diffusion as a function of the square root of time according to equation 4. Distance along the vein is divided into cells of length l, where l equals the distance the fluid travels in one timestep. To avoid infinite diffusion rates, each cell is given a very thin halo at the start of the simulation. At time 0 fluid with initial reactable hydrogen concentration C ∑ H + is 0

introduced into the first cell. Diffusion is calculated perpendicular to the vein for the

9

Figure 3: First stages of alteration in the first three finite difference cells located as shown in figure 2. The white area at the bottom of each cell represents the vein. The grey areas represent the distance perpendicular to the vein which the halo moves during a timestep. Three timesteps are shown. The stippled pattern at the top represents the unaltered host rock. At time = 0, the source fluid with composition C∑ H + is introduced into the first cell (a). During the timestep, reactable hydrogen is 0

allowed to diffuse out of the vein and into the host monzonite. Over one timestep, ρbG the halo grows by ∆z according to ∆z = , where j ∑ H + is given by equation 1. j ∑ H + ∆t The partially spent fluid is then moved into cell (b) and source fluid moved into cell (a). The halo of cell (b) advances one step and the halo of cell (a) advances a second step. As this process is repeated the vein halos in the cells grow and cells further along the vein develop halos.

10

specified timestep assuming a linear gradient to zero across the width of the halo. The reactable hydrogen diffused from the vein over the timestep is computed using equation 1. At the end the timestep, the halo width in the first cell is increased in accordance with the amount of reactable hydrogen diffused out of the vein over the timestep. Reactable hydrogen concentration in the vein is decreased by this same amount, and the fluid is moved to the next downstream cell. Vein fluid in the first cell is replaced with fluid with the initial reactable hydrogen concentration of the source, C∑ H + . At each subsequent timestep, diffusion is calculated using the 0

reactable hydrogen concentration in each cell, the reactable hydrogen concentration is reduced by the amount diffused across the halo, and the fluid is moved to the next downstream cell. Reactable hydrogen concentration in the first cell remains constant at C∑ H + throughout the simulation. This process continues until the maximum 0

duration of venting, tventing, is reached. This simulation shows that the steepness of the taper depends on the fluid flux as predicted by equation 7. The taper is always linear from the source to the end of the vein, verifying a key assumption made in the semi-analytical solution (figure 4).

Extending finite difference simulations to three dimensions and including pervasive alteration

In an actual porphyry system flow is radial from the source and veins near the source have halos so wide that they overlap. These complexities can be accommodated by transforming the length axis of the single vein model and stopping reaction when the halos merge. The vein in each cell is now considered part of a hemispheric shell containing many veins; the shells are nested around the fluid source 11

Figure 4: Linear single vein finite difference simulation of the formation of an alteration halo around one vein for EDM (top) and sericitic (bottom) alteration, showing the relationship between fluid flux and halo taper. The solid black lines show how the halos develop over time. The topmost line in each graph shows the shape of the alteration halo after the total time of halo formation, 4.5 years for EDM halos and 2.7 years for sericitic halos. Earlier profiles are given every 0.45, and 0.27 years, respectively. Lower fluxes yield a steeper taper, but all fluxes result in a linear decrease in alteration halo width with distance along the vein. Parameters are: C ∑ H + 0

= 3 moles/kg fluid, f = 0.6, φ = 1%, τ = 5, G = 4.68 moles reactable hydrogen/kg rock and temperature = 600 C for EDM alteration, G = 3.55 moles/kg rock and temperature = 400 C for sericitic alteration.

12

(figure 5). The volumes of the veins in all shells are the same. At the start of the simulation, fluid enters the veins inside a hemispheric volume extending from r = 0 to r0 (across the first cell). At the end of the first timestep the reactable hydrogen ion concentration in these veins is reduced by the diffusive losses across the halo, and the vein fluid is moved to the veins in the next concentric shell.

The calculation is the

same as in the single vein case just described. At each timestep, the same volume of magmatic fluid is moved from shell to shell just as it was moved from cell to cell in the single vein model. If vein spacing is constant in the porphyry and all veins contain equal volumes of fluid per unit area of vein, the radius of the outer boundary of each shell is equal to ro 3 n where is n is the number of the shell (n = 1 represents the innermost sphere surrounding the source, n = 2 represents the first shell, etc). The maximum halo width possible before vein halos merge is one half the average distance between fractures (figures 7 and 8). Once the maximum halo width is reached, the rock between two veins can no longer neutralize acid. None of the rock in the shells closer to the source can react with hydrogen either, so the concentration of reactable hydrogen becomes C ∑ H + , and this fluid is moved to the next shell out. 0

Using this procedure, the model can simulate the expansion of the pervasive alteration zone with time.

The radius of the pervasively altered zone, p, is marked by the first occurrence of a calculated halo width less than the maximum halo width defined by the vein spacing. From this radius outward the halo width narrows with distance from the source. The halo taper from this point to the radius at which halo width goes to

13

zero is predicted by the model and constrained from drill core measurements. The total mass of magmatic volatiles vented is constrained by the total volume of halo alteration, and the rate of venting is constrained by the halo taper just as in the single vein case. The taper constraint applies from the radius of the pervasive zone to the edge of the porphyry shell. Thus the spherical model is exactly the same as the single vein model except that the distance scale is distorted and a zone of pervasive alteration precedes the zone of halo taper (figure 5).

14

Figure 5: The transformation from the linear model with flow at a constant rate through a single vein to a spherical model with radial flow out from the center of the porphyry. In the single vein model, each cell is the same width, contains the same volume of vein fluid, and is evenly spaced along the length of the vein (a). The halo tapers linearly from x = 0 to x = xmax (b). To convert to the spherical model, each two-dimensional cell becomes a three dimensional shell nested around the center of the porphyry. Each shell contains the same volume of pore fluid; therefore, the outer edge of the shells get closer together with distance from the source (c). In the spherical transformation, the halos merge in the zone of pervasive alteration and there is a steep taper beyond this zone The average maximum halo width determined as described in figure 5 is set as a limit in the spherical model; once the halo width has reached this limit, reaction no longer occurs in that shell (d).

15

Parameters for the Butte, Montana Porphyry Copper System

Mineralogy and pre-main stage veining

The host rock of the Butte, Montana porphyry copper system is the 76 Ma Butte Quartz Monzonite, part of the Boulder Batholith (Geiger et al., 2002). The monzonite is coarse-grained and is mostly composed of plagioclase, quartz, Kfeldspar, biotite, hornblende, magnetite, and titanite (Brimhall, 1977). The deposit developed in three stages called the pre-main stage (66-64 Ma), the main stage (64-62 Ma), and the post-main stage (60 Ma). Each stage shows different mineralogy. Main stage and Post-main stage mineralization has been well described in the literature (Brimhall, 1977; Meyer et al., 1968; Sales and Meyer, 1949; among others) and will not be discussed here. Pre-main stage mineralization occurs in two porphyry domes, each of which shows zoned mineralization which is depicted in figure 6 by the magnetite vein zone. Pre-main stage veins with halos can be divided into types based on vein and halo mineralogy. The older veinlets are usually less than one centimeter wide and are dominated by quartz. They are surrounded by EDM alteration halos (Meyer, 1965 and Brimhall, 1977), which contain most or all of the following minerals: muscovite, andalusite, quartz, alkali feldspar, calcite, pyrite, and chalcopyrite (Brimhall, 1977). Cross cutting the EDM veins are quartz/pyrite veinlets which show sericitic alteration halos. These veinlets are also millimeters to centimeters wide, with halos typically only a few centimeters wide. In these halos, the feldspar, biotite, and hornblende from the Butte Quartz Monzonite have been altered to varying degrees to sericite. 16

Figure 6: Cross section modified from Reed et al 1997, showing the porphyry shells (defined by the magnetite vein zones) and the location of drill holes. The east block of the Continental Fault has been restored to approximate pre-faulted location by adding 1300 meters in depth. The fluid source is the solid black circle located roughly at the center of a circle defined by the Pittsmont dome. All drill core measurements in figures 6 and 8 are plotted as a function of radial distance from this location. The minimum radius of pervasive alteration, p, is shown as a solid black curve around the source. In calculating distance from the source, 1300 meters depth is added to measurements taken from drill hole 10 in order to approximate pre-faulted conditions.

17

Table 1 lists the parameters required to characterize a porphyry system according to the method described here and shows the values used for the Butte, Montana porphyry copper system. These parameters were approximated or taken directly from previous publications as described in the table, except for the maximum halo width, halo taper, vein spacing, and consumption of reactable hydrogen. These parameters are described below. Vein spacing

Figure 6 shows an east-west cross section through the pre-main stage alteration at Butte from Reed (2005, personal communication). We have displaced the easternmost block of the Continental fault downwards by ~1300 meters to approximate pre-faulted geometry (Dilles, 2006 personal communication). The fluid source is roughly the center a circle defined by the reconstructed Pittsmont dome. Figure 7 shows measured number of veinlets per meter length of core as a function of distance from the fluid source. We assume that the vein density is equal in all directions. With these assumptions, the number of veins per meter length core equals the area of veins per cubic meter of rock, and the distance between veins in three dimensions is the same as the distance between veins in the core. The vein spacing is variable, but on average it does not increase much with radius from the center of the porphyry; therefore, in the model presented here spacing is considered to be constant at ~0.023 meters. This value is plotted as a dashed line in figures 7 -9.

18

Figure 7: Pre-main stage vein density in the porphyry shell measured from drill core. Observed pervasive alteration where no distinct halos can be observed in the core is shaded in grey. The grey circles show measured vein density from drill core vs distance from the center of the porphyry shell The dashed line shows a best-fit line which represents the average vein density as a function of distance. Because the average density changes very little with distance, here it is considered to be constant at 22 veins/ meter. The maximum halo width is constant at one-half the distance 1 1 between two veins, ⋅ = 0.023 m. 2 22 veins / meter

19

Table 1 : Definition of terms and values used in the model Symbol

Description

Value

Units

Source

φ

host rock (monzonite) porosity

0.01

[no units]

estimated for monzonite

τ

halo rock tortuosity

5

[no units]

estimated for monzonite

ρf

density of hydrothermal source fluid

400 (EDM)

[kg/m3 fluid]

estimated for the temperatures at halo formation, based on the density of water

ρI

density of source intrusion

2600

[kg/m3 intrusion]

estimated for felsic magma

ρr

density of host rock

2712

[kg/m3 rock]

Geiger et al 2002

σ

vein density

22

[veins/m length core] =

measured from drill core

600 (sericitic)

[veins/m3 rock] A

C∑ H + 0

in the 2D linear model, the area of a representative surface around the fluid source used to estimate fluid volume

10.2

[km2]

R 4π ( ) 2 2

concentration of reactable hydrogen in the magmatic source fluid

3.0

[moles critical species/kg magmatic fluid]

table 7

20

Table 1 (Continued)

C∑ H +

concentration of reactable hydrogen in the vein as a function of distance along the vein (x) and elapsed time (t)

variable

[moles critical species/kg fluid]

calculated in the finite difference simulations

D

diffusivity of critical species in water

1.19 (EDM)

[m2/year]

sericitic: Geiger et al 2002

x ,t

0.5 (sericitic)

EDM: calculated using Arrhenius equation f

fraction of total hydrogen concentration in magmatic source fluid which drives diffusion

0.6

[no units]

estimated as described in this paper

dt

timestep- in the finite element simulations, the amount of time fracture fluid diffuses into matrix rock between advection steps

variable:

[years]

optimized for specific simulations

fV

mass fraction of volatiles (fluids) in intrusion

0.05

[no units]

estimated for felsic magma

G

“Consumption”: moles of reactable hydrogen which must be added to host rock to create observed alteration, per kg of rock

4.68 (EDM)

[moles critical species/kg rock]

calculated as described in this paper

2 ⋅ 10 −4 to 5 ⋅ 10 −4

3.55(sericitic)

21

Table 1 (Continued) Mf

total mass of magmatic source fluid

23 to 30

[billion tons]

calculated as described in this paper

MI

total mass of intrusion

3.8 ⋅ 1011 − 6.2 ⋅ 1011

tons

calculated as described in this paper

Qp

in the 3D finite difference simulation, the initial fluid flux at time = 0

80 -140

[tons fluid/meter length vein/yr]

calculated as described in this paper

QT

flux of magmatic source fluid through entire porphyry system, estimated from 2D single vein model

8.7 ⋅ 10 8 − 2.0 ⋅ 10 9

[tons fluid/yr]

calculated as described in this paper

QV

in the 2D single vein model, the flux of magmatic source fluid through one vein

11 - 15

[tons fluid/m length vein/yr]

calculated from semianalytic model or single-vein finite difference simulation

p

radius of zone of pervasive alteration

1250 - 1400

[ m]

estimated from drill core data

R

radius of the porphyry shell

1800 - 2000

[m]

estimated from drill core data

r0

in the 3D spherical model, the starting radius

variable: 200 -500

[m]

calculated as described in this paper

T

temperature at halo formation

600 (EDM)

[degrees Celsius]

EDM: Brimhall (1977)

400 (sericitic)

sericitic: Geiger et al (2002)

22

Table 1 (Continued) tventing

total duration of magmatic fluid venting/halo formation

6 - 13

[years]

calculated as described in this paper

Vftot

total fluid volume

38 - 76

[km3]

calculated as described in this paper

ℜ

ratio of the amount of reactable hydrogen in the source fluid to the amount required to alter the rock

2.6 (EDM) 2.0 (sericitic)

[kg magmatic fluid/kg altered rock]

calculated as described in this paper

x

in the two dimensional single vein model, measurement of distance along vein.

variable

[m]

calculated as part of the finite difference simulation as described in this paper

xmax

in the two dimensional single vein model, the position along length of vein where alteration halo width goes to zero (same as R in the spherical model)

1800 -2000

[m]

estimated from drill core data

z

width of alteration halo measured normal to vein

variable

[m]

calculated in the finite difference simulations as described in this paper

zmax

maximum width of alteration halo. zmax occurs at x = 0

0.035 (EDM)

[m]

measured from drill core

0.025 (sericitic)

23

The implication of this observation is that the rock properties of the Butte Quartz Monzonite control the fracture pattern, and the fractures in the monzonite determine vein distribution. Maximum halo width zmax, taper, and radius of pervasive alteration zone, p

Halo width measurements were taken from a skeleton collection of the Deep Diamond Drill Core from drill holes two, six, and ten stored at the Museum of Mining in Butte, Montana. Halos were measured from the edge of the fracture (defined by the area dominated by quartz) to the outermost edge of the alteration halo. We estimate the uncertainty on each halo measurement is + or - 1 mm. Measurements of EDM and sericitic alteration halos were kept separate because cross-cutting relationships and the different temperatures of formation suggest that these different vein types were formed from different expulsion events separated in time (Geiger et al., 2002). Halo widths were measured around only around veins 1-3 mm in width. Veins of this size represent about 70% of all veins observed in the drill core. The widest veins fill quickly because vein mineral deposition is proportional to flow rate through the vein, which depends on the cube or fourth power of vein width. Thus the wide veins probably formed early and filled quickly, so the bulk of flow is through fractures in the 1-3 mm range. Figure 8 shows the measured halo widths around veins 1 to 3 mm wide as a function of distance from the fluid source. The minimum radius of pervasive alteration, p, is defined by the region in drill core 2 in which no distinct halos are visible (figures 6 and 7).

Drill hole 10 is in its restored position as illustrated in

24

Figure 8: Measured and calculated alteration halo widths versus distance from the fluid source (figure 4) for expelled volatile fluid mass of 30 billion tons for EDM alteration and 23 billion tons for sericitic alteration. Circles represent halo widths measured from core. Maximum halo width determined from vein spacing is shown as horizontal dashed lines; the radius of pervasive alteration, p, is the distance to which the calculated halo width equals maximum halo width. Each curve represents a different duration of venting: in the EDM case, the curves shown are for venting duration of 6.5, 13.0, and 26.0 years. In the sericitic halo case, curves shown are for venting duration of 3.1, 6.2, and 12.4 years. The results which best match measured alteration halos are shown in bold. In this case, duration of venting is 13.0 years for EDM alteration and 6.2 years for sericitic alteration.

25

figure 6. The horizontal line at the top with the dashed extension shows the fracture spacing used in the finite difference calculation. The spacing is based on the average fracture density. Thus at the edge of the model zone of pervasive alteration half the vein halos have merged. The near vertical line shows the results of the finite difference models we describe later. Consumption of reactable hydrogen, G

The amount of reactable hydrogen added during alteration can be determined by comparing unaltered mineral assemblages to altered mineral assemblages. The method is described by Cathles (1991). Application to the Butte system is made here. Table 2 shows the mineral composition of the fresh host rock and alteration halos at Butte. The sericitic alteration halos show two distinct zones, the "grey sericitic" and the "sericite with remnant biotite". For example, the innermost 0.01 meters of halo show grey sericitic alteration, and the outermost 0.005 meters show sericite with remnant biotite alteration (Geiger et al 2002). We assume that this 2:1 relationship between the two zones is similar in every sericite vein. Therefore, the average sericite assemblage used is a 2:1 weighted average of grey sericitic assemblage and sericite with remnant biotite assemblage. Rather than using solid solutions, we describe the rock composition in terms of a set of minerals with end- member chemical compositions. The end member mineral assemblage must contain the same number of minerals as the basis species we use to determine alteration, and all basis species are contained within at least one mineral. Table 3 shows the set of minerals used to describe the two types of alteration halos.

26

Table 2: Mineral composition of fresh Butte Quartz Monzonite and the alteration halo assemblages

Mineral

Butte Quartz Monzonite [mineral mode] (Brimhall 1977)

Averaged Sericitic Alteration [mineral mode] (estimated from Geiger et al 2002)

EDM alteration [mineral mode] (estimated from Brimhall 1977)

Quartz

23

53

46

Orthoclase (K-spar)

22

17

6

Andesine (Plagioclase)

37

9

9

Andalusite

0

7

0

Sericite (Muscovite)

0

6

30

Anhydrite

0

3

0

Hornblende

9

0

1

Biotite

8

5

3

Pyrite and Magnetite

1

1

5

27

Table 3: Proxy mineral assemblages for fresh Butte Quartz Monzonite and alteration halo assemblages

Original Mineral

Mineral Proxy Representation

Butte Quartz Monzonite, [moles mineral/kg rock] (Brimhall 1977)

EDM alteration [moles mineral/kg rock] (estimated from Brimhall 1977)

Averaged Sericitic Alteration [moles mineral/kg rock] (estimated from Geiger et al 2002)

Quartz

Quartz

3.67

8.47

7.17

K-Spar

Microcline

0.74

0.57

0.20

Plagioclase

Anorthite

0.67

0.16

0.16

Plagioclase

Albite

0.68

0.16

0.17

Andalusite

Andalusite

0.00

0.50

0.00

Sericite or Muscovite

Muscovite

0.00

0.14

0.76

Anhydrite

Anhydrite

0.00

0.20

0.00

Hornblende

Tremolite

0.12

0.00

0.01

Biotite

Annite

0.09

0.06

0.03

Biotite

Daphnite

0.07

0.04

0.02

Magnetite

Magnetite

0.05

0.04

0.18

Pyrite

Pyrite

0.06

0.08

0.34

28

Table 4 shows the computed the difference between moles of minerals/kg rock in the unaltered rock and moles of minerals/kg rock in the halo assemblages. Here, minerals created during alteration have a positive sign, and minerals destroyed have a negative sign. The mineral alteration makes up column vector [A]. Table 5 shows the stoichiometric matrices for minerals created and destroyed during alteration. In a stoichiometric matrix [SM], each row represents a mineral (the same minerals in the same order as are represented in [A]), and each column represents a basis species. The stoichiometric matrix is based on dissolution reactions from table 5 by Geiger et al (2002). The moles of basis species added or lost from the fresh Butte Quartz Monzonite during alteration, column vector [B], equals the product of the transpose of the mineral stoichiometric matrix times and mineral alteration column vector (see Cathles 1991):

[ B ] = [ S M ][ A] T

(12)

Table 6 shows [B], the moles of basis species added or lost during alteration. Species with negative signs were lost from the monzonite during alteration; species with positive signs were added. In both types of alteration the amount of H+ added was larger in magnitude than the loss or gain of any other basis species. For this reason reactable hydrogen is considered to be the basis species which controls the alteration. For Butte, the values for the reactable hydrogen consumption are: G = 4.68 moles acid/kg rock for the EDM alteration, G = 3.55 moles acid/kg rock for the sericitic alteration.

29

Table 4: Calculated difference between alteration halo assemblages and Butte Quartz Monzonite (negative sign indicates minerals destroyed during alteration) Mineral

EDM alteration [moles mineral/kg rock]

Quartz

Averaged Sericitic Alteration [moles mineral/kg rock]

4.79

3.50

Microcline

-0.17

-0.54

Anorthite

-0.50

-0.50

Albite

-0.51

-0.51

Andalusite

0.50

0.00

Muscovite

0.14

0.76

Anhydrite

0.20

0.00

Tremolite

-0.12

-0.11

Annite

-0.03

-0.06

Daphnite

-0.03

-0.05

Magnetite

-0.01

0.13

0.02

0.28

Pyrite

30

Table 5: Stoichiometric matrices for alteration (based on mineral/fluid reactions listed in Geiger et al. 2002) a. EDM alteration SiO2 K+ Al+3 Na+ Ca+2 Fe+2 Mg+2 HS- H+ SO421 0 0 0 0 0 0 0 0 0 Quartz 3 1 1 0 0 0 0 0 -4 0 Microcline 2 0 2 0 1 0 0 0 -8 0 Anorthite 3 0 1 1 0 0 0 0 -4 0 Albite 1 0 2 0 0 0 0 0 -6 0 Andalusite 3 1 3 0 0 0 0 0 -10 0 Muscovite 0 0 0 0 1 0 0 0 -1 0 Anhydrite 8 0 0 0 2 0 5 0 -14 0 Tremolite 3 1 0 0 0 3 0 0 -10 0 Annite 3 0 2 0 0 5 0 0 -16 0 Daphnite 0 0 0 0 0 3 0 0 -4 0 Magnetite 0 0 0 0 0 1 0 1.75 0.25 0.25 Pyrite

HSO4- H20 0 0 0 2 0 4 0 2 0 3 0 6 1 0 0 8 0 6 0 12 0 2 0 -1

b. averaged sericitic alteration K+

SiO2

Al+3

Na+

Ca+2

Fe+2

Mg+2

HS- H+

H2 0

Quartz

1

0

0

0

0

0

0

0

0

0

Microcline

3

1

1

0

0

0

0

0

-4

2

Anorthite

2

0

2

0

1

0

0

0

-8

4

Albite

3

0

1

1

0

0

0

0

-4

2

Muscovite

3

1

3

0

0

0

0

0

-10

6

Tremolite

8

0

0

0

2

0

5

0

-14

8

Annite

3

1

0

0

0

3

0

0

-10

6

Daphnite

3

0

2

0

0

5

0

0

-16

12

Magnetite

0

0

0

0

0

3

0

0

-4

2

Pyrite

0

0

0

0

0

1

0

2

0.25

-1

31

Table 6: Basis species change during alteration compared to fresh Butte Quartz Monzonite (negative sign indicates basis species removed during alteration) Basis Species

EDM alteration [moles mineral/kg rock]

SiO2

Averaged Sericitic Alteration [moles basis species/kg rock]

1.49

2.46

K+

-0.06

0.34

Al+3

-0.32

-1.82

Na+

-0.51

-0.51

Ca+2

-0.55

-0.60

Fe+2

-0.24

1.50

Mg+2

-0.62

-0.24

HS-

0.03

-0.23

H+

4.68

3.55

< 0.01

0.00

0.20

0.00

-2.58

-1.01

SO42HSO4H2 0

32

Results

Fluid volume is determined from the volume of altered rock and the duration of venting is determined from finite difference simulations using parameters described above. (Parameter values are compiled in table 1). Duration of venting, fluid volume, and intrusion size are calculated separately for each type of alteration. The total mass of volatiles vented is determined from C∑ H + and the volume of 0

altered host monzonite. The computed halo taper is approximately linear from the edge of the zone of pervasive alteration to the porphyry radius (figure 8). Therefore, the volume of altered rock between p and R is about one-third the volume of all the rock in that zone. We estimate the total volume of altered rock as this volume plus the entire volume of rock within the calculated pervasive alteration zone. If the concentration of reactable hydrogen in the magamatic source fluid is 3.0 moles/kg fluid (table 7), the mass of magmatic fluid required to alter one kilogram of rock is 2.6 kilograms of fluid for EDM alteration and 2.0 kilograms of fluid for sericitic alteration. For an unaltered rock density ~2700 kg/m3, 3.0 ⋅ 1010 tons of magmatic fluid are required to produce the alteration observed at Butte if all alteration were EDM, and 2.3 ⋅ 1010 tons of magmatic fluid is required if all alteration were sericitic alteration. This mass of fluid would require a spherical source intrusion with a diameter of 7.6 or 7.0 km, assuming the intrusion expels 5 wt% volatiles. The duration of venting for fixed T, C∑ H + , and f is constrained by the halo 0

taper from the edge of the zone of pervasive alteration, p, to the porphyry radius, R. The duration of venting is determined by matching measured halo widths to widths predicted by the finite difference simulations. Duration of venting is adjusted until the

33

Table 7: Calculation of hydrogen yield from CHILLER calculations at 400 deg Celsius. From Mark Reed, personal communication H+ ClHCl H2S aq. SO2 aq. HSO4S(nat)

start 6.29E-02 3.08E-01 1.78E-01 2.61E-01

5000gBQM 7.81E-06 4.25E-01 2.72E-05 1.49E-02

Difference factor H+ yield 6.29E-02 1 6.29E-02 1.78E-01 2.46E-01

1 1

1.78E-01 2.46E-01

4.08E-01

1.26E-05

4.08E-01

1.75

7.15E-01

1.03E-02 2.28E-04 1.50E+00 0.00E+00

1.00E-02 1.50E+00

1 1.25

1.00E-02 1.88E+00 3.09E+00 moles of H+ gone by reaction

1.0 SO2 aq. + 1.0 H2O = 1.75 H+ + .75 SO4--

+

.25 HS-

S native + 1.0 H2O = .75 HS- + 1.25 H+ + .25 SO4--

34

calculated halo widths match as closely as possible the maximum halo widths. For initial acid concentration C∑ H + of 3 moles/kg fluid, diffusion concentration factor f of 0

0.6, EDM halo formation at 600 degrees C, and sericitic halo formation at 400 degrees C, the outer envelope of EDM halo widths is best matched if 3.0 ⋅ 1010 tons of magmatic fluid are vented over 13.0 years, and the outer edge of sericitic alteration halos are best matched if 2.3 ⋅ 1010 tons of magmatic fluid are vented over 6.2 years (figure 8). A summary of the results is shown in table 8.

35

Table 8: Results for the Butte, Montana porphyry copper system 89 - 141 tons magmatic fluid/m length vein/yr

Initial fluid flux Duration of halo formation

6 - 13 years

Mass of volatiles expelled

23 - 30 billion metric tons 176 - 232 km3

Volume of intrusion, if volatiles = 5 wt% Radius of intrusion, assuming a spherical intrusion

3.5 - 3.8 km

36

Discussion

The base case interpretation

The estimate of 23 to 30 billion tons of volatiles expelled is based on the volume of alteration and C∑ H + . The system is closed in the sense that at sufficient 0

distance the acidity of the magmatic fluids has been exhausted. We can estimate the total volume of alteration, and constrain the volume of vented volatiles, at least when halos are forming. This method is similar to previous methods, and it yields similar results. For example Dilles and Reed (2006 personal communication) estimate that at least ~20 billion tons of magmatic fluid was required to produce the observed Pittsmont Dome pre-main stage alteration. The greatest uncertainty in our estimate is in the concentration of reactable acid supplied by the magmatic volatiles. The rate of volatile expulsion is measured by the halo taper beyond the zone of pervasive alteration. Halo taper has not, to our knowledge, been interpreted before, although Steefel and Lichtner (1998a, 1998b) suggested it could be useful. Measured halo widths show that even within the model pervasive alteration zone, many veins have halo widths less than the average fracture spacing (figure 8). We interpret this to result from the fact that not all the fractures in the monzonite are major flow fractures. The fractures through which volatiles could move radially outward could be major flow conduits, but the fractures which cut these major conduit fractures and are perpendicular to the outward fluid flow direction would accommodate much weaker flow. The halos around these veins would therefore be much less developed. The veins that are radially directed develop the widest halos. Because the width of a halo

37

measures minimum venting duration, the radially directed halos are the ones that should be used to constrain the duration of venting. With this in mind the outer envelope of halo width was used to match calculated halo width profiles. The model curves indicate that the halo taper will be steep unless the duration of venting is very short. The distribution of data where the halo width approaches zero is sparse, but we know that EDM vein halo widths decrease to zero in a road cut about 2.3 km from our model fluid source. If the duration of venting were decreased until the halo taper just reached the road cut, the duration of venting would be much shorter (3.3 years if calculated from EDM alteration, and 2.6 years if calculated from sericitic alteration). By matching the calculated curves as closely as possible to the outer envelope of the measured halo widths, we obtain the longest duration of venting that can accommodate the known data. Parametric variations

We consider the base case the most likely interpretation, but had parameters been chosen differently, slightly different interpretations would have resulted. For example if the total reactable hydrogen concentration of the source were 6 moles H+ per kg rather than 3, half the mass of magmatic volatiles and an intrusion about 25% smaller in diameter would be required. The source concentration is probably less than 3 (table 7), so our base case intrusion diameter is probably a low estimate. Parameters that could affect our conclusions include: f, φ , τ , and D. The manner in which these parameters affect the volatile flux is shown by equation 8. Increasing f, φ or D, or decreasing τ will increase the venting rate and decrease the duration of venting. The effects of changes in f are shown in figure 9. Estimates of

38

φ and τ are average for monzonite, and we believe our estimates of D are reasonable within 10%. Most of the uncertainty is therefore in f. Figure 9 shows that as long as f ≥ 0.6 our conclusions are unaffected, but alteration could occur further from the zone

of pervasive alteration if f were significantly lower than 0.6. Total intrusion volume required by the Butte porphyry system

The total volume of volatile-rich intrusions needed to produce the entire Butte porphyry is much larger than what we estimate was required to produce the mineralization and alteration in the Pittsmont Dome. The intrusions must produce not only the Pittsmont dome alteration but also the Anaconda dome alteration (see figure 6). The Anaconda dome alteration volume is about the same at that associated with the Pittsmont dome, and so two intrusions ~7 km in diameter are required. Quartz/moly veins which have no alteration halos were formed after the EDM halos and before the sericite halos. Our fluid volume estimates are based on the volume of altered rock, so they do not include fluid expelled through the quartz/moly veins. The expulsion rate to form these veins must have been very rapid to avoid any halo development, and the volume of fluid which flowed through them could be equal or greater than the fluid expelled though the EDM or sericitic veins. Most of the mineralization occurred in the main stage, not the pre-main stage. If the volume of volatiles vented is proportional to the metals deposited, at least an order of magnitude more volatiles are required. Thus, the sum total of volatile-rich magmas that intruded the crust near the Butte porphyry system was probably more than 30 times that required by the Pittsmont Dome. An intrusion more than 21 km in diameter was probably needed to form all the mineralization at Butte.

39

Figure 9: Results of the spherical finite difference simulation of formation of alteration halos around veins for EDM (top) and sericitic (bottom) alteration, in which the fraction of reactable hydrogen concentration which drives diffusion, f was varied. Maximum halo width determined from vein spacing is shown as horizontal lines with dashed extensions. The solid black lines show how the halos develop over time. The topmost line in each graph shows the shape of the alteration halo after the total time of volatile venting, 13.0 years for EDM and 6.2 years for sericitic alteration. Earlier curves are at one-year intervals. Parameters are: C ∑ H + = 3 moles/kg fluid, f = 0.6, φ 0

= 1%, τ = 5, G = 4.68 moles reactable hydrogen/kg rock and temperature = 600 C for EDM alteration, G = 3.55 moles/kg rock and temperature = 400 C for sericitic alteration.

40

Implications

The volume of pervasive alteration indicates total fluid expulsion. The larger the volume of pervasive alteration, the greater the volume of volatiles expelled, and the greater the potential for large tonnages of mineralization. The fact that the alteration tapers rapidly away from the pervasive zone (figure 8) means that the pervasive alteration zone provides a good estimate of the total alteration. Unlike a volatile-rich intrusion that has lost its volatiles, pervasive alteration might be easily detected at depth by geophysical methods. Defining the volume of deep pervasive alteration using geophysical methods might be a good way to prioritize deep porphyry exploration targets. Changes in halo taper suggests that pre-main stage expulsion in the Pittsmont Dome system began slowly (EDM alteration), was followed by a period of very rapid expulsion (unaltered quartz/moly veins), and slowed towards the end of expulsion (sericitic alteration). A significant amount of time then must then have passed to allow the system to cool from the 600 °C, the temperature at which EDM halos were formed to 400 °C, the temperature at which sericitic halos were formed. A 2 km diameter body, for example, takes about ~8000 years to cool (Carslaw and Jaeger, 1959). The duration of volatile venting is much shorter than this, so the EDM and sericitic halos formed in very short expulsion events which were separated by a much larger amount of time. In fact, the venting rates we estimate are so short that they are best thought of as controlled eruptions, and control of the eruption is critical. If volatiles are expelled explosively, any mineralization will be dispersed across a broad landscape, as probably occurred during the eruption of Mt. Pinatubo in Mexico (Hattori, 1996).

41

Finally, it is interesting and perhaps useful that the essential characteristics of a porphyry system that are useful in exploration can be captured by so few variables. Just two observations are required: the mineralogy of altered versus unaltered rock and the diameter of pervasive alteration. The mineralogy is similar for most porphyrys, so the potential mineral tonnage of a porphyry is dependant primarily on volume of pervasively altered rock. The duration of venting is constrained by a third observation: the distance at which fractures cease to form relatively abundant veins with halos.

42

Summary and Conclusions

1. Finite difference simulations of vein halo formation suggest that, assuming that mineralogy of altered and unaltered rock is known, the size and the time to form a porphyry system can be estimated from just two observations: the radius of pervasive alteration and the distance at which vein halo width goes to zero. 2. Because halo width tapers linearly with distance from the zone of pervasive alteration, the volume of alteration is easy to estimate. The mineral potential of of a porphyry system should be related to the volume of magmatic volatiles expelled and thus is measured by the volume of alteration. Mapping the volume of pervasive alteration using geophysical techniques could provide a way to rank deep exploration targets. 3. The taper of alteration halos beyond the radius of pervasive alteration constrains the duration and rate of magmatic fluid venting. A longer duration/slower venting rate yields a steeper taper and a shorter duration/higher venting rate yields a shallower taper. 4. For Butte pre-main stage Pittsmont Dome alteration, at least 23 - 30 billion tons of magmatic fluid was vented in less than 20 years. Two short pulses of rapid volatile venting separated by thousands of years formed the EDM and sericite alteration halos. Volatile expulsion seems to have started slowly at the EDM stage, increased markedly when quartz/moly veins were formed, and then slowed again when the sericitic halos were formed.

43

APPENDIX

Matlab code for finite difference simulations

function [temp]= Execute(D,porob,porov,G,m,densb,H,taob,Ceq,FLUX) %% Calculates the position of the outermost edge of a reaction front normal %% to a vein, based on the diffusion of a single critical basis species. clear all %%------------Default Adjustable parameters---------------------------------------R = 1800; %m max radius of porphyry p = 1300; %m radius of pervasive alteration iterations = 1; porob = .01; %fract porov = 1;%fract H = .001; %m *Note! enter TOTAL diameter of vein- H is divided by 2 in the alogrithm to model as half-space taob = 5; %no units Ceq = 0; %mole/kg fluid Cv = 3; %moles h+/kg fluid DelCo = Cv-Ceq; %mole/m^3 densb = 2712 %kg/m^3 rock densw = 400 %kg/m^3 fluid f = .6 [EDMWidths,EDMDist,SRWidths,SRDist]= Widths(H) %%-------------Calc Altered Rock Volume---------------------------------Vporph = (2/3)*pi*R^3; %%m^3 Volume of porphyry VpervAlt = (2/3)*pi*p^3; %%m^3 Volume of pervasive alt ValtRest = (Vporph - VpervAlt)/3; %%m^3 Volume of alt, not pervasive VAlt = VpervAlt + ValtRest; %% m^3 total volume of altered rock

% %%------------------SrAdjustable Parameters----------------------------plotorder = 0; alttype = 'Sr'; temp = 600 %%deg C 44

taper = 1.6E-5; densw = (1 - (.001*temp))*(100^3/1000); D = 0.5*exp(-2551.02*((1/(273+temp))-(1/673)));%%m^2/yr G = 4.68; %%moles/kg rock WWR = (G/(DelCo))*(densb/densw) %%--------------------------Calculate EDM Fluid Volume-----------------------FLUIDVOLUME = VAlt*WWR %%Volume of fluid = vol of alt rock water/rock ratio FLUIDMASS = FLUIDVOLUME*densw; INTMASS = FLUIDMASS/0.05; INTVOLUME = (INTMASS/2600) INTRAD = ((3/(4*pi))*INTVOLUME)^(1/3) %--------------------EDM Parameters for plotting/advection-----------------noshells = 200; %timesteps %%number of timesteps/cells timesteps = 54167; % deltat = 2.4E-4; % maxtime = timesteps*deltat; % plotorder = plotorder + 1;% ro = R/(noshells^(1/3)); n = 2:1:noshells+1; XAXIS(1) = 0; XAXIS(n) = ro*((n-1).^(1/3)); cells = length(XAXIS) Vo = FLUIDVOLUME/timesteps; InFLUX = Vo/deltat %-------- Initial Conditions for Numeric Solution--------------------------DELTAC = zeros(1,cells); DELTAC(1) = DelCo; ZH = 0.001*ones(1,cells); %sqrt(D*deltat/(1+(G/(DelCo))))*ones(1,cells); DELTAZ = zeros(1,cells); DDELTAC = zeros(1,cells); T = zeros(1,cells); %%%-----NUMERIC SOLUTION FOR DECREASING VEIN CONC---------plotinterval = round(((timesteps)/10)) %how often to plot a timestep for t = 1:timesteps count = timesteps - t;

45

FLUIDVOLUME = t*Vo; iterations = max(1,double(int8(100000/t)));%number of iterations in a timestep [DELTAC,DELTAZ,ZH,spacing]= calcfront(t,D,porob,porov,G,taob,... Ceq,DDELTAC,DELTAC,DELTAZ,ZH,T,DelCo,deltat,cells,maxtime,H,... plotinterval,iterations,densb,densw,f,plotorder,noshells); %Call calc front (numeric calc loop) if (mod(count,plotinterval) == 0); count, if t == timesteps, [k]= plotf(XAXIS,ZH,H,D,maxtime,INTVOLUME,FLUIDVOLUME,plotorder,alttype,... t,deltat,spacing,INTRAD,noshells,EDMWidths,EDMDist,SRWidths,SRDist,InFLUX, DelCo,f); %call plots end end %end if end %%end for t %%------------------Sr Adjustable Parameters----------------------------plotorder = 2; alttype = 'Sr'; temp = 400 %%deg C taper = 8.0E-6; densw = (1 - (.001*temp))*(100^3/1000); D = 0.5*exp(-2551.02*((1/(273+temp))-(1/673)));%%m^2/yr G = 3.55; %%moles/kg rock WWR = (G/(DelCo))*(densb/densw) %%--------------------------Calculate Sr Fluid Volume-----------------------FLUIDVOLUME = VAlt*WWR; %%Volume of fluid = vol of alt rock water/rock ratio FLUIDMASS = FLUIDVOLUME*densw; INTMASS = FLUIDMASS/0.05; INTVOLUME = (INTMASS/2600); INTRAD = ((3/(4*pi))*INTVOLUME)^(1/3) %% ------------Calculated Parameters for plotting/advection-----------------noshells = 200 %timesteps %%number of timesteps/cells timesteps = 24900; deltat = 2.5E-4 % maxtime = timesteps*deltat

46

ro = R/(noshells^(1/3)); n = 2:1:noshells+1; XAXIS(1) = 0; XAXIS(n) = ro*((n-1).^(1/3)); cells = length(XAXIS) Vo = FLUIDVOLUME/timesteps; InFLUX = Vo*deltat; %%--------Initial Conditions for Numeric Solution--------------------------DELTAC = zeros(1,cells); DELTAC(1) = DelCo; ZH = 0.001*ones(1,cells); %sqrt(D*deltat/(1+(G/(DelCo))))*ones(1,cells); DELTAZ = zeros(1,cells); DDELTAC = zeros(1,cells); T = zeros(1,cells); %%-------NUMERIC SOLUTION FOR DECREASING VEIN CONC-----------plotinterval = round(((timesteps)/10)) %how often to plot a timestep for t = 1:timesteps count = timesteps - t; FLUIDVOLUME = t*Vo; iterations = max(1,double(int8(100000/t)));%number of iterations in a timestep [DELTAC,DELTAZ,ZH,spacing]= calcfront(t,D,porob,porov,G,taob,... Ceq,DDELTAC,DELTAC,DELTAZ,ZH,T,DelCo,deltat,cells,maxtime,H,... plotinterval,iterations,densb,densw,f,plotorder,noshells); %Call calc front (numeric calc loop) if (mod(count,plotinterval) == 0); count, if t == timesteps, [k]= plotf(XAXIS,ZH,H,D,maxtime,INTVOLUME,FLUIDVOLUME,plotorder,alttype,... t,deltat,spacing,INTRAD,noshells,EDMWidths,EDMDist,SRWidths,SRDist,InFLUX, DelCo,f); %call plots end end %end if end %%end for t end %%% end function

47

function [DELTAC,DELTAZ,ZH,spacing]= calcfront(t,D,porob,porov,G,taob,... Ceq,DDELTAC,DELTAC,DELTAZ,ZH,T,DelCo,deltat,cells,maxtime,H,... plotinterval,iterations,densb,densw,f,plotorder,noshells); %% calculates the position of a rection front normal to the vein in 2D spacing = 0.023*ones(1,length(ZH)); %%flat DDELTAC=(2*DELTAC*f*D*porob*deltat./(ZH*H*porov*taob)); for checkZ = 1:noshells+1; if (ZH(checkZ)> spacing(checkZ)); DDELTAC(checkZ) = 0; end end DELTAZ = (DDELTAC*densw/(2*densb*G)).*(H - (porob*DELTAZ)); ZH = ZH + DELTAZ; DELTAC = DELTAC - DDELTAC;

DELTAC(2:cells) = DELTAC(1:cells-1); DELTAC(1) = DelCo; T(t) = t; end %% end function

48

function [EDMWidths,EDMDist,SRWidths,SRDist]= Widths(H) H = 2; EDMT = 0.001*[20,10,16,1.5,1,1,5,5,2,1,35,35,5,12,5,3,10, 27,15,2,3,9,3,21,15,2,9,4,5,4,3,12,10,3]; EDMWidths = EDMT' EDMDist = [1470.889687,1469.382847,1422.086701,1422.086701, 1422.086701,1422.086701,1405.526913,1376.935284,1361.592446, 1298.166284,1271.739984,1271.739984,1167.712325,1097.748636, 1092.672653,1086.105215,1085.209787,1063.132766,1057.169469, 1057.169469,713.10754,660.5615976,650.2255135,926.8588173, 934.5820105,901.9528576,907.3300534,920.7447171,925.1998742, 941.9169355,944.285086,946.9977379,979.4297852,979.4297852]' SRT = 0.001*[0.5,5,2,2,9,15,10,25,14,25,5,10,20,10,15,20 7,4,2,10,10,12,4,4,10,3]; SRWidths = SRT' SRDist = [1564.680423,1546.576234,1507.365717, 1469.382847,1357.080722,1345.652841,1298.166284,1207.249215, 1207.249215,1188.673035,1188.673035,1085.209787,1081.628389, 1052.495064,982.3974372,936.9181713,905.6159942,905.6159942, 902.379947,902.379947,959.9581573,965.8071877,974.2172045, 980.1526616,1137.194232,1198.295189]' end %% end function

49

function [k]= plotf (XAXIS,ZH,H,D,maxtime, INTVOLUME, FLUIDVOLUME, plotorder,alttype,t,deltat,spacing,INTRAD,noshells,EDMWidths,EDMDist,SRWidths, SRDist,InFLUX,DelCo,f); %subplot(2,3,plotorder) plot(XAXIS,ZH,'-g') hold on % if (plotorder < 4) % subplot(2,3,plotorder) % plot(EDMDist,EDMWidths,'ok') % plot(XAXIS,spacing,'--') % else % subplot(2,3,plotorder) plot(SRDist,SRWidths,'ok') plot(XAXIS,spacing,'--') %end %hold off axis([0 2300 0 0.04]) xlabel('Distance from source (m)','FontSize',14) ylabel('Width of halo (m)','FontSize',14) %title%({[alttype,' alteration']});... %['Timestep = ',num2str(deltat),' years'] title({['Fluid Volume = ',num2str(FLUIDVOLUME/(1000^3)),' km^3 fluid'];... ['Intrusion Radius = ',num2str(INTRAD/(1000)),' km'];... ['t_{venting} = ',num2str((maxtime)),' years',' DelCo =',num2str(DelCo)]}) %['Initial Flux = ',num2str(round(InFLUX/1000)),' ton fluid/m length vein/yr']})%%;... pause(.0001) k = inFLUX;

end %%end function

50

function [G]= BasisChange(SM,A) %%B is the colum vector of the total change in basis species during %%alteration. SM is the Stoiciomentric Matrix. A is the alteration (in %%minerals. SM and A must be entered. %%For sericitic SM = [1 0 0 3 1 3 2 0 2 3 0 1 3 1 1 0 0 0 0 0 0 8 0 0 3 1 0 3 0 2

0 0 0 1 0 0 0 0 0 0

0 0 1 0 0 0 0 2 0 0

0 0 0 0 0 1 3 0 3 5

0 0 0 0 0 0 0 5 0 0

0 0 0 0 0 2 0 0 0 0

0 -10 -8 -4 -4 0.25 -4 -14 -10 -16

0, 6, 4, 2, 2, -1, 2, 8, 6, 12]

%%For EDM % SM = %[1 0 %3 1 %2 0 %3 0 %1 0 %3 1 %0 0 %8 0 %3 1 %3 0 %0 0 %0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 1 0 0 0 1 2 0 0 0 0

0 0 0 0 0 0 0 0 3 5 3 1

0 0 0 0 0 0 0 5 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 1.75

0 -4 -8 -4 -6 -10 -1 -14 -10 -16 -4 0.25

0 0 0 0 0 0 0 0 0 0 0 0.25

0 1 2 1 2 3 0 0 0 2 0 0

SMT = SM'; ASR = [6.25E-01,2.53E-01,-1.98E-01,-2.01E-01,1.33E-02, 3.63E-02,1.86E-02,-1.12E-02,-3.47E-02,-2.43E-02]

%%%%For EDM % AEDM = [10179.61,814.86,-1613.27, -1639.97, 67.51,731.34, % 397.88,-6.77, -139.28, -91.27];

51

0 0 0 0 0 0 1 0 0 0 0 0

0, 2, 4, 2, 3, 6, 0, 8, 6, 12, 2, -1]

% AEDM = [13082.55659,-459.8857471, -1377.744133, -1394.570755, % 1360.691145,389.0045396,545.3944469, -337.9002179, -92.71174544, % -72.26956548,-66.72713138, 125.2708785]

% %%%%%For Averaged SR ASR = [6.25E-01,2.53E-01,-1.98E-01,-2.01E-01,1.33E-02 3.63E-02,1.86E-02,-1.12E-02,-3.47E-02,-2.43E-02]; % BaSp = {'Si' 'K' 'Al' 'Na' 'Ca' 'Fe' 'Mg' 'HS' 'H' 'H20'} % BEDM = SMT*AEDM % GEDM = max(abs(BEDM))

BSR = SMT*ASR GSR = max(abs(BSR)) %% % BGSR = SMT*AGSR % GreyGSR = max(abs(BGSR)) end %%end function

52

REFERENCES

Brimhall, G.H (1977); Early fracture controlled disseminated mineralization at Butte, Montana. Economic Geology, v. 72, p 37-59. Carslaw, H.S. and J.C. Jaeger, 1959, Conduction of Heat in Solids, 2nd ed., Oxford University Press, London, 510 p. Cathles, L.M., 1979, Predictive capabilities of a finite difference model of copper leaching in low grade industrial sulfide waste dumps, Mathematical Geology, v. 11, p. 175 - 191. ——— 1991, The importance of vein salvaging in controlling the intensity and character of subsurface alteration in hydrothermal systems, Economic Geology, v. 86 p 446-471. Cathles, L.M., and J.A.Apps, 1975 A model of the dump leaching process that incorporates oxygen balance, heat balance, and air convection, Metallurgical Transactions, v. 6B, p. 617- 623. Geiger, S, R. Haggarty, J.H. Dilles, M.H. Reed, and S.K. Matthai, 2002. New insights from reactive transport modeling: the formation of the sericitic vein envelopes during early hydrothermal alteration at Butte, Montana: Geofluids, v. 2, p 185-201. Guilbert, J.M and C.F. Park, 1986: Ore Deposits, W.H Freeman and Company, New York, 985 p. Hattori, K., 1996, Occurrence and origin of sulfide and sulfate in the 1991 Mount Pinatubo eruption products, in Newhall, C.G., and Punongbayan, R.S., eds, Fire and Mud: Eruptions an lahars of Mount Pinatubo, Philippines: Quezon City, Philippines, Philippine Institute of Voloconalogy and Seismology, p 807-824. Henley, R.W. and Berger, B.R., 2000, Self-ordering and complexity in epizonal mineral deposits: Annual Reviews of Earth and Planetary Science, v. 28, p 699-719. Levenspiel, O., 1967, Chemical reaction engineering, John Wiley & Sons, New York p.501. Meyer, C., E.P. Sea, and C.C. Goddard,1968; Ore Deposits at Butte, Montana. In: Ore Deposits of the United States (1933-1967) (Ed. Ridge. D.J.) Graton-Sales Volume, 2, 1363-1416. American Institute of Mining, Metallurgical, and Petroleum Engineers, New York.

53

Richards, 2003; Tectono-magmatic precursors for porphyry Cu-(Mo-Au) deposit formation. Economic Geology, 98 pp 1515-1533. Sales, R.H. and C. Meyer (1949); Results for preliminary studies of vein formation at Butte, Montana. Economic Geology, 44, pp 465-484. Singer, D. A. (1995); World class base and precious metal deposits; a quantitative analysis, Economic Geology, 90, pp 88-104. Steefel, C.I., and P.C. Lichtner,1998a, Multicomponent reactive transport in discrete fractures. I. Controls on reaction front geometry, Journal of Hydrology, v. 209 p 186 199. ——— 1998b, Multicomponent reactive transport in discrete fractures. II. Infiltration of hyperalkaline groundwater at Maqarin, Jordan, a natural analogue site, Journal of Hydrology, v. 209 p 200 - 224.

54

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & Close