Chemical mass balance equations for open-system magma chamber processes that result in crystal zoning

Journal of Volcanology and Geothermal Research 374 (2019) 181–196

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

Chemical mass balance equations for open-system magma chamber processes that result in crystal zoning Koshi Nishimura Natural Science Laboratory, Toyo University, 5-28-20 Hakusan, Bunkyo-ku, Tokyo 112-8606, Japan

a r t i c l e

i n f o

Article history: Received 5 October 2018 Received in revised form 13 February 2019 Accepted 18 February 2019 Available online 10 March 2019 Keywords: Open-system magmatic process Crystal isotope stratigraphy Geochemical modeling Sr isotopes Magma chamber

a b s t r a c t Recent advances in micro-analysis have revealed that isotopic zoning profiles in igneous crystals reflect progressive changes in liquid composition during open-system magmatic processes. Consequently, the isotope ratio of the integrated (bulk) crystals and that of liquid are inevitably different. This study develops a mass balance model to describe the compositional relationships between the bulk crystals (total suspended crystals), liquid, magma (crystals plus liquid), and crystal rims during open-system magma chamber processes. The model incorporates the effects of concurrent magma recharge (or mixing), assimilation, magma extraction, and fractional crystallization. The mass balance differential equations for suspended crystals and liquid are solved simultaneously. The analytical solution of the equations gives a quantitative account of the evolution paths of trace elements and isotopes within total suspended crystals, liquid, magma, and crystal rims. The isotopic zoning profiles in phenocrysts vary markedly depending on the relative rates of recharge, assimilation, magma extraction, and fractional crystallization. The model can reproduce the 87Sr/86Sr zoning patterns in natural feldspars from different magma systems. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Until the 1990s, isotopic analyses of igneous rocks were largely restricted to powdered whole rock samples, except for the isotopic studies of separated mineral grains (e.g., Doe, 1970; Deloule et al., 1986). Whole rock isotopic analyses revealed isotopic heterogeneity at the scale of rock samples and provided evidence of open-system magma processes such as magma mixing (e.g., Halliday et al., 1980; Dickin and Exley, 1981; Arth et al., 1986; Holden et al., 1987) and combined assimilation and fractional crystallization (e.g., Taylor Jr., 1980; DePaolo, 1981; Hildreth and Moorbath, 1988; Davidson and Wilson, 1989). Studies that identified isotopic disequilibrium between whole rock (or matrix) samples and mineral separates (i.e., grain-scale disequilibrium) placed further constraints on the timings of assimilation and magma mixing (e.g., Cortini and van Calsteren, 1985; Geist et al., 1988; Mazzone and Grant, 1988; Simonetti and Bell, 1993). Recent micro-analyses of isotopes have revealed isotopic heterogeneity at the scale of sub-grains (e.g., Feldstein et al., 1994; Davidson and Tepley, 1997; Davidson et al., 1998, 2001, 2007; Knesel et al., 1999; Tepley et al., 1999, 2000; Waight et al., 2000; Ramos et al., 2005; Charlier et al., 2006, 2008). Isotopic variations recorded from core to rim of a single mineral grain reflect progressive changes in the composition of the liquid from which the mineral crystallized (e.g., Davidson et al., 1998, 2007). E-mail address: [email protected].

https://doi.org/10.1016/j.jvolgeores.2019.02.012 0377-0273/© 2019 Elsevier B.V. All rights reserved.

Consequently, the isotope ratios of liquid and of integrated (bulk) crystals are inevitably different in an open-system case in which the isotope composition of liquid changes over time (Davidson et al., 2007). However, the quantitative evolution paths of trace elements and isotopes within bulk crystals (total suspended crystals), liquid, and magma (crystals plus liquid) remain poorly understood because of their complexity. This study seeks to develop a mass balance model for trace element and isotope behavior during the opensystem magma chamber processes associated with the development of crystal zoning. The chemical evolution of an open-system magma chamber is basically controlled by competition between mixing and fractionation processes (Albarède, 1995). DePaolo (1981) developed a systematic geochemical model for simultaneous assimilation and fractional crystallization (AFC) that is the basis of modeling more complex open-system magma chamber processes (e.g., DePaolo, 1985; Hagen and Neumann, 1990; Spera and Bohrson, 2001, 2002, 2004; Nishimura, 2012, 2013; Bohrson et al., 2014). The growth of zoned crystals, however, is not included in the DePaolo formulation. Recently, Bohrson et al. (2014) provided a thermodynamic model (the Magma Chamber Simulator, MCS) for energy-constrained open-system magma chamber processes that enables calculation of major and trace element concentrations, isotopic ratios, masses, and thermal constraints for liquid + crystals ± fluid. This model assumes that when crystals form in equilibrium with liquid in the magma chamber, they are incrementally removed to the cumulate

182

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

reservoir and, once removed, they do not interact chemically or thermally with the magma chamber. Thus, in their model, zoned crystals do not develop. Nishimura (2012) developed a mass balance model describing trace element and isotopic zoning in crystals during simultaneous assimilation and fractional crystallization (AFC). The present study expands this model to more general and complex open-system magma chamber processes that include magma recharge (or mixing), assimilation, magma extraction, and fractional crystallization. This new model can be used to reproduce the trace element and isotopic evolution paths of total suspended crystals, liquid (matrix or glass), magma (whole rock), and crystal zoning during various open-system magma chamber processes.

concentration of suspended crystals, and the isotopic ratio of suspended crystals can be respectively expressed as _ aCa þ M _ i C i −M _ c DC l −M _ e ð1−ϕÞC l ¼ C l dM l þ Ml dC l ; M dt dt

_ a C a εa þ M _ i C i ε i −M _ c DC l ε l −M _ e ð1−ϕÞC l ε l ¼ C l ε l dMl þ Ml ε l dC l þ C l M l dε l ; M dt dt dt

ð2Þ _ s C x −M _ e ϕC x ¼ C x dMx þ Mx dC x ; _ c DC l −M M dt dt

The present study develops a geochemical model for chemical zoning of crystals in an open-system magma chamber in the style of the model development used by DePaolo (1981). Fig. 1 shows the schematic illustration of the general open-system magma chamber processes, including melt influx (recharge or mixing), the growth of zoned crystals, magma extraction, assimilation of anatectic wall-rock melt, and partial settling of crystals. Zoned crystals are assumed to be homogeneously distributed by vigorous convection in a magma chamber. The extracted (erupted) magma is therefore composed of liquid and suspended crystals. The crystals are partially separated from the magma body to the cumulate pile at the base of the magma chamber due to gravity, even if the convective velocity is much greater than the crystal settling velocity (Martin and Nokes, 1988). The mass balance differential equations for an element in liquid, the isotopic ratio of liquid, the elemental

ð3Þ

and _ s C x εx −M _ e ϕC x εx ¼ C x εx dMx þ Mx εx dC x þ C x M x dε x ; _ c DC l εl −M M dt dt dt

2. Model

ð1Þ

ð4Þ

_ a is the mass assimilation rate (mass/unit time), Ca is the where M element's concentration in the melt (liquid) derived from wall-rock melt_ i is the mass influx rate (mass/unit time), Ci is the element's concening, M

_ c is the mass crystallization rate (mass/unit tration in the injected liquid, M time), D is the bulk crystal/liquid partition coefficient for the element, Cl is _ e is the element's concentration in the liquid in the magma chamber, M

the mass extraction rate (mass/unit time), ϕ is the weight fraction of the suspended crystals in the magma chamber relative to the total mass of magma, Ml is the mass of the liquid in the magma chamber, εa is the isotope ratio of the wall-rock, εi is the isotope ratio of the injected liquid, _ s is the rate of εl is the isotope ratio of the liquid in the magma chamber, M gravitational separation of crystals from magma (mass/unit time), Cx is the mean concentration of the element in suspended crystals, Mx is the

Magma extraction M

e

Growth of zoned crystals Mc Assimilation of anatectic wall-rock melt M

a

Crystal separation M s Cumulate pile

Melt influx M

i

Fig. 1. Schematic illustration of an open-system magma chamber model incorporating the effects of simultaneous melt influx (recharge or mixing), the growth of zoned crystals, magma extraction, assimilation, and the partial settling of crystals. A capital ‘M’ with a subscript letter denotes the mass of the relevant process.

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

183

Table 1 Symbols used in the formulation. Symbol

Description

_a M _i M

Rate of assimilation of anatectic wall-rock melt

kg s−1

Rate of liquid influx (recharge or mixing)

kg s−1

_c M _e M

Rate of crystallization

kg s−1

Rate of magma extraction

kg s−1

_s M ra ri re Ml Mx Mm Mm0 Mc Ms ϕ Fm Ca Ci Cl Cx Cm C0l C0x Cm0 Cα−rim D Dα εa εi εl εx εm ε0l ε0x εm0 εα−rim h−o

Rate of crystal separation from magma by gravitational settling Ratio of rate of assimilation to rate of crystal separation Ratio of rate of liquid influx to rate of crystal separation Ratio of rate of magma extraction to rate of crystal separation Mass of liquid in magma body Mass of total suspended crystals Mass of magma (liquid + suspended crystals) Mass of initial magma Total crystallized mass Total mass of crystals separated from the magma by gravitational settling Weight fraction of suspended crystals relative to total magma mass Ratio of magma mass to initial magma mass Trace element concentration of anatectic melt (assimilant) Trace element concentration of injected liquid Trace element concentration of liquid in magma body Trace element concentration of total suspended crystals Trace element concentration of magma Trace element concentration of initial liquid in magma body Trace element concentration of initial suspended crystals Trace element concentration of initial magma Trace element concentration of rim of mineral α Bulk crystal/liquid partition coefficient Crystal/liquid partition coefficient of mineral α Isotope ratio of anatectic wall-rock melt Isotope ratio of injected liquid Isotope ratio of liquid in magma body Isotope ratio of total suspended crystals Isotope ratio of magma Isotope ratio of initial liquid in magma body Isotope ratio of initial suspended crystals Isotope ratio of initial magma Isotope ratio of rim of mineral α Constants defined in Eqs. (10)–(13) and (19)–(22)

kg s−1

