Chemical and isotopic effects in magmatic liquids evolving by AFC processes where the concentration in the assimilate changes by fractional melting (I). Analytical solutions

Abstract

Analytical solutions are proposed predicting the change in concentration and in isotopic composition of trace elements in magmatic liquids evolving by concurrent processes of fractional crystallization and wallrock assimilation (AFC), in the case the concentration of the trace element in the assimilate is assumed to change according to Shaw’s (1970) model of fractional melting of the wallrock. This particular AFC model, was first suggested and numerically approached by Spera and Bohrson (2001) in their energy-constrained AFC model (EC-AFC), without however giving analytical solutions.

The solutions here proposed are ‘natural extensions’ of the simpler solutions given by DePaolo (1981) and by Taylor and Sheppard (1986), turning to them if the distribution coefficient of the trace element between bulk residual wallrock and instantaneous melt is 1. However, as also shown by Bohrson and Spera’s calculation examples, these solutions generally predict concentration and isotopic evolution paths which are very different from those predicted by DePaolo’s and Taylor and Sheppard’s model. This is because, in addition to the trace element distribution coefficient between the residual solid and the instantaneous melt produced by melting of the wallrock, the analytical solutions contain a new parameter, which is the mass ratio between the initial magmatic liquid and the wallrock involved in the melting process.

Introduction

In the last five decades, petrological and geochemical studies of magmatic rocks have provided more and more evidences of chemical and isotopic changes in magmatic systems by melting and assimilation of rocks in magma chambers.

Changes in concentration and in isotopic ratios of trace elements in magmatic liquids evolving by concurrent processes of fractional crystallization and wallrock assimilation (AFC) have been quantitatively investigated by Minster and Allègre (1978), DePaolo (1981, 1985), Neumann et al. (1954) and Taylor and Sheppard (1986). Some AFC related aspects, i.e., relative assimilation, degree of contamination, thermal constraints and melting degree, were subsequently approached by Aitcheson and Forrest (1994), Bohrson and Spera (2001), Cavazzini (1996, 2010), Marsh (1988), Spera and Bohrson (2001) and Thompson et al. (2002).

Neumann et al. (1954) were the first authors to calculate the change in concentration of a trace element in a crystallizing magmatic liquid in the case the latent heat of crystallization causes re-adsorption of early-formed plagioclase crystals. In their mathematical treatment, these authors used for the first time a parameter P representing the instantaneous mass ratio between re-adsorption of the plagioclase phase and crystallization of the magmatic liquid.

As sometimes happens, this first quantitative approach seems to have been forgotten in the years 1955–1980, so that DePaolo (1981), and subsequently Taylor and Sheppard (1986), proposed analytical solutions for the evolution of the concentration of a trace element in the magmatic liquid which are exactly equivalent to those proposed by Neumann et al. (1954), because all these authors developed their mathematical models starting from identical conceptual premises (the ratio between the rate of assimilation and the rate of crystallization is called r instead of P, the distribution coefficient of the trace element is called D instead of K etc.; see eqs. (6)–(10) in Neumann et al., 1954). However, nor DePaolo (1981), nor Taylor and Sheppard (1986) provided an equation for the particular case where D = 1 – r, whereas Neumann et al. (1954) did (see eq. (10) in Neumann et al. (1954), and see also eqs. (22)–(24) in Cavazzini (1996), for the calculation of the degree of contamination in the particular case where D = 1 – r).

In all these first quantitative models, one important assumption is that the concentration of the trace element in the material which is assimilated by the magmatic liquid is constant as the melting of the wallrock proceeds. This assumption is very useful from the point of view of the calculation, because the starting differential equations are much easier to solve.

However, geochemical modelling of rock melting (Cavazzini, 2000; Hertogen and Gijbels, 1976; Shaw, 1970, 1978, 1979) indicates that the chemical composition of the very small mass of melt which is produced in a very small time-span (i.e., the ‘instantaneous’ melt) progressively changes as the wallrock melting degree increases. Assuming that mass dL of instantaneous melt is ‘assimilated’ (i.e., incorporated in the magmatic liquid, and thoroughly dispersed within the liquid to give a perfectly homogeneous liquid) as it forms, the model of Shaw (1970) for the instantaneous melt produced in a process of fractional melting may give a good description of the change in the chemical composition of the assimilated material, because this model of melting assumes that the amount of melt produced in a very small time-span is removed from the source rock as it forms.

If Shaw’s fractional melting model is considered, the concentration of a trace element in the instantaneous melt not only may be very different from the concentration in the initial wallrock, but also it may change significantly as the melting proceeds. For this reason, in general, chemical and isotopic evolution paths of a magmatic liquid which evolves by concurrent processes of assimilation and fractional crystallization are expected to be very different from those predicted by the quantitative models of DePaolo (1981) and of Taylor and Sheppard (1986). This aspect was approached by Spera and Bohrson (2001) and Bohrson and Spera (2001) in their energy-constrained model of magmatic evolution by combined wallrock assimilation and fractional crystallization (EC-AFC). These authors calculated the chemical and the isotopic changes in magmatic liquids assuming the concentration of the trace element in the instantaneous assimilate to change according to Shaw’s (1970) model of fractional melting. However, their description is computational (numerical), and they did not provide analytical solutions as Neumann et al. (1954), DePaolo (1981) and Taylor and Sheppard (1986) did.

Analytical solutions, however, would be particularly desirable, not only because they would be easy-to-use (as the solutions proposed by DePaolo (1981) and by Taylor and Sheppard (1986)), but also because they could be readily used in further quantitative modelling, and the aim of this study is proposing analytical solutions to calculate the evolution of the concentration and of the isotopic ratios of trace elements in magmatic liquids evolving by AFC processes in the case the assimilated material is the instantaneous melt produced by fractional melting of the wallrock according to Shaw’s (1970) model. Therefore, in what follows, calculations are illustrated that result from a work of ‘interfacing’ the equations of Shaw (1970) for the instantaneous melt in a fractional partial melting of a rock with DePaolo’s (1981) (or Taylor and Sheppard’s (1986), they are equivalent) AFC differential equations for the concentration and the isotopic ratio of a trace element in the magmatic liquid. The reader will find that the analytical solutions we propose are ‘natural extensions’ of those proposed by DePaolo (1981) and by Taylor and Sheppard (1986).

The solutions we propose in this study correspond, with minor improvements, to the solutions used by Cavazzini (2010) to predict the concentration vs. isotopic ratio patterns for Sr, Nd and the wallrock melting degrees calculated by Bohrson and Spera’s (2001) EC-AFC model. In this study, we present the differential equations and their analytical solutions, whereas in a next paper we will show i) how these solutions work in predicting concentration vs. isotopic ratio patterns of trace elements in magmatic liquids evolving by AFC processes where the concentration of the trace element in the instantaneous assimilate is assumed to change according to Shaw’s model of fractional melting (e.g., Sr and Nd evolution patterns in the ‘Upper crustal’ and in the ‘Lower crustal’ cases illustrated in Figs. 3a, 4a and 5 in Bohrson and Spera (2001); see also Cavazzini (2010)); ii) how they work in predicting the wallrock melting degree in the two cases cited above (Bohrson and Spera, 2001; Cavazzini, 2010); and iii) how they can be used to significantly constrain interpretive scenarios in AFC evolution contexts (Cavazzini, 2010).

Previous Work

Shaw’s model of fractional melting

Shaw’s (1970) equation for the concentration of a trace element in the melt produced at an instant during a fractional melting process is:

  

$$C_{ist}^l=\frac{C_{s,0}}{D_0}(1-f{)}^{\frac{1}{D_0}-1}$$(1)

in the case of modal melting, and

  

$$C_{ist}^l=\frac{C_{s,0}}{D_0}{(1-\frac{P}{D_0}f)}^{\frac{1}{P}-1}$$(2)

in the case of non-modal melting, where:

$C_{ist}^l$ is the concentration of the element in the mass of melt produced in a very small time-span dt;

$f$ is the melting degree of the wallrock at that instant, i.e. the ratio between the total mass L of melt produced from the beginning of the melting process and the initial mass $W_0$ of rock which is involved in the melting process;

$C_{s,0}$ is the concentration of the element in the wall rock before the melting process begins;

$D_0$ is the coefficient which, at any instant, describes the distribution of the trace element between bulk residual solid and melt produced at that instant. (This coefficient has a constant value in the case of a modal melting process, whereas it represents a value at initial conditions in the case of a non-modal melting process).

$P$ is a parameter used by Shaw (1970) in the mathematical description of a non-modal process. It can be briefly described as the coefficient that at any instant describes the distribution of the element between the bulk solid and the melt produced at that instant, but referred to the fictitious solid which melts modally ‘within’ the bulk solid which actually melts non-modally.

DePaolo’s model of magmatic assimilation-fractional crystallization

DePaolo’s (1981) (or Taylor and Sheppard’s (1986), they are equivalent) AFC differential equations are reported below. Two cases can be distinguished depending on if ratio r (the ratio between the rate of assimilated mass and the rate of crystallized mass) is different from 1 or not. r ≠ 1: In the most general case, the change in the concentration of a trace element is given by the solution of the following differential equation:

  

$$\frac{dC}{d\ln F}+zC=\left(\frac{r}{r-1}\right)C_a.$$(3)

The change in an isotopic ratio of the element (e.g., ⁸⁷Sr/⁸⁶Sr, ¹⁴³Nd/¹⁴⁴Nd, …) is given by the solution of the following differential equation:

  

$$\frac{d\epsilon }{d\ln F}=\left(\frac{r}{r-1}\right)\frac{C_a}{C}({\epsilon }_a-\epsilon ),$$(4)

where:

$C$ is the concentration of the element in the magmatic liquid at an instant (i.e., the solution of eq. (3));

$C_a$ is the concentration of the element in the assimilated material at that instant;

$F$ is the mass fraction of magmatic liquid (i.e., the ratio between the mass of liquid at that instant and the mass of liquid at the beginning of the crystallization-assimilation process);

$r$ is the mass ratio between assimilation and crystallization at that instant;

$\epsilon$ is the value of an isotopic ratio of the trace element in the liquid at that instant;

${\epsilon }_a$ is the value of the same isotopic ratio in the assimilated material at that instant;

$z=(r+D-1)/(r-1)$, is a parameter where

$D$ is the instantaneous value of the distribution coefficient of the element between the bulk crystallizing assemblage and the magmatic liquid.

r = 1: In this case, the mass fraction of liquid cannot be used as the independent variable (DePaolo (1981)). The differential equation which gives the change in concentration is:

  

$$\frac{dC}{dt}=\frac{r_a}{M}(C_a-DC),$$(5)

and the change in the value of an isotopic ratio is given by the differential equation:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\frac{r_a}{M}\frac{C_a}{C}dt,$$(6)

where (DePaolo, 1981):

$r_a$ is the instantaneous rate of assimilation; and

$M$ is the mass of the magmatic liquid (a constant over the process).

Relationship Between the Mass Fraction of Magmatic Liquid and the Wallrock Melting Degree

By observing the equations above, it is clear that direct substitutions of eqs. (1) and (2) into eqs. (3)–(6) preclude immediate solutions, because Shaw’s equations give the concentration of the element in the instantaneous melt as a function of the melting degree of the wallrock involved in the melting process, whereas in the AFC differential equations of DePaolo (1981) the changes in the value of the concentration and of the isotopic ratio of the element in the magmatic liquid are given in terms of the change in the liquid mass fraction, F, (case of r ≠ 1), or in terms of the change in the total assimilated mass a (case of r = 1).

Thus, “interfacing” Shaw’s (1970) equations for fractional melting and DePaolo’s (1981) AFC differential equations is necessary, and requires providing a relationship between the melting degree of the wallrock and the mass fraction of liquid (r ≠ 1), or the total assimilated mass (r = 1).

From a physical point of view, a relationship between the melting degree of the wallrock and the fraction of residual magmatic liquid is suggested by heat-balance considerations, because melting is driven by the latent heat of crystallization of the magmatic liquid (Bohrson and Spera, 2001; Bowen, 1928; Spera and Bohrson, 2001; Taylor, 1980).

In calculating such mathematical relationship, it cannot be forgotten that using Shaw’s (1970) equations implies that the assumptions involved in Shaw’s melting model will be maintained. If it does not, the treatment would lose internal consistency.