total mass of suspended crystals, and εx is the mean isotopic ratio of suspended crystals (Table 1). One of the commonest forms of fractional crystallization is gravitational separation of crystals that are in instantaneous chemical equilibrium with the liquid as they form, but are immediately removed from contact with the liquid. In that case, the crystallized mass Mc is equivalent to the gravitationally separated mass Ms. However, in the present model, Mc is generally different from Ms because of the assumption of the partial settling of suspended crystals. The growth of zoned crystals also causes fractional crystallization (i.e., slow diffusivity in crystals causes the failure of the early formed crystal interior to react with the liquid that remains). Partial settling of suspended crystals has no effect on liquid mass and liquid composition; therefore, the mass balance equations for liquid do not include the mass of settled crystals Ms (Eqs. (1) and (2)). In contrast, the mass and elemental compositions of suspended crystals are affected by the partial settling of crystals as well as by crystallization. Therefore, the mass balance equations for suspended crystals include both Mc and Ms (Eqs. (3) and (4)). The differential equations for crystal zoning during AFC (Nishimura, 2012) corre_ i ¼ 0 and M _ e ¼ 0 . DePaolo (1981) spond to the extreme case at M defined the constant ratio of the rate of assimilation to the rate of crystal _ a =M _ s ) in order to analytically solve the differential separation (r a ¼ M equations for AFC. Similarly, the present study defines the constant _i ratio of the rate of liquid influx to the rate of crystal separation (r i ¼ M

_ s ) and that of the rate of magma (liquid + suspended crystals) ex=M _ e =M _ s ). It should be traction to the rate of crystal separation (r e ¼ M noted that the ratios ra, ri, and re can vary during the open-system evolution of a magma body (e.g., the ra variations reported by Bohrson and

Unit

kg kg kg kg kg kg

ppmw ppmw ppmw ppmw ppmw ppmw ppmw ppmw ppmw

Spera, 2001). However, constant ratios can be used as an approximation for a short crystallization period and can be estimated from variations in trace elements and isotopes, as for the AFC model (e.g., Powell, 1984; Mantovani and Hawkesworth, 1990; Young et al., 1992; Aitcheson and Forrest, 1994; Roberts and Clemens, 1995; Caffe et al., 2002). Changes in the ratios can be incorporated into the calculations as described in Section 3.3. The instantaneous rate of change in the mass of a magma body, affected simultaneously by assimilation, liquid magma influx, magma (liquid + suspended crystals) extraction, and gravitational separation of crystals, is given as dMm _ s: ¼ ðr a þ r i −r e −1ÞM dt

ð5Þ

_aþ where Mm is the mass of magma body. When ra + ri − re = 1 (i.e., M _ _ _ Mi ¼ Ms þ Me ), added material (assimilant + injected liquid) is balanced by an equivalent mass of subtracted material (gravitationally separated crystals + extracted magma), meaning that the magma mass is constant. The analytical solutions to the system of differential equations (Eqs. (1)–(4)) for an element can be separated conveniently into two cases: one for ra + ri − re ≠ 1, and another for ra + ri − re = 1 (the detailed derivation can be found in the Supplementary material S1). Assuming that ϕ, D, Ca, εa, Ci, and εi are constant, the solutions can be written as follows: For ra + ri − re ≠ 1 (general case),
 C l ¼ h=i þ C 0l −h=i F m −i ;

ð6Þ

184

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196


 j þ iC 0l ε0l − j F m −i
 εl ¼ ; h þ iC 0l −h F m −i

ð7Þ

and o¼

ð1−ϕÞðϕr e þ 1ÞD : ðϕr e þ 1ÞϕD−ð1−ϕÞ

ð22Þ



1 0 kD C 0l −h=i kD C 0l −h=i h h 0 −i @ A F m −k ; ð8Þ F m þ Cx − D þ C x ¼ D− i i i−k i−k

For all cases, the trace element concentration (Cm) and isotopic ratio (εm) of the magma (crystals plus liquid) can be respectively expressed as

and

and

2 j εx ¼ 4 D− i




 kD C 0l ε 0l − j=i i−k

1 3 0 kD C 0l ε0l − j=i j 0 −i 0 A F m −k 5=C x ; F m þ @C x ε x − D þ i i−k




ð9Þ where C 0l is the initial concentration of the element in the liquid, ε 0l is the initial isotope ratio of the liquid, Cx0 is the initial concentration of the element in the suspended crystals, εx0 is the initial isotope ratio of the suspended crystals, Fm is the ratio of magma mass to the initial magma mass (Fm ≡ Mm/Mm0),

εm ¼

ð1−ϕÞC l εl þ ϕC x εx : Cm

ð10Þ

i¼

ðϕr a þ ϕr i þ 1−ϕÞD þ ð1−ϕÞr e þ 1; ð1−ϕÞðr a þ r i −r e −1Þ

ð11Þ

j¼

ra C a εa þ ri C i εi ; ð1−ϕÞðr a þ ri −r e −1Þ

ð12Þ

ϕr a þ ϕr i þ 1−ϕ ; ϕðr a þ r i −r e −1Þ

ð13Þ

and

ð23Þ

ð24Þ

The trace element concentration of the rim of a crystal α can be expressed as C α−rim ¼ Dα C l ;

ð25Þ

where Dα is the crystal/melt partition coefficient. The isotopic composition of the crystal rim is identical to that of the liquid: εα−rim ¼ εl :

ra C a þ ri C i ; h¼ ð1−ϕÞðra þ r i −re −1Þ

k¼

C m ¼ ð1−ϕÞC l þ ϕC x ;

ð26Þ

The initial concentrations of magma, liquid, and suspended crystals must be consistent with the mass balance. The present study assumes that the initial suspended crystals grew via surface equilibrium crystallization (i.e., fractional crystallization by incomplete chemical reaction, producing zoned crystals) from a liquid with an initial magma composition Cm0 (Nishimura, 2009, 2012, 2013) except in the case that the liquid or crystal maintains a constant composition through time (hereafter referred to as compositional steady-state). The initial liquid composition can be derived from the Rayleigh fractionation equation (Rayleigh, 1896) as C 0l ¼ C 0m ð1−ϕÞðD−1Þ :

ð27Þ

From the mass balance, C0x can be expressed as

and for ra + ri − re = 1,
 Ms C l ¼ l þ emMm C 0l −l ;

ð14Þ

C 0x ¼

h
i Ms εl ¼ n þ emMm C 0l ε0l −n =C l ;

ð15Þ




ϕre þ1 Ms  ϕre þ1 Ms Ms C x ¼ lD þ o C 0l −l e− ϕ Mm −emMm þ C 0x −lD e− ϕ Mm ;

ð16Þ

The isotopic composition of the initial suspended crystal is assumed to be identical to that of the initial liquid (ε0x = ε0l ). If we consider an initial magma chamber with a steady-state composition, however, suspended crystals are not chemically zoned because the melt composition remains constant with time. In this case, the initial compositions of liquid, total suspended crystals and whole rock can be set to the steady-state compositions as described in the following section (Eqs. (29)–(32)).

and "
 ϕre þ1 ε x ¼ nD þ C 0x ε 0x −nD e− ϕ

Ms Mm

þ

ϕr e þ1 ðϕr e þ 1ÞðC 0l ε 0l −nÞ D
m Ms e Mm −e− ϕ mϕ þ ϕr e þ 1

Ms Mm

# 

=C x ;

ð17Þ

Z

l¼

t

_ s dt; M

0

ra C a þ ri C i ; ðϕre þ 1ÞD þ ð1−ϕÞre

m¼−

n¼

ð28Þ

3. Results and discussion 3.1. Trace element and isotopic evolution of the whole rock, liquid, total suspended crystals, and crystal rims

where Ms ¼

1−ð1−ϕÞD 0 Cm: ϕ

ðϕr e þ 1ÞD þ ð1−ϕÞr e ; 1−ϕ

ra C a εa þ ri C i εi ; ðϕr e þ 1ÞD þ ð1−ϕÞr e

ð18Þ

ð19Þ

ð20Þ

ð21Þ

Figs. 2–5 show examples of the evolution of the trace element and isotope ratio (87Sr/86Sr) of the liquid, total suspended crystals, magma (whole rock), and crystal rims for several values of ra, ri, re and D. Petrographic observations of natural volcanic rocks suggest that basaltic and andesitic magmas contain suspended crystals from a few wt% up to ~60 wt%, and rhyolitic magmas contain suspended crystals from a few wt% to ~20 wt% (Brophy, 1991). The weight fraction of suspended crystals ϕ is assumed here to be 0.2 for all cases. Under the assumption of constant ϕ, increments of crystals over ϕ are instantaneously separated from convecting magma to the cumulate pile throughout the crystallization process. The crystals being separated are assumed to have the same

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

D = 10

185

D = 0.1

10

1000

(a)

(b)

Liquid

1

100

0.1

10

ra = 0 ri = 0 re = 0

Total suspended crystals Magma (whole rock)

Relative concentration normalized by Cm0

0.01

1

D = 10

0.001 10

D = 0.1

0.1 1000

(c)

(d)

1

100

0.1

10

0.01

1

0.001 10

0.1 1000

ra = 0.2 ri = 0 re = 0

(e)

(f)

ra = 0.2 ri = 0.4 re = 0

100 1 10 0.1 1 0.01 10

0.1 1000

(g)

(h)

1

100

0.1

10

0.01

1

0.001

ra = 0.2 ri = 0 re = 0.6

0.1 1

0.8

0.6

0.4

0.2

0

1

0.8

0.6

0.4

0.2

0

Fm Time Fig. 2. Evolution of trace element concentrations in liquid, total suspended crystals, and magma (whole rock) for ra + ri − re b 1. Each concentration is normalized by the initial concentration of magma, Cm0. Trajectories are shown for ϕ = 0.2, Ca/Cm0 = 0.25, ra = 0.2 and different values of ri, re, and D (see Table 1 for descriptions of the parameters).

composition as those remaining in suspension. The effects of ϕ on the trace element and isotope evolution paths have been described for AFC (Nishimura, 2012). The relative trace element concentration of anatectic wall-rock melt normalized by that of the initial magma Ca/Cm0 is assumed to be 0.25 for all cases. The robustness of the assumption of a constant assimilant composition Ca is discussed in the following section. For simplicity, the recharged magma is assumed to have the same composition as the original unfractionated magma. The relative mass of magma remaining (Fm) for ra + ri − re b 1 decreases as the open-system process proceeds, whereas Fm for ra + ri − re N 1 increases as the process proceeds (an Excel spreadsheet for this calculation can be found in the Supplemen_ aþM _ i bM _s tary material S2). This is because when ra + ri − re b 1 (i.e., M

_ e), the mass of added material (assimilant + injected liquid) is smaller þM than the mass of subtracted material (gravitationally separated crystals + extracted magma), meaning that the magma mass decreases. When ra + _ aþM _ i NM _ sþM _ e), the mass of added material is greater ri − re N 1 (i.e., M than the mass of subtracted material, meaning that the magma mass increases.

For a compatible trace element (e.g., D = 10), the evolution paths of the whole rock are strongly affected by suspended crystals and are therefore nearly parallel to those of suspended crystals (left panels in Figs. 2 and 3). For an incompatible trace element (e.g., D = 0.1), the evolution paths of whole rock become close to those of liquid (right panels in Figs. 2 and 3). The results of simple fractional crystallization (ra = ri = re = 0) with 20 wt% suspended crystals are compared with the results when assimilation occurs. Given that the initial liquid had already experienced 20 wt% fractional crystallization (Eqs. (27) and (28)), Fm is related to the fraction of liquid remaining F (Ml/Mmo; as generally used in the Rayleigh equation) as F = 0.8Fm (i.e., F = (1 − ϕ)Fm). When ra = ri = re = 0, the equation for the liquid composition (Eq. (6)) becomes identical to the Rayleigh fractionation equation. Calculated liquid lines of descent (dotted lines in Fig. 2a, b) are therefore identical to those calculated from the Rayleigh fractionation equation. In this case (ra = ri = re = 0), the equation for the total suspended crystal composition (Eq. (8)) becomes identical to the ZCIS equation of Nishimura (2009). The compatible trace element concentrations of the liquid, total suspended crystals, and whole rock decrease continuously with

186

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

D = 10

D = 0.1

10

10

Liquid Total suspended crystals Magma (whole rock)

Relative concentration normalized by Cm0

(a)

1

1

0.1 10

0.1 10

(c)

(b)

ra = 0.2 ri = 2 re = 0

(d)

1

1

0.1 10

0.1 10

ra = 0.2 ri = 4 re = 0

(e)

(f)

1

ra = 0.2 ri = 4 re = 2

1

0.1 1

10

0.1 100 1

10

100

Fm Time Fig. 3. As for Fig. 2 but for ra + ri − re N 1.

decreasing Fm (Fig. 2a). When assimilation occurs and ra + ri − re b 1, the compatible trace element concentrations of liquid, total suspended crystals and whole rock reach a compositional steady state (Fig. 2c, e, g). In contrast, there is no steady state for the incompatible trace element concentrations regardless of whether assimilation occurs (right panels in Fig. 2). When ra + ri − re N 1 and assimilation occurs, steady states are attained for both the compatible and incompatible elements (Fig. 3). In the case of no recharge and no extraction (Fig. 2c, d), the evolution paths of liquid, total suspended crystals, and whole rock are identical to those obtained from the revised AFC (AIFC) model incorporating crystal zoning (Nishimura, 2012). The addition of magma recharge affects the steady state concentrations of liquid, total suspended crystals, and whole rock (Fig. 2e), whereas magma extraction does not affect the steady state concentrations (Fig. 2g, Fig. 3e, f). The evolution paths of liquid, total suspended crystals, and whole rock for ra + ri − re = 1 can be found in the Supplementary material S1. The value of the strontium distribution coefficient DSr is mainly a function of the amount of plagioclase and its anorthite (An) content. As a magma body evolves and crystallizes more sodic plagioclase, Sr becomes increasingly more compatible in the bulk solid (Blundy and Wood, 1991). The value of DSr for tholeiitic-type basaltic magma may be much less than 1, unless plagioclase is the main crystallizing phase. In contrast, the value of DSr for silicic magma may be much larger than 1 because Sr enters sodic plagioclase and/or sanidine. For example, the value of DSr for high-silica magma of the Bishop Tuff, California, USA, has been estimated as 7.04 ± 0.44 (Anderson et al., 2000). The modeling

results of 87Sr/86Sr variations for cases with DSr = 0.1 and DSr = 10 are shown in Figs. 4 and 5 to emphasize the difference between the trajectories for incompatible and compatible elements. The initial magma was assumed to have Sr = 400 ppm and 87Sr/86Sr = 0.7030; the assimilant (anatectic wall-rock melt) was assumed to have Sr = 100 ppm and 87 Sr/86Sr = 0.7200. These values are identical to those used in fig. 3 in DePaolo (1981). When DSr = 10, the 87Sr/86Sr paths of the whole rock are close to those of the total suspended crystals (left panels in Figs. 4 and 5), whereas when DSr = 0.1 the 87Sr/86Sr paths of whole rock become almost identical to those of liquid or crystal rims (right panels in Figs. 4 and 5). For ra + ri − re b 1, the 87Sr/86Sr ratios of liquid (or crystal rims), total suspended crystals, and whole rock reach steady states only when DSr is large (Fig. 4). For ra + ri − re N 1, the 87Sr/86Sr ratios of liquid (or crystal rim), total suspended crystals and whole rock reach steady states regardless of the DSr value (Fig. 5). The steady-state trace element concentrations and isotopic ratios of the liquid and total suspended crystals for ra + ri − re ≠ 1 can be expressed as C l ≈ h=i;

ð29Þ

C x ≈ ðh=iÞD;

ð30Þ

C m ≈ ðh=iÞðϕD þ 1−ϕÞ;

ð31Þ

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

D = 10

187

D = 0.1

0.7250

0.7045 Liquid / Crystal rims

(a) 0.7200

(b)

Total suspended crystals

0.7040

ra = 0.2 ri = 0 re = 0

Magma (whole rock)

0.7150 0.7035 0.7100 0.7030

0.7050 0.7000

0.7025

0.7250

0.7045

(c)

87

Sr/86Sr

0.7200

(d)

ra = 0.2 ri = 0.4 re = 0

0.7040

0.7150 0.7035 0.7100 0.7030

0.7050 0.7000

0.7025

0.7250

0.7045

(e)

0.7200

(f)

ra = 0.2 ri = 0 re = 0.6

0.7040

0.7150 0.7035 0.7100 0.7030

0.7050 0.7000

0.7025 1

0.8

0.6

0.4

0.2

1

0

0.8

0.6

0.4

0.2

0

Fm Time Fig. 4. Evolution of the isotope ratio (87Sr/86Sr) of liquid, crystal rims, total suspended crystals, and magma (whole rock) for ra + ri − re b 1. Trajectories are shown for ϕ = 0.2, Ca/Cm0 = 0.25, ra = 0.2, and different values of ri, re and D (see Table 1 for descriptions of the parameters).

and εl ≈ εx ≈ εm ≈ j=h:

ð32Þ

When ra + ri − re b 1 and D is large, then i (Eq. (11)) and k (Eq. (13)) may take large negative values. In such a case, Fm−i in Eqs. (6)–(9) and Fm−k in Eqs. (8) and (9) show a rapid reduction with decreasing Fm, and the steady state compositions are obtained as Eqs. (29)–(32). When ra + ri − re N 1, then i and k may take positive values regardless of D, and Fm increases as the process advances. In this case, Fm−i and Fm−k show a rapid reduction with increasing Fm, and the steady state compositions become identical to Eqs. (29)–(32), which can also be expressed as Cl ≈

ra C a þ ri C i ; ðr a þ r i C i −1ÞðϕD þ 1−ϕÞ þ D

ð33Þ

Cx ≈

ðr a C a þ r i C i ÞD ; ðr a þ r i C i −1ÞðϕD þ 1−ϕÞ þ D

ð34Þ

Cm ≈

ðr a C a þ r i C i ÞðϕD þ 1−ϕÞ ; ðr a þ r i C i −1ÞðϕD þ 1−ϕÞ þ D

ð35Þ

and εl ≈ εx ≈ εm ≈

ra C a εa þ ri C i εi : ra C a þ ri C i

ð36Þ

The steady state trace element compositions of liquid, total suspended crystals, and whole rock are functions of ra, ri, Ca, Ci, ϕ, and D. The steady

state isotopic ratio of the liquid, total suspended crystals, and whole rock can be expressed as a common function of ra, ri, Ca, Ci, εa, and εi. These equations show that magma extraction does not influence the steady state compositions; however, it does affect the Fm values at the point that the steady state compositions are achieved (Figs. 2–5). The steady state trace element and isotopic composition for ra + ri − re = 1 can be found in Supplementary material S1. The bulk isotope ratio of crystals shows an increase after assimilation, but the effect is buffered compared with the dramatic change in the isotope ratio of the liquid, because there is already a large volume of crystal with a relatively low bulk isotope ratio (Davidson et al., 2007). When the magma reaches a steady state composition, the texture of the magma or rock experiences the three stages shown schematically in Fig. 6. Before the liquid composition reaches a steady state, the crystal rim composition changes continuously as the magma evolves (Stage 1 in Fig. 6). After the liquid composition reaches a steady state, growth with a steady state composition occurs at crystal margins (Stage 2 in Fig. 6). The zoned crystals are continuously settled out from the magma body and are reduced in number. After the total suspended crystals reach a steady state composition, an isotopically homogeneous magma is obtained (Stage 3 in Fig. 6). Trace element and isotopic variations recorded from core to rim of a single crystal grain reflect progressive changes in the liquid composition in response to open-system magmatic processes (e.g., Davidson et al., 1998, 2007). Fig. 7 shows the 87Sr/86Sr ratios and Sr concentrations from core to rim of a plagioclase phenocryst calculated for Dsr (bulk partition coefficient) = 10, DPlSr (plagioclase/melt partition coefficient) =

188

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

D = 10 0.7035

D = 0.1 (a)

0.7035

0.7034

0.7034

0.7033

0.7033

0.7032

0.7032

0.7031

0.7031

0.7030

0.7030

0.7029

0.7029

(b)

ra = 0.2 ri = 2 re = 0

Liquid / Crystal rims

0.7035

Sr/86Sr

Magma (whole rock)

0.7035

(c)

0.7034

87

Total suspended crystals

0.7033

0.7033

0.7032

0.7032

0.7031

0.7031

0.7030

0.7030

0.7029

0.7029

0.7035

(e)

0.7034

ra = 0.2 ri = 4 re = 0

0.7035

(f)

0.7034

0.7033

0.7033

0.7032

0.7032

0.7031

0.7031

0.7030

(d)

0.7034

ra = 0.2 ri = 4 re = 2

0.7030

0.7029 1

10

0.7029 100 1

10

100

Fm Time Fig. 5. As for Fig. 4 but for ra + ri − re N 1.

10, ra = 0.2, and several values of ri (the values of other compositional parameters are the same as in Figs. 4 and 5). During crystal growth, the rim composition evolves towards a steady state composition. The chemical evolution trend recorded within the plagioclase phenocryst varies markedly in response to the ri value (Fig. 7). This result suggests that the trace element and isotopic variations preserved within a single crystal provide a key to understanding the values of physical parameters such as the magma influx rate and the assimilation rate. Chemical evolution trends observed in natural plagioclase phenocrysts, however, often show more complex patterns (e.g., Tepley et al., 2000), suggesting changes in parameters such as ri and D. 3.2. Treatment of anatectic wall-rock melt DePaolo (1981) assumed a constant composition for anatectic wallrock melt to analytically solve the differential equations for AFC. Inversions of the DePaolo (1981) equations have been performed to investigate the compositions of the parental magma and contaminant, amount of assimilation, and relative rates of assimilation and crystal fractionation (e.g., Powell, 1984; Mantovani and Hawkesworth, 1990; Young et al., 1992; Aitcheson and Forrest, 1994; Roberts and Clemens, 1995; Caffe et al., 2002). The present model also assumes a constant wallrock melt composition and this assumption yields the compositional steady states (Figs. 2–7). In contrast, Spera and Bohrson's energy-constrained open-system magma chamber models (Bohrson and Spera, 2001, 2003, 2007; Spera and Bohrson, 2001, 2002, 2004) quantitatively addressed how changes

in the composition of wall-rock melt affect magma trace element and isotopic compositions. In their ECRAFC model, for an incompatible element in the wall-rock for which the bulk D is ≪1, and as the wall-rock is partially melted, the availability (mass) of that element is depleted in the first few steps of assimilation by fractional melting; the degree to which this happens is a function of the concentration and the bulk distribution coefficient D in the wall-rock (Bohrson and Spera, 2001, 2003, 2007; Spera and Bohrson, 2001, 2002, 2004). This causes a reduction in the rate of isotopic change as AFC proceeds because the mass of the trace element from the wall-rock quickly drops to almost 0. Hence the mass balance for the isotope ratio is heavily weighted to that of the magma chamber, meaning it acts almost like a closed system in which the radiogenic isotope ratio is constant. These features are opposite to those in Fig. 4, in which the ‘compatible’ elements show a flattening pattern. Whether the composition of assimilated wall-rock melt is constant or variable depends on the assimilation mechanisms. If the melt is continuously removed from the wall-rock to magma, the melt composition changes with time. For example, the Magma Chamber Simulator (MCS) formulation of Bohrson et al. (2014) assumes that when the melt fraction in the wall-rock exceeds a critical melt volume fraction (0.04–0.12), partial-melt is instantaneously removed or added to the magma by percolative flow. This model tracks the melting conditions of the wall-rock and describes the effects of changing partial-melt composition on the trace element and isotopic composition of the magma. In contrast, if a highly molten wall-rock (magma) with a constant temperature is mixed into the main magma body, the wall-rock melt composition will be near-constant. In the heat and mass transfer

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

189

Isotope ratio

Liquid / Crystal rims Total suspended crystals Magma (whole rock)

Time

Isotope ratio

Stage 1

Stage 2

Stage 3

High

Low

Crystal zoning

Growth of crystal margins with a steady state composition

Isotopically homogeneous magma

Fig. 6. Schematic illustration of the relationship between the isotopic heterogeneity of magma and the isotopic evolution paths of liquid, crystal rims, total suspended crystals, and whole rock. The composition of wall-rock melt is assumed to be constant. Crystals are homogeneously distributed by vigorous convection and are partially separated from the magma body to the cumulate pile at the base of the magma chamber due to gravity. The crystals being separated are assumed to have the same composition as those remaining in suspension. Crystal nucleation and growth occur continuously. During stage 1, crystal rims show a continuous change in isotope ratio. During stage 2, growth zones with a steady-state isotope ratio develop at crystal margins. The zoned crystals with early formed cores are continuously settled out from the magma body and gradually decrease in abundance. During stage 3, isotopically homogeneous magma develops. See text for details.

model of Huppert and Sparks (1988), when the partial melt in a roofrock exceeds a critical melt fraction (30%–50%; e.g., Shaw, 1980; Marsh, 1981), the connectedness of the crystalline framework is destroyed and the rock behaves as a liquid with suspended crystals (i.e., a magma). The temperature at this critical melt fraction is called the effective fusion temperature (EFT) (Huppert and Sparks, 1988). In 0.706

ri = 0

Sr/86Sr

0.705

Steady state

Rim

the case of melting at the roof of a magma chamber, the heat and mass are efficiently conveyed by vigorous convection above the EFT. The anatectic roof-rock magma will be mixed into the main magma body unless the density of the anatectic magma at a critical melt fraction at the EFT is lower than that of the main magma. The boundary between the partially molten roof-rock and magma chamber has a fixed EFT and moves upward as melting proceeds (Huppert and Sparks, 1988; Koyaguchi and Kaneko, 1999, 2000; Nishimura et al., 2005). In this case, the partial melt, and xenocrysts and/or xenoliths, from the roofrock at the EFT can be mixed into the resident magma. The present mass balance model can be applied in this case, assuming negligible chemical interaction between the magma and xenocrysts and/or xenoliths.

87

3.3. Application of the model to successive recharge and eruption events 0.3

0.704

0.7

Core

1 1.5

2

2.5

0.703 0

100

200

300

400

500

600

700

Sr (ppm) Fig. 7. 87Sr/86Sr ratios and Sr concentrations from core to rim of a plagioclase phenocryst calculated for Dsr (bulk partition coefficient) = 10, DPlSr (plagioclase/melt partition coefficient) = 10, ra = 0.2, and several values of ri (other compositional parameter values are as used in Figs. 4 and 5). During crystal growth, the rim composition evolves towards a steady state composition.

Successive recharge and eruption events (Fig. 8) can be modeled by the repeated use of the analytical solutions (Eqs. (6)–(9) and (14)– (17)) that are applicable for each crystallization interval of constant ra, ri, re, ϕ and D. At the point of changes in these parameters, the initial compositions of liquid and total suspended crystals are reset and the calculation is restarted. This calculation technique is principally applicable to any complex changes in constant parameter values (Shaw, 2006; Nishimura, 2009). The effect of the partial dissolution of crystals on the mass balance is neglected for simplification. The simplified schematic diagrams in Fig. 8a show successive recharge and eruption (extraction) events. Fresh magma with the original

190

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

(a)

AFC (assimilation and fractional crystallization)

AFC + Recharge

AFC + (Rapid) recharge + Extraction

(b)

0.720

Liquid / Crystal rims Total suspended crystals Magma (whole rock)

0.710

87

Sr/86Sr

0.715

0.705 0.700

Fm

1.0

ra = 0.6 ri = 0 re = 0

0.8

(c) ra = 0.6 ri = 0.5 re = 0

ra = 0.6 ri = 5 re = 5

0.6 0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

Mc/Mm0 Fig. 8. Evolution of the isotope ratio (87Sr/86Sr) of liquid, crystal rims, total suspended crystals, and magma (whole rock) during successive recharge and eruption events. Changes in the relative mass of magma remaining (Fm) are also shown. Mc/Mm0 is the total crystallized mass normalized by the initial magma mass. See text for details.

magma composition enters the magma chamber undergoing AFC and triggers an eruption. The 87Sr/86Sr ratios of liquid, crystal rims, total suspended crystals, and whole rock evolve towards the steady state ratios for each stage (Fig. 8b). The 87Sr/86Sr ratio of magma also depends on ra, ϕ and D, which here are assumed to be constant to focus on the variations in recharge and eruption rates. The values of ra, ri and re are shown in Fig. 8c and the values of other parameters are the same as in Fig. 7. The total crystallized mass normalized by the initial magma mass Mc/Mm0 (plotted as the horizontal axis in Fig. 8b, c) irreversibly increases as the open-system process advances. The relative mass of magma remaining (Fm) for ra + ri − re b 1 decreases with increasing Mc/Mm0, whereas Fm for ra + ri − re N 1 increases with increasing Mc/Mm0 (Fig. 8c). The 87Sr/86Sr ratio of the total suspended crystals has a small rate of change compared with the dramatic change in the 87Sr/86Sr ratio of liquid (Fig. 8b), because there is already a large volume of zoned crystal (Davidson et al., 2007). The 87Sr/86Sr ratio of the liquid (or crystal rims) rapidly increases towards a steady state composition during the AFC stage and is then followed by two sharp decreases to the steady state compositions determined by changes in the recharge rate under constant ra (Eq. (36)). The liquid or crystal rims maintain a higher 87Sr/86Sr ratio than the total suspended crystals before the onset of recharge. However, this relationship is promptly reversed by the recharge event (Fig. 8b). This result suggests that the isotopic relationship between liquid and total suspended crystals plays a role in determining the presence or absence of recharge events. 3.4. Reproduction of natural zoning patterns Davidson et al. (2007) presented typical 87Sr/86Sr zoning patterns of feldspars from contrasting magma systems: El Chichón (McGee and

Tilling, 1983; Tilling et al., 1984; McGee et al., 1987; Espíndola et al., 2000; Tepley et al., 2000; Andrews et al., 2008; Tilling, 2009; Arce et al., 2014, 2015), Ngauruhoe (Steiner, 1958; Hobden, 1997; Hobden et al., 1999, 2002), and Taylor Creek (Duffield and Dalrymple, 1990; Duffield and du Bray, 1990; Duffield and Ruiz, 1992, 1994; Wittke et al., 1996; Knesel et al., 1999) volcanoes (Fig. 9a–c). The magmatic processes for these three cases were constrained by combining the zoning data with textural observations and other geochemical data (e.g., Knesel et al., 1999; Tepley et al., 2000; Davidson et al., 2001, 2007). At El Chichón Volcano (Mexico), recharge events are recorded as the dissolution surfaces and/or significant compositional breaks within plagioclase phenocrysts from trachyandesite lavas of the 1982 and pre-1982 eruptions (Fig. 9a). Tepley et al. (2000) suggested that Sr-isotope variations across the plagioclase phenocrysts most probably record crystal growth in a magma chamber into which magma with a lower 87Sr/86Sr ratio and a higher Sr concentration was periodically added. Crystal residence times (~500 years) estimated from diffusion modeling indicate that the magma recharge cycles are shorter than the eruptive occurrence interval (Andrews et al., 2008). At Ngauruhoe Volcano (New Zealand), Sr-isotope data across plagioclase phenocrysts from basaltic andesite lavas reveal crystal cores with lower 87Sr/86Sr and crystal rims with higher 87Sr/86Sr (Fig. 9b). This observation is consistent with the later growth of crystals in more assimilated magma. Davidson et al. (2007) proposed that dissolution surfaces in crystals indicate transfer between magmas rather than continuous assimilation. Finally, at Taylor Creek (New Mexico, USA), continuous core-to-rim variations in 87Sr/86Sr are observed in sanidine phenocrysts from a rhyolite lava dome (Fig. 9c). The phenocrysts contain no distinct dissolution surfaces. Zoning data for several sanidine phenocrysts of different size show the same overall pattern, with low 87Sr/86Sr in the core, increasing to a maximum value in the mantle zone and then decreasing towards the rim (Knesel et al., 1999; Davidson et al., 2007). Knesel et al. (1999) suggested that (1) the 87Sr/86Sr profiles cannot be produced by diffusion; (2) high-87Sr/86Sr crust was effectively assimilated at an early stage of magmatic evolution, followed by a renewed influx of low-87Sr/86Sr, low-Sr rhyolite; and (3) crystal nucleation and growth occurred throughout much, if not all, of the open-system evolution of the magma. Davidson et al. (2007) concluded that the continuous variation in 87Sr/86Sr (Fig. 9c) records growth in a magma that gradually changed in composition rather than punctuated events. Core-to-rim variations in 87Sr/86Sr and Fm are calculated using the present mass balance equations (Fig. 9d–i). Based on the above descriptions of previous studies, periodic magma recharge (Model 1 in Fig. 10), crystal transfer from lower to upper magma chambers (Model 2 in Fig. 10), and continuous change in the magma recharge rate (Model 3 in Fig. 10) are assumed for El Chichón, Ngauruhoe, and Taylor Creek, respectively. The model calculations are performed by the repeated use of the mass balance equations as described in the previous section. Bulk and feldspar/melt distribution coefficients for Sr are assumed to be constant for simplification (The detailed values of changing distribution coefficients and crystal fraction should be incorporated in future models). The parameter values used for the calculation are listed in Table 2. The phenocryst contents of El Chichón trachyandesite, Ngauruhoe basaltic andesite, and Taylor Creek rhyolite are ~20 vol% (Arce et al., 2015), 25–37 vol% (Hobden, 1997), and 12–38 vol% (Duffield and Ruiz, 1992), respectively. In the model calculations, the weight fraction of suspended crystals is assumed to be 0.2 (20 wt%) in all cases. This paper employs the cubic root of Mc/Mm0 (total crystallized mass normalized by initial magma mass) as the horizontal axis of crystal zoning profile (Fig. 9d–i) because crystal length is proportional to the cubic root of the mass ratio. Assuming that crystal nucleation and growth occur continuously in the open-system magma chamber (Cashman and Marsh, 1988; Marsh, 1988a), zoning profiles were calculated for feldspars that nucleate at Mc/Mm0 = 0. Later nucleation and growth result in

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

El Chichón

Ngauruhoe

(a)

191

Taylor Creek

(c)

87

Sr/86Sr

(b)

Core 0.7053

Rim

(d)

0.7051

Model 1

Sr/86Sr

Rim

Rim

(e)

Model 2

0.717

(f)

Model 3

0.716 0.715

0.7047

0.714

0.706

Upper magma

0.7045

0.713

0.7055

0.7043

0.712 0.711

0.7041

0.705

0.71

Lower magma

0.7039

0.709

0.7045

0.7037 0 1.2

0.2

0.4

0.8

0

1

0

0

1

0.6

40

(g)

1.1

20

0.9

Fm

Core 0.718

0.7065

0.7049

87

Core 0.707

35

ri = 0

20

1

2

3

2

(h)

0

0.6

20

0.4

1.2

0.2

0

33

Mc / M

0.8

0 m

(i)

1

6

1.4

ri = 1.5 Upper magma

10 5

0.6

0.8

5

Lower magma

15

0.3 0.4

0.7

1.6

20

0.7

0.2

0.6

1.8

25

0

0.5

30

0.8

0.5

0.708 0.4

4

1.3 1.5

1

ri = 0.3 0

1

33

2

Mc / M

0 m

3

0.8 0.4

3

ri = 1.2

0.5

0.6

33

2

0.7

Mc / M

0.8

0 m

Fig. 9. Reproduction of natural zoning patterns. (a–c) Typical 87Sr/86Sr zoning patterns of feldspars from different magma systems, taken from Davidson et al. (2007). (a) Plagioclase from El Chichón volcano. 87Sr/86Sr ranges from 0.7040 to 0.7052 and the typical distance from core to rim is 2.5–3.5 mm (Tepley et al., 2000). (b) Plagioclase from Ngauruhoe volcano. 87Sr/86Sr ranges from 0.7050 to 0.7060 and the typical distance from core to rim is 0.6–1.2 mm (Davidson et al., 2007). (c) Sanidine from Taylor Creek volcano. 87Sr/86Sr ranges from 0.7050 to 0.7180 and the typical distance from core to rim is 2.5–3.5 mm (Knesel et al., 1999). (d–f) Calculated zoning patterns. (g–i) Calculated ratio of magma mass to the initial magma mass (Fm). See text for details.

smaller crystals with evolved core compositions as observed in Taylor Creek Rhyolite (Knesel et al., 1999). All three cases yield similar core-to-rim variations in 87Sr/86Sr to those found in natural crystals (Fig. 9d–f). Calculated results for El Chichón reveal that recharge events cause sharp compositional gaps and increasing Fm values (Fig. 9d, g). In this model, however, Fm decreases markedly during periods of no recharge and reaches 0.3 at crystal rims (Fig. 9g). The constant (homogeneous) core compositions of Ngauruhoe (Fig. 9b) suggest that the plagioclase crystallized from a magma chamber in a compositional steady state (or a closed-system state). The calculated compositional profile (Fig. 9e) shows crystal

transfer from a lower magma chamber in a compositional steady state to an upper magma chamber. The inflation rate of the lower magma chamber is greater than that of upper magma chamber (Fig. 9h). This modeling can be applied to more complex polybaric crystallization along vertically extended plumbing systems (e.g., Hamada et al., 2014; Ubide and Kamber, 2018). The calculated results for Taylor Creek show that 87Sr/86Sr in crystallizing rims increases at constant ri (1.2) and then decreases with increasing ri (Fig. 9f, i). The mass of the magma increases with increasing ri. This zoning pattern has also been interpreted as a record of the compositional change that results from the assimilation of wall-rock melt derived from biotite- and feldspar-

192

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

Model 1

Time

AFC AFC + Recharge (assimilation and fractional crystallization)

AFC

Model 2

AFC + Recharge

Model 3

Steady state composition Change in recharge rate Fig. 10. Schematic illustrations of three open-system magma chamber models: periodic magma recharge (Model 1), crystal transfer from lower to upper magma chambers (Model 2), and continuous change in the magma recharge rate (Model 3). Assimilation and fractional crystallization (AFC) processes are assumed in all models.

bearing roof-rocks; the wall-rock melt initially comes mostly from high-87Sr/86Sr biotite, and then comes mostly from low-87Sr/86Sr feldspar (Duffield and Ruiz, 1994). It should be noted that other interpretations are possible for the all three cases, and each model might need revision based on further petrological and geochemical studies. Core-to-rim variations in 87Sr/86Sr for El Chichón and Ngauruhoe can also be explained by the Model 3 (Fig. 10). A stepwise increase (decrease) in ri value results in a stepwise decrease (increase) in 87Sr/86Sr (Fig. 11). For both cases, Fm increases as the process advances. The dissolution surfaces observed in the plagioclase from El Chichón can also be explained by an increase in the recharge rate when using Model 3 (the relationships among Model 1, Model 3, dissolution surfaces and the resolution of analysis are described in the following section). Results for Model 1 and Model 3 can clearly be discriminated from each other Table 2 Parameter values used in model calculations. El Ngauruhoe Chichón

Taylor Creek

El Chichón

Ngauruhoe

Model 1 Model 2

Model 3

Model 3

Model 3

Lower Upper magma magma ra ri re DSr DfeldsparSr ϕ Cm0 Sr εm0 (87Sr/86Sr) CSr a εa (87Sr/86Sr) CSri εi (87Sr/86Sr)

0.2 0, 20 0 0.5 1.5 0.2 850 0.7050

1 1.5

0.95 1.2, 1.3–6.7 0 7a 11a 0.2 10 0.7050d

0.9 0.18, 0.3, 0.45, 0.8 0 0.5 1.5 0.2 850 0.7050

1 1.5, 0.3

200 0.7080c

10 0.8000d

300 0.7070b

200 0.7080c

400 432.4⁎ 0.7040 0.7050⁎

10 0.7040

1500 0.7030

400 0.7040

0.3 0 0.5 0.5 1.5 1.5 0.2 0.2 389 432.4⁎ 0.7050 0.7050⁎

300 200 0.7070b 0.7080c 1500 0.7030

1 0.3

0 0.5 1.5 0.2 389 0.7050

References: a. Anderson et al. (2000), b. Tepley et al. (2000), c. Snelling (2003), d. Duffield and Ruiz (1994). ⁎ Steady state composition of liquid in the lower magma chamber.

with reference to a diagram showing core-to-rim variations in Sr versus Sr/86Sr (Fig. 12). Model 1 causes a stepwise core-to-rim trend and highly Sr-rich rim that are not observed in the El Chichón plagioclase phenocrysts. When using Model 3, changes in the steady state composition produce a smoothly curving core-to-rim trend similar to that in an El Chichón plagioclase phenocryst (Fig. 12). 87

3.5. Diffusional relaxation of a compositional step The core-to-rim variation of a rapidly diffusing element can change during storage at magmatic temperatures, even if the compositional breaks are initially recorded. The degree of diffusive relaxation of Sr and Mg profiles in plagioclase phenocrysts can be estimated using a simple one-dimensional diffusion equation in infinite space. The diffusion coefficients of both elements are largely dependent on the anorthite (An) content, as well as temperature (e.g., Cherniak and Watson, 1992, 1994; Costa et al., 2003). Initial compositional steps of Sr (from 500 to 400 ppm) and Mg (from 200 to 100 ppm) were assumed in a plagioclase phenocryst with homogeneous An contents of An93 and An23 (Fig. 13). Although natural plagioclase phenocrysts generally have An zoning profiles, the simple assumption of constant An content was used here to obtain an analytical solution of the diffusion equation at constant diffusivity (details of the diffusion modeling can be found in Supplementary material S1). Recent studies of diffusion in feldspars suggest that many feldspars from various types of volcanic rock have magma residence times ranging from a few days to a few decades (e.g., Costa et al., 2003, 2010; Saunders et al., 2010; Druitt et al., 2012; Moore et al., 2014; Till et al., 2015; Singer et al., 2016). The plagioclase phenocrysts are assumed to be kept at a magmatic temperature of 1000 °C for 0.1 and 10 yr. In plagioclase phenocrysts larger than 0.5 mm in size, the concentration step in Sr is recognizable after the diffusion time for all cases, although diffusion slightly blurs the initial compositional boundary (Fig. 13). In contrast, the rapid diffusion of Mg produces moderately sloping, continuous Mg variations in the plagioclase with An23 (oligoclase). This feature becomes more pronounced at high temperature because of the increase in the diffusion coefficient. Furthermore, crystals from rhyolitic supereruptions often show long diffusion timescales of 100–1000 yr (Wark et al., 2007; Saunders et al.,

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

El Chichón 0.7053

Ngauruhoe 0.7070

(a)

0.7051

Sr/86Sr

(b)

Model 3

Model 3

0.7065

0.7049

87

193

0.7047 0.7060

0.7045 0.7043

0.7055

0.7041 0.7039

0.7050

0.7037 0.7045

0.7035 0 350

1

2

3

4

5

6

7

8

80

(c)

1

2

3

4

5

6

5

6

(d)

70

300

60

250

Fm

0

9

ri = 0.8

200

50

ri = 0.3

40

150

30

100

20

0.45 50

0.18

1.5

10

0.3

0

0 0

1

2

3

3

4

5

6

7

8

M c / M m0

9

0

1

2

3

3

4

M c / M m0

Fig. 11. 87Sr/86Sr zoning patterns and Fm values for El Chichón and Ngauruhoe calculated using Model 3.

2010; Druitt et al., 2012; Allan et al., 2013). In those cases, Sr diffusion profiles might develop in plagioclase. Therefore, when applying the present mass balance model to a rapidly diffusing element, detailed constraints on the initial core-to-rim profiles (e.g., Costa et al., 2010; Druitt et al., 2012; Moore et al., 2014; Till et al., 2015) are required to compare with the model results. 3.6. Resolution dependence If the time interval of periodic magma recharge is too short (e.g., a few days) to be recorded in crystal growth zoning, the core-to-rim

Fig. 12. Comparison of the core-to-rim chemical variations for Model 1 (broken line) and Model 3 (solid line) on an 87Sr/86Sr ratio versus Sr diagram. Solid squares indicate the microanalysis data across a single plagioclase crystal from El Chichón (sample 74–26 Pl; Tepley et al., 2000).

chemical variations will show apparently continuous recharge. The solid line in Fig. 14 shows a schematic illustration of mass changes due to periodic recharge events. Dissolution surfaces are formed soon after a large magma influx, as observed at El Chichón (Tepley et al., 2000; Davidson et al., 2007; Andrews et al., 2008). During gradual mixing, however, crystal growth takes place without dissolution, as observed at Taylor Creek (Knesel et al., 1999: Davidson et al., 2007). The typical sampling scale for microanalysis can be estimated as ~50 μm from an examination of the distribution of micromill pits (Davidson et al., 2007). Core-to-rim isotopic profiles of crystals, therefore, show ~50 μm moving averages. Using typical plagioclase growth rates of 10−8 to 10−10 cm s−1 (Cashman, 1993; Larsen, 2005; Armienti, 2008; Pupier et al., 2008; Ruprecht et al., 2008), such a crystal would take between ~2 months and 160 years to grow to ~50 μm in size. When nonrecharge periods are sufficiently longer than the time taken for ~50 μm of crystal growth, the growth zones that are formed during nonrecharge periods can be recognized in core-to-rim profiles. In this case, the core-to-rim chemical variations show evidence of periodic recharge (Model 1 in Fig. 9). When non-recharge periods are nearly equal to or shorter than the time taken for ~50 μm of crystal growth, apparently continuous recharge is recorded because of the scale of the moving average of long-term recharged mass (dashed line in Fig. 14). In this case, the core-to-rim chemical variations could reflect changes in the recharge rate (Model 3 in Fig. 9). It should be noted that other interpretations might exist for this type of core-to-rim variation. If ri is estimated from the modeling of core-to-rim chemical varia_ s , where M _ i can be obtained as M _ i ¼ ri M _s tions, the recharge rate M can be calculated from the basal area of the magma chamber × crystal fraction × crystal settling velocity. The settling velocity can be estimated from the melt viscosity, the crystal–melt density contrast, the crystal fraction, and crystal size (e.g., Richardson and Zaki, 1954; Gray and

194

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

Fig. 13. Initial patterns of Sr and Mg zoning (solid lines) and calculated distributions of Sr (dotted lines) and Mg (dashed lines) after diffusion at 1000 °C for various time periods. The diffusion coefficient of Mg in plagioclase is D = 2.92 × 10−4.1XAn−3.1exp[−266,000 / (RT)] (Costa et al., 2003). The diffusion coefficients of Sr in anorthite (An93) and in oligoclase (An23) are D = 3.85 × 10−6exp[−330,000 / (RT)] (Cherniak and Watson, 1992) and D = 8.43 × 10−7exp[−273,000 / (RT)] (Cherniak and Watson, 1994), respectively, where D is the element diffusion coefficient in m2 s−1, XAn is the mole fraction of anorthite, T is the temperature in kelvin and R is the ideal gas constant. Sr and Mg profiles in anorthite (An93) are shown after 0.1 years (a) and 10 years (b), and in oligoclase (An23) after 0.1 years (c) and 10 years (d).

Integrated mass of recharged magma

Crain, 1969; Fujii, 1974; Marsh, 1988b; Martin and Nokes, 1988; Nishimura, 2006; Verhoeven and Schmalzl, 2009). Estimates of a past recharge rates provide a key to the prediction of eruptions, especially when the amount of magma chamber expansion necessary for eruption is poorly constrained (e.g., the Yellowstone caldera-forming eruption, USA). If the mean crystal residence time (τ) is obtained from the diffusion analysis of fast-diffusing elements, the mean crystal growth rate (G) can

Partial dissolution of crystals

Partial dissolution of crystals

be calculated from the crystal size distribution (CSD). In an open-system magma chamber with a steady state of magma replenishment and withdrawal, continuous nucleation and growth result in an exponential crystal-size distribution, as is common in many lavas (e.g., Cashman and Marsh, 1988; Marsh, 1988a, 1998). The slope on a plot of the natural log of population density versus crystal size is −1 / (Gτ), from which the growth rate G can be found if the crystal residence time τ is independently known. However, the residence times estimated by diffusion chronometry have large uncertainties of ~±0.4 log unit (Kanno and Ishibashi, 2017), which are then propagated through to uncertainties on the crystal growth rate. Nevertheless, the growth rate estimation can provide important insights into the timescale of recharge events from crystal zoning. For a magma system with periodic magma recharge, the combination of the estimates of the past recharge rate and period could have implications for monitoring future activity.

4. Conclusion

Actual mass change Decipherable mass change when non-recharge period = the crystallization time to grow to 50µm in size Time Fig. 14. Schematic illustration of the relationship between actual and decipherable mass changes during periodic magma-recharge events. See text for details.

The present study formulated the mass balance equations describing trace element and isotopic evolution during open-system magma chamber processes associated with the development of crystal zoning. These equations can be applied to successive events including episodic magma recharge and extraction. Several 87Sr/86Sr zoning patterns of feldspars from contrasting magma systems were well reproduced by using these equations. The magma recharge rate has a strong effect on crystal core-to-rim profiles of trace element concentrations and isotopic ratios. When non-recharge periods are longer than the time required for the crystal to grow by ~50 μm (i.e., the sampling scale of crystals during microanalysis), core-to-rim profiles may show episodic recharge events. In contrast, when non-recharge periods are nearly equal to or shorter than the time required for the crystal to grow by ~50 μm, the core-to-rim profiles may show apparently continuous recharge. The core-to-rim profiles can also be affected by other processes such as assimilation, eruption, and element diffusion in crystals. Nevertheless, records of trace element and isotopic evolution preserved within crystals

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

might provide a key to estimating the recharge rate leading up to past eruptions. Acknowledgements W. Bohrson and T. Ubide are thanked for their constructive reviews. Discussions with T. Kawamoto, T. Shibata, M. Yoshikawa, T. Kuritani, M. Hamada, and K. Takemura contributed substantially to the content of this paper. This work was supported by a Toyo University Individual Research Grant (2017–2018). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.jvolgeores.2019.02.012. References Aitcheson, S.J., Forrest, A.H., 1994. Quantification of crustal contamination in open magmatic systems. J. Petrol. 35, 461–488. Albarède, F., 1995. Introduction to Geochemical Modeling. 543. Cambridge University Press, Cambridge. Allan, A.S., Morgan, D.J., Wilson, C.J., Millet, M.A., 2013. From mush to eruption in centuries: assembly of the super-sized Oruanui magma body. Contrib. Mineral. Petrol. 166, 143–164. Anderson, A.T., Davis, A.M., Lu, F., 2000. Evolution of Bishop Tuff rhyolitic magma based on melt and magnetite inclusions and zoned phenocrysts. J. Petrol. 41, 449–473. Andrews, B.J., Gardner, J.E., Housh, T.B., 2008. Repeated recharge, assimilation, and hybridization in magmas erupted from El Chichón as recorded by plagioclase and amphibole phenocrysts. J. Volcanol. Geotherm. Res. 175, 415–426. Arce, J.L., Walker, J., Keppie, J.D., 2014. Petrology of two contrasting Mexican volcanoes, the Chiapanecan (El Chichón) and Central American (Tacaná) volcanic belts: the result of rift- versus subduction-related volcanism. Int. Geol. Rev. 56, 501–524. Arce, J.L., Walker, J., Keppie, J.D., 2015. Petrology and geochemistry of El Chichón and Tacaná: two active, yet contrasting Mexican volcanoes. In: Scolamacchia, T., Macías, J. (Eds.), Active Volcanoes of Chiapas (Mexico): El Chichón and Tacaná. Active Volcanoes of the World. Springer, Berlin, Heidelberg, pp. 25–43. Armienti, P., 2008. Decryption of igneous rock textures: crystal size distribution tools. Rev. Mineral. Geochem. 69, 623–650. Arth, J.G., Zen, E., Sellers, G., Hammarstrom, J., 1986. High initial Sr isotopic ratios and evidence for magma mixing in the Pioneer batholith of southwest Montana. J. Geol. 94, 419–430. Blundy, J.D., Wood, B.J., 1991. Crystal-chemical controls on the partitioning of Sr and Ba between plagioclase feldspar, silicate melts, and hydrothermal solutions. Geochim. Cosmochim. Acta 55, 193–209. Bohrson, W.A., Spera, F.J., 2001. Energy-constrained open-system magmatic processes II: application of energy-constrained assimilation–fractional crystallization (EC-AFC) model to magmatic systems. J. Petrol. 42, 1019–1041. Bohrson, W.A., Spera, F.J., 2003. Energy-constrained open-system magmatic processes IV: geochemical, thermal and mass consequences of energy-constrained recharge, assimilation and fractional crystallization (EC-RAFC). Geochem. Geophys. Geosyst. 4. https://doi.org/10.1029/2002GC00316. Bohrson, W.A., Spera, F.J., 2007. Energy-constrained recharge, assimilation and fractional crystallization (EC-RAχFC): a Visual Basic computer code for calculating trace element and isotope variations of open-system magmatic systems. Geochem. Geophys. Geosyst. 8, Q11003. https://doi.org/10.1029/2007GC001781. Bohrson, W.A., Spera, F.J., Ghiorso, M.S., Brown, G.A., Creamer, J.B., Mayfield, A., 2014. Thermodynamic model for energy-constrained open-system evolution of crustal magma bodies undergoing simultaneous recharge, assimilation and crystallization: the magma chamber simulator. J. Petrol. 55, 1685–1717. Brophy, J.G., 1991. Composition gaps, critical crystallinity, and fractional crystallization in orogenic (calc-alkaline) magmatic systems. Contrib. Mineral. Petrol. 109, 173–182. Caffe, P.J., Trumbull, R.B., Coira, B.L., Romer, R.L., 2002. Petrogenesis of early Neogene magmatism in the northern Puna; implications for magma genesis and crustal processes in the central Andean Plateau. J. Petrol. 43, 907–942. Cashman, K.V., 1993. Relationship between plagioclase crystallization and cooling rate in basaltic melts. Contrib. Mineral. Petrol. 113, 126–142. Cashman, K.V., Marsh, B.D., 1988. Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization: 2. Makaopuhi lava lake. Contrib. Mineral. Petrol. 99, 292–305. Charlier, B.L.A., Ginibre, C., Morgan, D., Nowell, G.M., Pearson, D.G., Davidson, J.P., Ottley, C.J., 2006. Methods for the microsampling and high-precision analysis of strontium and rubidium isotopes at single crystal scale for petrological and geochronological applications. Chem. Geol. 232, 114–133. Charlier, B.L.A., Wilson, C.J.N., Davidson, J.P., 2008. Rapid open-system assembly of a large silicic magma body: time-resolved evidence from cored plagioclase crystals in the Oruanui eruption deposits, New Zealand. Contrib. Mineral. Petrol. 156, 799–813. Cherniak, D.J., Watson, E.B., 1992. A study of strontium diffusion in K-feldspar, Na-K feldspar and anorthite using ion implantation and Rutherford backscattering spectroscopy. Earth Planet. Sci. Lett. 113, 411–425.

195

Cherniak, D.J., Watson, E.B., 1994. A study of strontium diffusion in plagioclase using Rutherford backscattering spectroscopy. Geochim. Cosmochim. Acta 58, 5179–5190. Cortini, M., van Calsteren, P.W.C., 1985. Lead isotope differences between whole-rock and phenocrysts in recent lavas from southern Italy. Nature 314, 343–345. Costa, F., Chakraborty, S., Dohmen, R., 2003. Diffusion coupling between trace and major elements and a model for calculation of magma residence times using plagioclase. Geochim. Cosmochim. Acta 67, 2189–2200. Costa, F., Coogan, L.A., Chakraborty, S., 2010. The time scales of magma mixing and mingling involving primitive melts and melt–mush interaction at mid-ocean ridges. Contrib. Mineral. Petrol. 159, 371–387. Davidson, J.P., Tepley III, F.J., 1997. Recharge in volcanic systems: evidence from isotope profiles of phenocrysts. Science 275, 826–829. Davidson, J.P., Wilson, I.R., 1989. Evolution of an alkali basalt–trachyte suite from Jebel Marra volcano, Sudan, through assimilation and fractional crystallization. Earth Planet. Sci. Lett. 95, 141–160. Davidson, J.P., Tepley III, F.J., Knesel, K.M., 1998. Isotopic fingerprinting may provide insights into evolution of magmatic systems. EOS Trans. Am. Geophys. Union 79 185 (189), 193. Davidson, J.P., Tepley III, F.J., Palacz, Z., Meffan-Main, S., 2001. Magma recharge, contamination and residence times revealed by in situ laser ablation isotopic analysis of feldspar in volcanic rocks. Earth Planet. Sci. Lett. 184, 427–442. Davidson, J.P., Morgan, D.J., Charlier, B.L.A., Harlou, R., Hora, J.M., 2007. Microsampling and isotopic analysis of igneous rocks: implications for the study of magmatic systems. Annu. Rev. Earth Planet. Sci. 35, 273–311. Deloule, E., Allegre, C.J., Doe, B.R., 1986. Lead and sulfur isotope microstratigraphy in galena crystals from Mississippi valley-type deposits. Econ. Geol. 81 (6), 1307–1321. DePaolo, D.J., 1981. Trace element and isotopic effects of combined wallrock assimilation and fractional crystallization. Earth Planet. Sci. Lett. 53, 189–202. DePaolo, D.J., 1985. Isotopic studies of processes in mafic magma chambers: I. The Kiglapait Intrusion, Labrador. J. Petrol. (4), 925–995. Dickin, A.P., Exley, R.A., 1981. Isotopic and geochemical evidence for magma mixing in the petrogenesis of the Coire Uaigneich Granophyre, Isle of Skye, N.W. Scotland. Contrib. Mineral. Petrol. 76, 98–108. Doe, B.R., 1970. B.R DoeLead Isotopes. Springer Verlag, Berlin, p. 137. Druitt, T.H., Costa, F., Deloule, E., Dungan, M., Scaillet, B., 2012. Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. Nature 482, 77–80. Duffield, W.A., Dalrymple, G.B., 1990. The Taylor Creek Rhyolite of New Mexico: a rapidly emplaced field of lava domes and flows. Bull. Volcanol. 52, 475–487. Duffield, W.A., du Bray, E.A., 1990. Temperature, size, and depth of the magma reservoir for the Taylor Creek Rhyolite, New Mexico. Am. Mineral. 75, 1059–1070. Duffield, W.A., Ruiz, J., 1992. Compositional gradients in large volume reservoirs of silicic magma as evidenced by ignimbrites versus Taylor Creek Rhyolite lava domes. Contrib. Mineral. Petrol. 110, 192–210. Duffield, W.A., Ruiz, J., 1994. A model to help explain Sr-isotope disequilibrium between feldspar phenocrysts and melt in large-volume silicic magma systems. US Geol. Surv. Open-File Rep. 94–423, 1–9. Espíndola, J.M., Macías, J.L., Tilling, R.I., Sheridan, M.F., 2000. Volcanic history of El Chichón Volcano (Chiapas, Mexico) during the Holocene, and its impact on human activity. Bull. Volcanol. 62, 90–104. Feldstein, S.N., Halliday, A.N., Davies, G.R., Hall, C.M., 1994. Isotope and chemical microsampling: constraints on the history of an S-type rhyolite, San Vincenzo, Tuscany, Italy. Geochim. Cosmochim. Acta 58, 943–958. Fujii, T., 1974. Crystal settling in a sill. Lithos 7, 133–137. Geist, D.J., Myers, J.D., Frost, C.D., 1988. Megacryst-bulk rock isotopic disequilibrium as an indicator of contamination process: the Edgcumbe Volcanic Field, SE Alaska. Contrib. Mineral. Petrol. 99, 105–112. Gray, N.H., Crain, I.K., 1969. Crystal settling in sills: a model for suspension settling. Can. J. Earth Sci. 6, 1211–1216. Hagen, H., Neumann, E.-R., 1990. Modeling of trace-element distribution in magma chambers using open-system models. Comput. Geosci. 16, 549–556. Halliday, A.N., Stephens, W.E., Hannon, R.S., 1980. Rb-Sr and O isotopic relationships in 3 zoned Caledonian granitic plutons, Southern Uplands, Scotland: evidence for varied sources and hybridisation of magmas. J. Geol. Soc. Lond. 137, 329–348. Hamada, M., Okayama, Y., Kaneko, T., Yasuda, A., Fujii, T., 2014. Polybaric crystallization differentiation of H2O-saturated island arc low-K tholeiite magmas: a case study of the Izu-Oshima volcano in the Izu arc. Earth Planets Space 66, 1–15. Hildreth, H., Moorbath, S., 1988. Crustal contributions to arc magmatism in the Andes of Central Chile. Contrib. Mineral. Petrol. 98, 455–489. Hobden, B.J., 1997. Modelling Magmatic Trends in Time and Space: Eruptive and Magmatic History of Tongariro Volcanic Complex, New Zealand. PhD thesis, University of Canterbury. Hobden, B.J., Houghton, B.F., Davidson, J.P., Weaver, S.D., 1999. Small and short-lived magma batches at composite volcanoes: time windows at Tongariro volcano, New Zealand. J. Geol. Soc. Lond. 156, 865–868. Hobden, B.J., Houghton, B.F., Nairn, I.A., 2002. Growth of a young, frequently active composite cone: Ngauruhoe volcano, New Zealand. Bull. Volcanol. 64, 392–409. Holden, P., Halliday, A.N., Stephens, W.E., 1987. Neodymium and strontium isotope content of microdiorite enclaves points to mantle input to granitoid production. Nature 330, 53–55. Huppert, H.E., Sparks, R.S.J., 1988. The generation of granitic magmas by intrusion of basalt into continental crust. J. Petrol. 29, 599–642. Kanno, T., Ishibashi, H., 2017. Pre-eruptive timescales estimated by diffusion chronometry of phenocryst minerals: a review. Geosci. Repts. Shizuoka Univ. 44, 47–64 (in Japanese with English abstract). Knesel, K.M., Davidson, J.P., Duffield, W.A., 1999. Evolution of silicic magma through assimilation and subsequent recharge: evidence from Sr isotopes in sanidine phenocrysts, Taylor Creek Rhyolite, NM. J. Petrol. 40, 773–786.

196

K. Nishimura / Journal of Volcanology and Geothermal Research 374 (2019) 181–196

Koyaguchi, T., Kaneko, K., 1999. A two-stage thermal evolution model of magmas in continental crust. J. Petrol. 40, 241–254. Koyaguchi, T., Kaneko, K., 2000. Thermal evolution of silicic magma chambers after basalt replenishments. Trans. R. Soc. Edinb. Earth Sci. 91, 47–60. Larsen, J.F., 2005. Experimental study of plagioclase rim growth around anorthite seed crystals in rhyodacitic melt. Am. Mineral. 90, 417–427. Mantovani, M.S.M., Hawkesworth, C.J., 1990. An inversion approach to assimilation and fractional crystallisation processes. Contrib. Mineral. Petrol. 105, 289–302. Marsh, B.D., 1981. On the crystallinity, probability of occurrences and rheology of lava and magma. Contrib. Mineral. Petrol. 78, 85–98. Marsh, B.D., 1988a. Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization: 1. Theory. Contrib. Mineral. Petrol. 99, 277–291. Marsh, B.D., 1988b. Crystal capture, sorting and retention in convecting magma. Geol. Soc. Am. Bull. 100, 1720–1737. Marsh, B.D., 1998. On the interpretation of crystal size distributions in magmatic systems. J. Petrol. 39, 553–599. Martin, D., Nokes, R., 1988. Crystal settling in a vigorously convecting magma chamber. Nature 332, 534–536. Mazzone, P., Grant, N.K., 1988. Mineralogical and isotopic evidence for phenocryst-matrix disequilibrium in the Garner Mountain andesite. Contrib. Mineral. Petrol. 99, 267–272. McGee, J.J., Tilling, R.I., 1983. 1982 and pre 1982 magmatic products of El Chichón volcano, Chiapas. Mexico. EOS Trans. Am. Geophys. Union 64, 893. McGee, J.J., Tilling, R.I., Duffield, W.A., 1987. Petrologic characteristics of the 1982 and pre 1982 eruptive products of El Chichón volcano, Chiapas. Mexico. Geofis. Int. 26-1, 85–108. Moore, A., Coogan, L.A., Costa, F., Perfit, M.R., 2014. Primitive melt replenishment and crystal-mush disaggregation in the weeks preceding the 2005–2006 eruption 9° 50′ N, EPR. Earth Planet. Sci. Lett. 403, 15–26. Nishimura, K., 2006. Numerical modeling of trace element behavior during crystal settling and reequilibration in high-silica magma bodies. J. Geophys. Res. 111, B08201. https://doi.org/10.1029/2005JB003844. Nishimura, K., 2009. A trace-element geochemical model for imperfect fractional crystallization associated with the development of crystal zoning. Geochim. Cosmochim. Acta 73, 2142–2149. Nishimura, K., 2012. A mathematical model of trace element and isotopic behavior during simultaneous assimilation and imperfect fractional crystallization. Contrib. Mineral. Petrol. 164, 427–440. Nishimura, K., 2013. AIFCCalc: an Excel spreadsheet for modeling simultaneous assimilation and imperfect fractional crystallization. Comput. Geosci. 51, 410–414. Nishimura, K., Kawamoto, T., Kobayashi, T., Sugimoto, T., Yamashita, S., 2005. Melt inclusion analysis of the Unzen 1991–1995 dacite: implications for crystallization processes of dacite magma. Bull. Volcanol. 67, 648–662. Powell, R., 1984. Inversion of the assimilation and fractional crystallization (AFC) equations; characterization of contaminants from isotope and trace element relationships in volcanic suites. J. Geol. Soc. Lond. 141, 447–452. Pupier, E., Duchene, S., Toplis, M.J., 2008. Experimental quantification of plagioclase crystal size distribution during cooling of a basaltic liquid. Contrib. Mineral. Petrol. 155, 555–570. Ramos, F.C., Wolff, J.A., Tollstrup, D.L., 2005. Sr isotope disequilibrium in Columbia River flood basalts: evidence for rapid shallow-level open-system processes. Geology 33, 457–460. Rayleigh, J.W.S., 1896. Theoretical considerations respecting the separation of gases by diffusion and similar processes. Philos. Mag. 42, 77–107. Richardson, J.F., Zaki, W.N., 1954. Sedimentation and fluidization. Part I. Trans. Inst. Chem. Eng. 32, 35–53. Roberts, M.P., Clemens, J.D., 1995. Feasibility of AFC models for the petrogenesis of calkalkaline magma series. Contrib. Mineral. Petrol. 121, 139–147. Ruprecht, P., Bergantz, G.W., Dufek, J., 2008. Modeling of gas-driven magmatic overturn: tracking of phenocryst dispersal and gathering during magma mixing. Geochem. Geophys. Geosyst. 9, Q07017. https://doi.org/10.1029/2008GC002022.

Saunders, K.E., Morgan, D.J., Baker, J.A., Wysoczanski, R.J., 2010. The magmatic evolution of the Whakamaru supereruption, New Zealand, constrained by a microanalytical study of plagioclase and quartz. J. Petrol. 51, 2465–2488. Shaw, H.R., 1980. The fracture mechanisms of magma transport from the mantle to the surface. In: Hargraves, R.B. (Ed.), Physics of Magmatic Processes. Princeton University Press, Princeton, NJ, pp. 201–264. Shaw, D.M., 2006. Trace Elements in Magmas: A Theoretical Treatment. Cambridge University Press, Cambridge. Simonetti, A., Bell, K., 1993. Isotopic disequilibrium in clinopyroxenes from nephelenitic lavas, Napak Volcano, eastern Uganda. Geology 21, 243–246. Singer, B.S., Costa, F., Herrin, J.S., Hildreth, W., Fierstein, J., 2016. The timing of compositionally-zoned magma reservoirs and mafic ‘priming’ weeks before the 1912 Novarupta-Katmai rhyolite eruption. Earth Planet. Sci. Lett. 451, 125–137. Snelling, A.A., 2003. The relevance of Rb-Sr, Sm-Nd and Pb-Pb isotope systematics to elucidation of the genesis and history of recent andesite flows at Mt Ngauruhoe, New Zealand, and the implications for radioisotopic dating. In: Ivey, R.L. Jr. (Eds.), Proceedings of the Fifth International Conference on Creationism, 285–303. Spera, F.J., Bohrson, W.A., 2001. Energy-constrained open-system magmatic processes I: general model and energy constrained assimilation and fractional crystallization (EC-AFC) formulation. J. Petrol. 42, 999–1018. Spera, F.J., Bohrson, W.A., 2002. Energy-constrained open-system magmatic processes 3. Energy-Constrained Recharge, Assimilation, and Fractional Crystallization (ECRAFC). Geochem. Geophys. Geosyst. 3 (12), 8001. https://doi.org/10.1029/ 2002GC000315. Spera, F.J., Bohrson, W.A., 2004. Open-system magma chamber evolution: an energyconstrained geochemical model incorporating the effects of concurrent eruption, recharge, variable assimilation and fractional crystallization (EC-E′RAχFC). J. Petrol. 45, 2459–2480. Steiner, A., 1958. Petrogenetic implications of the 1954 Ngauruhoe lava and its xenoliths. N. Z. J. Geol. Geophys. 1, 325–363. Taylor Jr., H.P., 1980. The effects of assimilation of country rocks by magmas on 18O/16O and 87Sr/86Sr systematics. Earth Planet. Sci. Lett. 47, 243–254. Tepley III, F.J., Davidson, J.P., Clynne, M.A., 1999. Magmatic interactions as recorded in plagioclase phenocrysts of Chaos Crags, Lassen Volcanic Center, California. J. Petrol. 40, 787–806. Tepley III, F.J., Davidson, J.P., Tilling, R.I., Arth, J.G., 2000. Magma mixing, recharge and eruption histories recorded in plagioclase phenocrysts from El Chichón Volcano, Mexico. J. Petrol. 41, 1397–1411. Till, C.B., Vazquez, J.A., Boyce, J.W., 2015. Months between rejuvenation and volcanic eruption at Yellowstone caldera, Wyoming. Geology 43, 695–698. Tilling, R.I., 2009. El Chichón's “surprise” eruption in 1982: lessons for reducing volcano risk. Geofis. Int. 48-1, 3–19. Tilling, R.I., Rubin, M., Sigurdsson, H., Carey, S., Duffield, W.A., Rose, W.I., 1984. Holocene eruptive activity of El Chichón Volcano, Chiapas, Mexico. Science 224, 747–749. Ubide, T., Kamber, B.S., 2018. Volcanic crystals as time capsules of eruption history. Nat. Commun. 9, 326. Verhoeven, J., Schmalzl, J., 2009. A numerical method for investigating crystal settling in convecting magma chambers. Geochem. Geophys. Geosyst. 10, Q12007. https://doi. org/10.1029/2009GC002509. Waight, T.E., Maas, R., Nicholls, I.A., 2000. Fingerprinting feldspar phenocrysts using crystal isotopic composition stratigraphy: implications for crystal transfer and magma mingling in S-type granites. Contrib. Mineral. Petrol. 139, 227–239. Wark, D.A., Hildreth, W., Spear, F.S., Cherniak, D.J., Watson, E.B., 2007. Pre-eruption recharge of the Bishop magma system. Geology 35, 235–238. Wittke, J.H., Duffield, W.A., Caron, J., 1996. Roof-rock contamination of Taylor Creek Rhyolite, New Mexico, as recorded in hornblende phenocrysts and biotite xenocrysts. Am. Mineral. 81, 135–140. Young, E.D., Wooden, J.L., Shieh, Y.-N., Farber, D., 1992. Geochemical evolution of Jurassic diorites from the Bristol Lake region, California, USA, and the role of assimilation. Contrib. Mineral. Petrol. 110, 68–86.