The mass fraction of magmatic liquid at an instant is:

  

$$F=\frac{M}{M_0}$$(7)

where M is the mass of liquid at that instant, and M₀ is the mass of liquid at the beginning of the process of concurrent crystallization and assimilation.

Differentiating (7) we obtain:

  

$$dM=M_{0}dF.$$(8)

Due to concurrency of processes of assimilation and crystallization, we can also write the change of the mass of magmatic liquid in a very small time-span as:

  

$$dM=da-dc$$(9)

where da, dc are the mass of liquid added by assimilation and the mass of liquid removed by crystallization in that instant of time, respectively. Thus, using ratio r between the rate of assimilation and the rate of crystallization in that time instant, eq. (9) can be written as:

  

$$dM=\left(\frac{r-1}{r}\right)da.$$(10)

Now (Shaw, 1970), the wallrock melting degree is defined as:

  

$$f=\frac{L}{W_0}$$(11)

where L is, at an instant, the total (accumulate) mass of melt produced from the beginning of the melting process, and W₀ is the mass of wallrock involved in the melting process.

Differentiating (11), and assuming that mass dL of melt produced in a very small time-span is totally and instantaneously assimilated (i.e., da = dL), eqs. (8) and (10) can be used to write:

  

$$M_{0}dF=\left(\frac{r-1}{r}\right)W_{0}df.$$(12)

Integrating eq. (12) on a finite time-span of melting of the wallrock and of crystallization of the magmatic liquid, with r as a constant (as done by Neumann et al. (1954), DePaolo (1981) and Taylor and Sheppard (1986)) and W₀ as a constant (this condition is required for internal consistency with the calculations in Shaw’s (1970) melting model) we obtain:

  

$$f=\frac{M_0}{W_0}\left(\frac{r}{r-1}\right)(F-1).$$(13)

Note that eq. (13) cannot be used in the particular case r = 1, because it involves an indeterminate. As shown by DePaolo (1981), and illustrated above, in the case of r = 1 the mass fraction of magmatic liquid cannot be the independent variable, because the mass of magmatic liquid does not change during the process.

Case r ≠ 1: Concentration

Modal melting

In the case of modal melting, substituting eqs. (1) and (13) into eq. (3) we obtain the following differential equation:

  

$$\frac{dC}{d\ln F}+zC=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}{[1-\mu \left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{D_0}-1}.$$(14)

in which μ = M₀/W₀. The solution of this differential equation is

  

$$C=[C_0-\frac{a}{z}{\left(\frac{1}{b+1}\right)}^{-d}•{}_2{F}_1(-d,z;z+1;\frac{b}{b+1})]F^{-z}+\frac{a}{z}(1+b-bF{)}^d{(1+\frac{bF}{b+1})}^{-d}•{}_2{F}_1(-d,z;z+1;\frac{bF}{b+1})$$(15)

where

  

$$a=\frac{C_a}{D_0}\left(\frac{r}{r-1}\right)b=\mu \left(\frac{r}{r-1}\right)d=\frac{1}{D_0}-1.$$(16)

Solution (15) involves two ₂F₁ hypergeometric functions. A hypergeometric function is a power series

  

$${\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+{\beta }_3{x}^3+\dots =\sum _{n=0}^{+\infty }{\beta }_n{x}^n,$$

in which β₀ = 1, and the ratio between the coefficients of two subsequent powers xⁿ⁺¹ and xⁿ is a rational function of n:

  

$$\frac{{\beta }_{n+1}}{{\beta }_n}=\frac{A(n)}{B(n)}.$$

For example, in the case of the second hypergeometric function, we have:

  

$${}_2{F}_1(-d,z;z+1;\frac{bF}{b+1})=1+\left(\frac{-dz}{z+1}\right)\frac{bF}{b+1}+\frac{[-d(1-d)][z(1+z)]}{(z+1)(1+z+1)}{\left(\frac{bF}{b+1}\right)}^2+\frac{[-d(1-d)(2-d)][z(1+z)(2+z)]}{(z+1)(1+z+1)(2+z+1)}{\left(\frac{bF}{b+1}\right)}^3+\dots$$

In synthesis:

  

$${}_2{F}_1(-d,z;z+1;\frac{bF}{b+1})=\sum _{n=0}^{+\infty }\frac{(-d{)}_n(z{)}_n}{(z+1{)}_n}{\left(\frac{bF}{b+1}\right)}^n,$$

where (–d)_n, (z)_n e (z+1)_n are called ‘Pochhammer’s symbols’.

Solution (15) can be used to exactly calculate the evolution of the concentration of the trace element. However, an hypergeometric function is a series, and we do not have an analytical solution in the mathematical sense. For this reason, we worked to approximate parts of eq. (14), and two possible ways can be found, which are illustrated in the following as Mod. 1 and Mod. 2. The first approximation gives an analytical solution in the common mathematical sense. The second approximation gives a series. In this case, from the point of view of the accuracy, the calculation is ‘perfectible’, and it can be considered substantially accurate if the series is truncated after the term of the fifth order. The accuracy of the values of concentration and isotopic ratio calculated by these two solutions substantially depend on the behavior of the trace element in the process of magmatic crystallization on one side and in the process of melting of the wallrock on the other. In the case of Mod. 1, the differential equations for the concentration and for the isotopic ratios are both quite easily solved.

Mod. 1 considers the binomial expansion:

  

$${[1+(F-1)]}^{-\mu \left(\frac{r}{r-1}\right)}=1+[-\mu \left(\frac{r}{r-1}\right)](F-1)+\frac{1}{2!}\{-\mu \left(\frac{r}{r-1}\right)[-\mu \left(\frac{r}{r-1}\right)-1]\}(F-1{)}^2+\dots .$$(17)

If the conditions are such that good accuracy is maintained by truncating this series after the term of the first-order, eq. (14) is simplified to:

  

$$\frac{dC}{d\ln F}+zC=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}{\left(e^{\ln F}\right)}^{-\mu \left(\frac{r}{r-1}\right)(\frac{1}{D_0}-1)}.$$(18)

This differential equation is solved by the integrating factor $e^{{\int }_{}^{}zd\ln F}=e^{z\ln F}$:

  

$$C•e^{z\ln F}=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}{\int }_{}^{}{e}^{(q+z)\ln F}d\ln F+Cost$$(19)

where Cost is the essential arbitrary constant, and

  

$$q=-\mu \left(\frac{r}{r-1}\right)(\frac{1}{D_0}-1).$$(20)

Solving the integral, we calculate the value of the arbitrary constant by the initial conditions (lnF = 0):

  

$$Cost=C_0-\frac{C_{s,0}}{D_0}\left(\frac{r}{r-1}\right)\frac{1}{q+z}$$(21)

where C₀ is the concentration of the element in the magmatic liquid at the beginning of the process. Finally, substituting (21) in (19), we re-arrange to:

  

$$C=C_0{F}^{-z}+\frac{r{C}_{s,0}}{\mu (D_0-1)r+D_0(r+D-1)}(F^q-F^{-z}).$$(22)

Note that eq. (22) does not involve an indeterminate for the particular case z = 0 (D = 1 – r). An indeterminate occurs for z = 0 in the corresponding equation proposed by Neumann et al. (1954), DePaolo (1981) and Taylor and Sheppard (1986), and for this particular case Neumann et al. (1954) gave a solution. Instead, in eq. (22) it is sufficient to substitute z = 0 and D = 1 – r, obtaining the simplified formula:

  

$$C=C_0+\frac{C_{s,0}}{\mu (D_0-1)}(F^q-1).$$(23)

Mod. 2 considers that since the wallrock melting degree is in the range 0–1, we are sure that the power at the second member in eq. (14) can be expressed as a variant of Newton’s binomial formula:

  

$${[1-\mu \left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{D_0}-1}=1+[-\mu \left(\frac{r}{r-1}\right)(F-1)](\frac{1}{D_0}-1)+\frac{1}{2!}{[-\mu \left(\frac{r}{r-1}\right)(F-1)]}^2(\frac{1}{D_0}-1)(\frac{1}{D_0}-2)+\dots .$$(24)

We can approximate series (24) by truncation, for example after the fifth-order term. Then, after solution of the resulting differential equation, we obtain:

  

$$C=C_0{F}^{-z}+\frac{G}{5+z}(F^5-F^{-z})+\frac{H}{4+z}(F^4-F^{-z})+\frac{I}{3+z}(F^3-F^{-z})+\frac{J}{2+z}(F^2-F^{-z})+\frac{K}{1+z}(F^{}-F^{-z})+\frac{L}{z}(1-F^{-z})$$(25)

where:

  

$$G=Q•E$$

  

$$H=Q(D-5E)$$

  

$$I=Q(C-4D+10E)$$

  

$$J=Q(B-3C+6D-10E)$$

  

$$K=Q(A-2B+3C-4D+5E)$$

  

$$L=Q(1-A+B-C+D-E)$$

and

  

$$Q=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}$$

  

$$A=(\frac{1}{D_0}-1)[-\mu \left(\frac{r}{r-1}\right)]$$

  

$$B=\frac{1}{2!}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2){[-\mu \left(\frac{r}{r-1}\right)]}^2$$

  

$$C=\frac{1}{3!}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2)(\frac{1}{D_0}-3){[-\mu \left(\frac{r}{r-1}\right)]}^3$$

  

$$D=\frac{1}{4!}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2)(\frac{1}{D_0}-3)(\frac{1}{D_0}-4){[-\mu \left(\frac{r}{r-1}\right)]}^4$$

  

$$E=\frac{1}{5!}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2)(\frac{1}{D_0}-3)(\frac{1}{D_0}-4)(\frac{1}{D_0}-5){[-\mu \left(\frac{r}{r-1}\right)]}^5$$

We emphasize that Mod. 1 (eq. (22)) and Mod. 2 (eq. (25)) are not equivalent, because they are the results of different premises. This means that in some cases Mod. 1 may describe the evolutive paths in the magmatic liquid more accurately than Mod. 2, and vice-versa, depending on the respective deviations of the integrations using eq. (17) or eq. (24) from the integration of differential equation (14). Truncation of series (24) after the term of the fifth order generally provides accurate calculations of the concentration of the trace element in the magmatic liquid over a large extent of the evolution process.

For example, assuming μ = 5, r = 0.1, and D = 2, in the case the trace element behaves as an incompatible element in partial melting, Mod. 1 is much more accurate than Mod. 2 in describing the evolution of the concentration in the magmatic liquid. This is because in Mod. 2 the approximation of the integrand function assumes low values of the exponent 1/D₀ – 1. However, in the case D₀ is low, for example in the range 0.05–0.3, the exponent is too much high for the approximation to be accurate.

Non-modal melting

In non-modal melting, substituting eqs. (2) and (13) into eq. (3), we have

  

$$\frac{dC}{d\ln F}+zC=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}{[1-\frac{P}{D_0}\mu \left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{P}-1}.$$(26)

As we have shown in the case of modal melting, to solve this differential equation we can follow two possible ways. We first consider the binomial expansion:

  

$${[1+(F-1)]}^{-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)}=1+[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)](F-1)+\frac{1}{2!}\{-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)-1]\}(F-1{)}^2+.\dots$$(27)

If the conditions are such that sufficient accuracy is maintained by truncating this series after the term of the first-order, eq. (26) can be approximated to:

  

$$\frac{dC}{d\ln F}+zC=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}{\left(e^{\ln F}\right)}^{-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)(\frac{1}{P}-1)}.$$(28)

As shown in the case of modal melting, eq. (28) is solved using the integrating factor $e^{{\int }_{}^{}zd\ln F}=e^{z\ln F}$. After calculation of the essential arbitrary constant, calling

  

$$q=-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)(\frac{1}{P}-1)$$(29)

we re-arrange to:

  

$$C=C_0{F}^{-z}+\frac{r{C}_{s,0}}{\mu (P-1)r+D_0(r+D-1)}(F^q-F^{-z})$$(30)

which is structurally identical to eq. (22).

The second possible way, is analogous to that we have proposed above in the case of a process of modal melting. Since the wallrock melting degree is in the range 0–1, if ratio P/D₀ ≤ 1 or if it is >1 but not much different from 1, the power at the second member in eq. (26) can be expressed in terms of a Taylor’s series:

  

$${[1-\frac{P}{D_0}\mu \left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{P}-1}=1+(\frac{1}{P}-1)[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)(F-1)]+\frac{1}{2!}{(\frac{1}{P}-1)(\frac{1}{P}-2)[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)(F-1)]}^2+\dots .$$(31)

Series (31) can be approximated by truncating it after the fifth-order term, for good accuracy. Therefore, the integration of the differential equation gives:

  

$$C=C_0{F}^{-z}+\frac{G}{5+z}(F^5-F^{-z})+\frac{H}{4+z}(F^4-F^{-z})+\frac{I}{3+z}(F^3-F^{-z})+\frac{J}{2+z}(F^2-F^{-z})+\frac{K}{1+z}(F^{}-F^{-z})+\frac{L}{z}(1-F^{-z})$$(32)

where:

  

$$G=Q•E$$

  

$$H=Q(D-5E)$$

  

$$I=Q(C-4D+10E)$$

  

$$J=Q(B-3C+6D-10E)$$

  

$$K=Q(A-2B+3C-4D+5E)$$

  

$$L=Q(1-A+B-C+D-E)$$

and

  

$$Q=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}$$

  

$$A=(\frac{1}{P}-1)[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)]$$

  

$$B=\frac{1}{2!}(\frac{1}{P}-1)(\frac{1}{P}-2){[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)]}^2$$

  

$$C=\frac{1}{3!}(\frac{1}{P}-1)(\frac{1}{P}-2)(\frac{1}{P}-3){[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)]}^3$$

  

$$D=\frac{1}{4!}(\frac{1}{P}-1)(\frac{1}{P}-2)(\frac{1}{P}-3)(\frac{1}{P}-4){[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)]}^4$$

  

$$E=\frac{1}{5!}(\frac{1}{P}-1)(\frac{1}{P}-2)(\frac{1}{P}-3)(\frac{1}{P}-4)(\frac{1}{P}-5){[-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)]}^5$$

Eq. (32) is structurally identical to eq. (25). Difference is due to different meaning of D₀ in the modal melting equation with respect to the non-modal melting one. In Shaw’s (1970) model, solutions for non-modal melting processes are fundamentally given in terms of modal melting of fictitious solids identified ‘within’ the solids which are subjected to non-modal melting (i.e., modal melting can be considered as a particular case of non-modal melting where P = D = D₀).

Case r ≠ 1: Isotopic Ratios

Now equations are available which give the change in the concentration of the trace element in the magmatic liquid during the process. Thus, we can use these equations inside differential equation (4) to calculate the evolution of the value of an isotopic ratio of the element.

In the following calculations, it is assumed that the value of the isotopic ratio in the material which is assimilated does not change as the melting degree increases. Cavazzini (2000) has shown that this occurs if (a) the proportions of the mineral phases which melt do not change during the process; (b) the chemical fractionation coefficients of the trace element in the mineral phases which melt do not change during the process; and, if the isotopic composition of the trace element changes over time due to processes of radioactive decay of other nuclides, (c) if the duration of the process is such that no significant change in the isotopic composition of the trace element occurs in the mineral phases which participate to melting (depending on respective radioactive parent/radiogenic daughter ratios).

Modal melting

In the case of a process of modal melting, substituting eqs. (1), (13) and (22) into eq. (4), we obtain the following differential equation:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\frac{{[1-\mu \left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{D_0}-1}}{C_0{F}^{-z}+a(F^q-F^{-z})}d\ln F$$(33)

where

  

$$a=\frac{r{C}_{s,0}}{\mu (D_0-1)r+D_0(r+D-1)}$$(34)

and q is as in eq. (20).

We can solve eq. (33) if the term within the square brackets on the right side of the equation is a very good approximation of the function

  

$${[1+(F-1)]}^{-\mu \left(\frac{r}{r-1}\right)}.$$

In this case, the equation simplifies to:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\frac{F^q}{C_0{F}^{-z}+a(F^q-F^{-z})}d\ln F$$(35)

and can be re-arranged to:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\frac{1}{a+(C_0-a)e^{-(q+z)\ln F}}d\ln F.$$(36)

This equation can be integrated to:

  

$$-\ln \left({\epsilon }_a-\epsilon \right)+\ln \left({\epsilon }_a-{\epsilon }_0\right)=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\left\{\begin{array}{c}\frac{1}{a(q+z)}\ln \left[a{e}^{(q+z)\ln F}+(C_0-a)\right]+\\ -\frac{1}{a(q+z)}\ln \left[a{e}^{(q+z)\ln 1}+(C_0-a)\right]\end{array}\right\}$$(37)

and opportunely re-arranged to:

  

$$\frac{{\epsilon }_a-\epsilon }{{\epsilon }_a-{\epsilon }_0}=\frac{C_0{e}^{-(q+z)\ln F}}{a+(C_0-a)e^{-(q+z)\ln F}},$$(38)

and, finally, to

  

$$\epsilon ={\epsilon }_a-({\epsilon }_a-{\epsilon }_0)\frac{C_0}{C_0+a(F^{q+z}-1)}.$$(39)

Eq. (39) corresponds to a degree of contamination of the trace element (Cavazzini, 1996) given by the following formula:

  

$${\delta }_c=\frac{\epsilon -{\epsilon }_0}{{\epsilon }_a-{\epsilon }_0}=1-\left(\frac{C_0}{C_0{F}^{-z}+a(F^q-F^{-z})}\right)F^{-z}.$$(40)

As in the case of eq. (22), also eq. (39) does not involve an indeterminate for z = 0. It simplifies to:

  

$$\epsilon ={\epsilon }_a-({\epsilon }_a-{\epsilon }_0)\frac{C_0}{C_0+\frac{C_{s,0}}{\mu (D_0-1)}(F^q-1)}$$(41)

Non-modal melting

In the case of a process of non-modal melting, substituting eqs. (2), (13) and (30) into eq. (4), we obtain the following differential equation:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\frac{{[1-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)(F-1)]}^{\frac{1}{P}-1}}{C_0{F}^{-z}+a(F^q-F^{-z})}d\ln F$$(42)

where a is different than in (34):

  

$$a=\frac{r{C}_{s,0}}{\mu (P-1)r+D_0(r+D-1)},$$(43)

and q is as in eq. (29).

We can solve eq. (42) if the term within the square brackets on the right side of the equation is a very good approximation of the function

  

$${[1+(F-1)]}^{-\mu \frac{P}{D_0}\left(\frac{r}{r-1}\right)}.$$

In this case, the equation simplifies to:

  

$$\frac{d\epsilon }{{\epsilon }_a-\epsilon }=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\frac{F^q}{C_0{F}^{-z}+a(F^q-F^{-z})}d\ln F$$(44)

which is structurally identical to (35) and can be integrated to

  

$$-\ln \left({\epsilon }_a-\epsilon \right)+\ln \left({\epsilon }_a-{\epsilon }_0\right)=\left(\frac{r}{r-1}\right)\frac{C_{s,0}}{D_0}\left\{\begin{array}{c}\frac{1}{a(q+z)}\ln \left[a{e}^{(q+z)\ln F}+(C_0-a)\right]+\\ -\frac{1}{a(q+z)}\ln \left[a{e}^{(q+z)\ln 1}+(C_0-a)\right]\end{array}\right\},$$(45)

finally giving a solution which is structurally identical to the solution obtained for the case of modal melting:

  

$$\epsilon ={\epsilon }_a-({\epsilon }_a-{\epsilon }_0)\frac{C_0}{C_0+a(F^{q+z}-1)}.$$(46)

Case r = 1: Concentration

In this particular case, we can provide useful analytical solutions for the evolution of the concentration of a trace element, but no satisfactory solution to predict the change in an isotopic ratio of a trace element.

If r = 1, eq. (13) cannot be used. We need an equation where the wallrock melting degree is expressed as a function of a different parameter. As the independent variable, we can choose the total amount of assimilated mass or, that is the same, the ratio between the total assimilated mass and the mass of the magmatic liquid, because the latter is a constant. This was the approach of DePaolo (1981), and since it is useful to provide equations easily comparable with the equations of DePaolo, in the following the same parameter is adopted.

If the mass of melt produced in a short time-span is completely assimilated, then da = dL and integrating we have:

  

$$f=\frac{a}{W_0}.$$(47)

Modal melting

Since $r_a=\frac{da}{dt}$, differential eq. (5) can be written as:

  

$$\frac{dC}{da}+\frac{D}{M}C=\frac{C_a}{M}$$(48)

and substituting (1):

  

$$\frac{dC}{da}+\frac{D}{M}C=\frac{1}{M}\frac{C_{s,0}}{D_0}{(1-\frac{a}{W_0})}^{\frac{1}{D_0}-1}.$$(49)

Eq. (49) can be solved by mean of the integrating factor $e^{{\int }_{}^{}\frac{D}{M}da}=e^{\frac{D}{M}a}$:

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}{\int }_{}^{}{(1-\frac{a}{W_0})}^{\frac{1}{D_0}-1}•e^{\frac{D}{M}a}da+Cost$$(50)

Since f is in the range 0–1, we can use Taylor’s binomial expansion for the quantity

  

$${(1-\frac{a}{W_0})}^{\frac{1}{D_0}-1},$$

and if the conditions exist such that this series can be usefully approximated by truncating it to the term of the second order, we have:

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}•{\int }_{}^{}\{[1+(\frac{1}{D_0}-1)(-\frac{a}{W_0})+\frac{1}{2!}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2){(-\frac{a}{W_0})}^2]•e^{\frac{D}{M}a}\}da+Cost$$(51)

This equation can be changed to:

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}\left\{\begin{array}{c}{\displaystyle \int e^{\frac{D}{M}a}da+\left(\frac{1}{D_0}-1\right)\left(-\frac{1}{W_0}\right)}{\displaystyle \int \left(a•e^{\frac{D}{M}a}\right)da}\\ +\frac{1}{2}\left(\frac{1}{D_0}-1\right)\left(\frac{1}{D_0}-2\right){\left(-\frac{1}{W_0}\right)}^2{\displaystyle \int \left(a^2•e^{\frac{D}{M}a}\right)da}\end{array}\right\}+Cost$$(52)

Solving the integrals, we finally have:

  

$$C=C_0{e}^{-\frac{D}{M}a}+\frac{C_{s,0}}{D_{0}D}[1-e^{-\frac{D}{M}a}-(\frac{1}{D_0}-1)\left(\frac{a}{M}\right)\mu +\frac{1}{2}(\frac{1}{D_0}-1)(\frac{1}{D_0}-2){\left(\frac{a}{M}\right)}^2{\mu }^2]+\frac{C_{s,0}}{D_0{D}^2}\mu (\frac{1}{D_0}-1)[1-e^{-\frac{D}{M}a}-(\frac{1}{D_0}-2)\left(\frac{a}{M}\right)\mu ]+\frac{C_{s,0}}{D_0{D}^3}{\mu }^2(\frac{1}{D_0}-1)(\frac{1}{D_0}-2)(1-e^{-\frac{D}{M}a}),$$(53)

which corresponds to the instantaneous degree of contamination:

  

$${\delta }_c=1-\frac{C_0}{C}{e}^{\frac{-Da}{M}}.$$(54)

For more accurate solutions, the binomial series can be progressively truncated after the term of the third, or of the fourth, or of the fifth order, etc.

Non-modal melting

Substituting equation (2) into eq. (48):

  

$$\frac{dC}{da}+\frac{D}{M}C=\frac{1}{M}\frac{C_{s,0}}{D_0}{(1-\frac{P}{D_0}\frac{a}{W_0})}^{\frac{1}{P}-1}.$$(55)

As in the case of eq. (49), eq. (55) can be solved by mean of the integrating factor $e^{{\int }_{}^{}\frac{D}{M}da}=e^{\frac{D}{M}a}:$

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}{\int }_{}^{}{(1-\frac{P}{D_0}\frac{a}{W_0})}^{\frac{1}{P}-1}•e^{\frac{D}{M}a}da+Cost,$$(56)

and if the conditions are such that $\left(\frac{P}{D_0}\right)\left(\frac{a}{W_0}\right)$ is between 0 and 1 and the Taylor’s series is well approximated although truncated to the term of the second order, we have:

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}•{\int }_{}^{}[1+(\frac{1}{P}-1)(-\frac{P}{D_0}\frac{a}{W_0})+\frac{1}{2!}(\frac{1}{P}-1)(\frac{1}{P}-2){(-\frac{P}{D_0}\frac{a}{W_0})}^2]•e^{\frac{D}{M}a}da+Cost.$$(57)

This equation can be changed to:

  

$$C•e^{\frac{D}{M}a}=\frac{1}{M}\frac{C_{s,0}}{D_0}\left\{\begin{array}{c}{\displaystyle \int e^{\frac{D}{M}a}da+\left(\frac{1}{P}-1\right)\left(-\frac{P}{D_0{W}_0}\right)}{\displaystyle \int \left(a•e^{\frac{D}{M}a}\right)da}\\ +\frac{1}{2}\left(\frac{1}{P}-1\right)\left(\frac{1}{P}-2\right){\left(-\frac{P}{D_0{W}_0}\right)}^2{\displaystyle \int \left(a^2•e^{\frac{D}{M}a}\right)da}\end{array}\right\}+Cost$$(58)

and solved to give:

  

$$C=C_0{e}^{-\frac{D}{M}a}+\frac{C_{s,0}}{D_{0}D}[1-e^{-\frac{D}{M}a}-(\frac{1}{P}-1)\left(\frac{P}{D_0}\right)\left(\frac{a}{M}\right)\mu +\frac{1}{2}(\frac{1}{P}-1)(\frac{1}{P}-2){\left(\frac{P}{D_0}\right)}^2{\left(\frac{a}{M}\right)}^2{\mu }^2]+\frac{C_{s,0}}{D_0{D}^2}\mu (\frac{1}{P}-1)\left(\frac{P}{D_0}\right)[1-e^{-\frac{D}{M}a}-(\frac{1}{P}-2)\left(\frac{P}{D_0}\right)\left(\frac{a}{M}\right)\mu ]+\frac{C_{s,0}}{D_0{D}^3}{\mu }^2(\frac{1}{P}-1)(\frac{1}{P}-2){\left(\frac{P}{D_0}\right)}^2(1-e^{-\frac{D}{M}a}).$$(59)

As in the case of (53), for more accurate solutions, the Taylor’s series can be truncated after the term of the third, fourth, fifth order, etc.

Discussion and Calculation Examples

As in the case of DePaolo’s (1981) and of Taylor and Sheppard’s (1986) AFC model, the mathematical treatment presented in this study describes an ideal process of contamination of the liquid phase in a magma by pure assimilation, i.e., assuming that the melt generated by melting of wallrock is perfectly dispersed thoroughly within the magmatic liquid to give a homogeneous liquid. It is reasonably expected that this can occur only for low melting degrees of the wallrock, because for higher melting degrees contamination in a magma does not just occur by assimilation of rock material in the liquid phase. Rock fragments would be incorporated in the magma to give complex chemical and isotopic patterns.

These solutions are natural extensions of the solutions calculated by DePaolo (1981) and by Taylor and Sheppard (1986) in their simpler analytical model. As a matter of fact, eq. (22) (modal melting) and eq. (30) (non-modal melting) turn to eq. (6a) of DePaolo (1981) for $D_0$ = 1 and $D_0$= P = 1, respectively, for any value of D. Eq. (23) turns to an indeterminate for $D_0$= 1 as occurs for z = 0 in DePaolo’s (1981) mathematical treatment. Eq. (39) (modal melting) and eq. (46) (non-modal melting) both change to eq. (11a) of DePaolo for $D_0$ = 1 and for $D_0$ = P = 1, respectively, for any value of D. And eqs. (41) and (46) turn to indeterminates for z = 0 as in DePaolo’s mathematical treatment if $D_0$ = 1 and $D_0$= P = 1, respectively.

In Fig. 1 and Fig. 2, graphical examples show the differences existing between Sr concentration vs. isotopic composition evolution paths of the magmatic liquid calculated by solutions proposed in this study and the respective ‘traditional’ DePaolo’s (1981) or Taylor and Sheppard’s (1986) solutions. The calculations assume modal melting of the wallrock, and the following conditions: Sr concentration in the initial magmatic liquid: 100 ppm; ⁸⁷Sr/⁸⁶Sr isotopic ratio in the initial magmatic liquid: 0.704; Sr concentration in the wallrock which is subject to melting: 500 ppm; ⁸⁷Sr/⁸⁶Sr isotopic ratio in the rock subject to melting: 0.720; Sr bulk partition coefficient between the crystal assemblage that forms from the magmatic liquid and the magmatic liquid: 2; Sr partition coefficient between wallrock and melt in partial melting: 0.5; the mass ratio between assimilation and crystallization (parameter r): 0.15. In Fig. 1, the mass ratio μ between the initial magmatic liquid and the wallrock involved in the melting process is assumed to be 5, whereas in Fig. 2 this ratio is assumed to be 1.

Fig. 1.

Ideal Sr concentration vs. isotopic composition evolution paths in a magmatic liquid evolving by concurrent fractional crystallization and assimilation. The calculation assumes modal melting of the wallrock. The evolution path predicted by DePaolo’s (1981) model and the paths predicted by Mod. 1 and Mod. 2, this work, are shown. The residual liquid mass fraction is in the range 1–0.01, and the wallrock melting degree (Mod. 1 and Mod. 2) in the range 0–87%. See text for further explanations.

Fig. 2.

Ideal Sr concentration vs. isotopic composition evolution paths in a magmatic liquid evolving by concurrent fractional crystallization and assimilation. The calculation assumes modal melting of the wallrock. The evolution path predicted by DePaolo’s (1981) model and the paths predicted by Mod. 1 and Mod. 2, this work, are shown. The residual liquid mass fraction is in the range 1–0.01, and the wallrock melting degree (Mod. 1 and Mod. 2) in the range 0–17%. See text for further explanations.

Both in Fig. 1 and in Fig. 2, are shown the evolution paths of the magmatic liquid predicted by the model of DePaolo (1981) (pink curve) and those predicted by Mod. 1 and Mod. 2 (eq. (22) and eq. (25), blue curve and black curve, respectively). As already extensively shown by Bohrson and Spera (2001), the evolution paths predicted by an AFC process where the concentration in the assimilate changes according to Shaw’s (1970) fractional melting model of the wallrock are often very different from those predicted by DePaolo’s model.

These analytical solutions make it possible to ascertain if a chemical and isotopic data set can be interpreted in terms of a process of concurrent fractional crystallization and wallrock assimilation where the wallrock is assimilated according to Shaw’s fractional melting model, as assumed by Spera and Bohrson (2001) in their energy-constrained AFC model (see eq. (8) in Spera and Bohrson (2001)). This is because in these solutions, along with parameter r, a new parameter, μ, is fundamental in an energy-constrained process. As a matter of fact, this parameter identifies the mass ratio between magmatic liquid—which generates heat by its own crystallization—and wallrock involved in the melting process—which, owing to that amount of heat, melts.

Parameters r and μ are interdependent because they calculate the melting degree of the wallrock at any step (eq. (13)), and the melting degree of the wallrock can never exceed 1. These parameters can be calculated by Bohrson and Spera’s (2001) energy-constrained data, and therefore it is possible to check if the solutions here proposed match the calculations of Bohrson and Spera (2001). If so, these solutions can be confidently used in place of Bohrson and Spera’s calculation.

Conclusions

Analytical solutions are proposed which predict the evolution of the concentration and of the isotopic composition of a trace element in magmatic liquids which evolve by processes of assimilation-fractional crystallization where the concentration of the trace element in the instantaneous assimilate changes according to Shaw’s (1970) model of fractional melting of the wallrock, and the isotopic composition of the trace element in the assimilate is constant. This is the model used by Spera and Bohrson (2001) and Bohrson and Spera (2001) in their energy-constrained model of assimilation and fractional crystallization (EC-AFC).

The solutions here proposed accurately predict the evolution of the concentration, and of the isotopic composition of a trace element within the limits which have been discussed in the text, limits corresponding to simplifications of parts of the respective starting differential equations.

Acknowledgments

This work has been financially supported by Consiglio Nazionale delle Ricerche, Istituto di Geoscienze e Georisorse, Italy. Dr. D. Roccato and Prof. E. M. Piccirillo are gratefully acknowledged for valuable discussions.

References

•  Aitcheson,  S. J. and  Forrest,  A. H. (1994) Quantification of crustal contamination in open magmatic systems. J. Petrol. 35, 461–488. https://doi.org/10.1093/petrology/35.2.461
•  Bohrson,  W. A. and  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. https://doi.org/10.1093/petrology/42.5.1019
•  Bowen,  N. L. (1928) The Evolution of the Igneous Rocks. Oxford Univ. Press, 332 pp.
•  Cavazzini,  G. (1996) Degrees of contamination in magmas evolving by assimilation-fractional crystallization. Geochim. Cosmochim. Acta 60, 2049–2052. https://doi.org/10.1016/0016-7037(96)00089-0
•  Cavazzini,  G. (2000) Modelling the Sr isotopic evolution of melts produced by fractional melting. Geochem. J. 34, 395–406. https://doi.org/10.2343/geochemj.34.395
•  Cavazzini,  G. (2010) Calculations of wallrock melting degree in AFC processes (abstract). Proc. SIMP (It. Soc. Min. Petrol.) Conf. 89th, 133. https://publications.cnr.it/doc/185184
•  DePaolo,  D. J. (1981) Trace element and isotopic effects of combined wallrock assimilation and fractional crystallization. Earth Planet. Sci. Lett. 53, 189–202. https://doi.org/10.1016/0012-821X(81)90153-9
•  DePaolo,  D. J. (1985) Isotopic studies of processes in mafic magma chambers: I. The Kiglapait intrusion, Labrador. J. Petrol. 26, 925–951. https://doi.org/10.1093/petrology/26.4.925
•  Hertogen,  J. and  Gijbels,  R. (1976) Calculation of trace element fractionation during partial melting. Geochim. Cosmochim. Acta 40, 313–322. https://doi.org/10.1016/0016-7037(76)90209-X
•  Marsh,  B. D. (1988) Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization. Contr. Mineral. and Petrol. 99, 277–291. https://doi.org/10.1007/BF00375362
•  Minster,  J. F. and  Allégre,  C. J. (1978) Systematic use of trace elements in igneous processes. Contr. Mineral. and Petrol. 68, 37–52. https://doi.org/10.1007/bf00375445
•  Neumann,  H.,  Mead,  J. and  Vitaliano,  C. J. (1954) Trace element variation during fractional crystallization as calculated from the distribution law. Geochim. Cosmochim. Acta 6, 90–99. https://doi.org/10.1016/0016-7037(54)90018-6
•  Shaw,  D. M. (1970) Trace element fractionation during anatexis. Geochim. Cosmochim. Acta 34, 237–243. https://doi.org/10.1016/0016-7037(70)90009-8
•  Shaw,  D. M. (1978) Trace element behaviour during anatexis in the presence of a fluid phase. Geochim. Cosmochim. Acta 42, 933–943. https://doi.org/10.1016/0016-7037(78)90104-7
•  Shaw,  D. M. (1979) Trace element melting models. Phys. Chem. Earth 11, 577–586. https://doi.org/10.1016/0079-1946(79)90055-7
•  Spera,  F. J. and  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. https://doi.org/10.1093/petrology/42.5.999
•  Taylor,  H. P. (1980) The effects of assimilation of country rocks by magmas on ¹⁸O/¹⁶O and ⁸⁷Sr/⁸⁶Sr systematics in igneous rocks. Earth Planet. Sci. Lett. 47, 243–254. https://doi.org/10.1016/0012-821X(80)90040-0
•  Taylor,  H. P. and  Sheppard,  S. M. F. (1986) Igneous rocks: I. Processes of isotopic fractionation and isotope systematics. Stable Isotopes in High Temperature Geological Processes 8 (Valley, J. W., Taylor, H. P. and O’Neil, J. R., eds.), 227–272, De Gruyter. https://doi.org/10.1515/9781501508936-013
•  Thompson,  A. B.,  Matile,  L. and  Ulmer,  P. (2002) Some thermal constraints on crustal assimilation during fractionation of hydrous, mantle-derived magmas with examples from central Alpine batholiths. J. Petrol. 43, 403–422. https://doi.org/10.1093/petrology/43.3.403