Lecture 13 Box Models and Mixing Theory Origin

Lecture 13: Box Models and Mixing Theory, Origin of Atmosphere and Hydrosphere • Questions – How do geochemists interpret compositional data in terms of mixing of a small number of discrete components? – How do geochemists use the idea of uniform reservoirs and fluxes to understand the cycling of elements through the Earth? – What can such mixing calculations and box models tell us about the history of the Earth’s fluid envelope? • Tools – Conservation of Atoms • Reading: – Albarède, Chapters 2 and 5 1

Mixing and Box Models: Motivation • Much of geochemistry consists of measuring the compositions of natural samples and trying to make sense of the measurements. • What sort of conclusions can usefully be drawn from such data? – What source or sources contributed to the observed samples? – What process or processes acted on the source to produce the observed samples? • The fundamental physical law used to trace matter from sources through processes to products is conservation of atoms (modified by radioactive decays) – Every process is either a differentiation, which takes a uniform source and makes two or more distinct products (e. g. , core formation, partial melting, fractional crystallization…) – Or a mixing process, which combines two or more distinct sources to make a range of products with intermediate compositions or, in the extreme, a uniform complete mixture. 2

Mixing Theory • Geochemistry often tries to model variations in measured composition as the result of mixing of a small number of components [N. B. a different use of component from thermodynamic usage] or end members – This reduces highly multivariate data to a few manageable dimensions – It allows identification of the end members with particular source or fluxes, hence a meaningful interpretation of data – Many geochemical processes are easily understood in terms of mixing or unmixing: • river water + ocean water = mixing • primary liquid – fractionated crystals = unmixing • We will work out the mixing relations for several spaces: – Element-element plots – Element-ratio plots (including elemental and isotope ratios) – Ratio-ratio plots 3

Mixing Theory • Mixing is simplest to see and understand when there are only two end members: Binary mixing • For concreteness, instead of a bunch of general symbols, let’s do all this with two end members, Archean gneiss (a) and basalt (b), with the following compositions: • The same relationships will apply for mixing of major elements as for trace elements. • The same relationships will apply for ratios of major elements, ratios of trace element concentrations, and ratios of isotopes. – I’ll leave it up to you to generalize the specific equations from this lecture, rather than the other way around. 4

• • • Binary Mixing I: element-element This is the simplest case. Binary mixing in concentration space always generates lines. Let mixtures be generated with mass fraction fa of end member a and fb of end member b, such that fa + fb = 1. Then for two species, say Sr and Nd, we have conservation of atoms and mass in the form [Sr]mix = fb[Sr]b + (1 – fb)[Sr]a [Nd]mix = fb[Nd]b + (1 – fb)[Nd]a This can be written [Sr]mix – [Sr]a = fb([Sr]b – [Sr]a) [Nd]mix – [Nd]a = fb([Nd]b – [Nd]a) Dividing these two equations gives the equation of the mixing relationship in ([Sr], [Nd]) space: 5

Binary Mixing I: element-element • This is the equation of a line with slope • Passing through points ([Nd]a, [Sr]a) and ([Nd]b, [Sr]b) 6

Binary Mixing I: element-element • If you know the compositions of the end members, you can solve for fb using the lever rule: • If you don’t know the end members, but only the data, what can you learn from a graph showing a linear correlation? – You can infer that if generated by mixing there are only two end members, otherwise the data would fill a triangle – You can infer that both end members lie on the mixing line, outside the extreme range of the data on both ends if they must have positive amounts (additive mixing) • Because mixing is linear in concentration space, you can use linear least squares analysis to interpret data in any number of dimensions for any number of end members. 7

Binary Mixing II: element-ratio • In geochemistry we very frequently work with ratios, either isotope ratios or ratios of concentrations. – Sometimes a ratio is all you can measure accurately – Sometimes ratios have significance where concentrations are more or less arbitrary (example: during fractionation of olivine from a basalt, [Sr] will change because it is incompatible and the amount of liquid is decreasing, but [Sr]/[Nd] will not) • You might think that the 87 Sr/86 Sr ratio of mixtures could be obtained in the same way as for [Sr]: (87 Sr/86 Sr)mix = fb (87 Sr/86 Sr) b + (1 – fb) (87 Sr/86 Sr) a • You would be wrong! • The isotope ratio of the mixture is going to be weighted by the concentration of Sr in each end member – More generally, the weighting of ratios in the mixture is controlled by the denominator of the ratio, in this case 86 Sr. 8

Binary Mixing II: element-ratio • Let’s do it right: we have [87 Sr]mix = fb [87 Sr] b + (1 – fb) [87 Sr] a [86 Sr]mix = fb [86 Sr] b + (1 – fb) [86 Sr] a • Taking the ratio of these, we have • And substituting [87 Sr] = (87 Sr/86 Sr)[86 Sr] for a and b: • Because differences in isotope ratios are much smaller than differences in concentrations, we can approximate this using [Sr] instead of [86 Sr] as the weighting factors. 9

Binary Mixing II: element-ratio • Now let’s consider plotting the isotope ratio of the mixture against an elemental concentration, perhaps [Nd]. • If we eliminate fb between the mixing equation for (87 Sr/86 Sr) and [Nd], we obtain the following equation (using the approximation Sr ~ 86 Sr): • What curve has general equation Ax +Bxy +Cy + D = 0? 10

Binary Mixing II: element-ratio • In general, element-ratio mixing generates a hyperbola. – The only case in which it is linear is B = [Sr]b–[Sr]a = 0 • The index of curvature r = [Sr]b/[Sr]a tells you “how hyperbolic” the hyperbola is going to be. 11

Binary Mixing III: inverse element-ratio • Although there is nothing special about the element-ratio case A/B vs. B (thus (87 Sr/86 Sr) vs. [Sr] is still hyperbolic), there is an especially useful test for mixing if you plot A/B vs. 1/B (in this case, (87 Sr/86 Sr) vs. 1/[Sr]). – Going back to our hyperbolic equation, if we replace [Nd]mix, [Nd]a and [Nd]b with [Sr]mix, [Sr]a and [Sr]b we have: • Or, , which is a line. 12

Binary Mixing III: inverse element-ratio • Mixing in A/B vs. 1/B space always generates a line. • The value of r = [Sr]b/[Sr]a now controls how hyperbolic the spacing of equal increments of mixing fraction are along the line. Since linear correlation is easy to calculate, it is much easier to test whether data are consistent with binary mixing in this space than in ratio-element space. 13

Binary Mixing IV: ratio-ratio • Our final case is plots of ratios against ratios, whether isotope ratios, trace element ratios, or major element ratios. • For example, let’s do (87 Sr/86 Sr) vs. (143 Nd/144 Nd). • We now have two equations of the same form: Where once again for the particular case of small variations in isotope ratios I have weighted by concentration rather than by the stable denominator isotope. 14

Binary Mixing IV: ratio-ratio • This time eliminating fb between the mixing equations gives • Still a hyperbola. Now term B gives the curvature index r = ([Sr]a/[Sr]b)/([Nd]a/[Nd]b) = ([Sr]a/[Nd]a)/([Sr]b/[Nd]b) 15

Binary Mixing IV: ratio-ratio This looks different from the element-ratio hyperbola because now the spacing of equal increments of mixing fraction is no longer regular. Given an array of ratio -ratio data, you can constrain the curvature parameter as well as the isotope ratios of the end members. 16

Binary Mixing IV: a practical example . 703 (87 Sr/86 Sr) . 704 We’d like to know whether the lower mantle has high 3 He/4 He because it is rich in 3 He or because it is poor in 4 He. Degassing at eruption makes this question very hard to answer. Assuming Loihi represents the plume or lower mantle component, mixing with an upper mantle MORB-like component richest at Mauna Loa, what does this plot say about the (4 He/Sr) ratios of the upper and lower mantles? 17

Box Models • Consider the geochemical cycling of an element such as Carbon. How do we go from observations to understanding? • In reality we see continuous variation of composition everywhere we go, and it takes an infinite amount of data to completely characterize the world…every air sample or rock sample is different…how then to gain understanding of the essentials? 18

Box Models • We simplify by constructing a Box Model: a set of reservoirs that can be considered homogenous without giving away too much useful information. • We consider all the fluxes that move the element of interest from one reservoir to another; we can associate these with processes as shown by the labeled arrows. 19

Box Models • Quantifying a box model then consists of measuring some concentrations and fluxes and solving for the rest. • According to this diagram, is Ca in steady-state in the ocean? – Do inputs equal outputs? 20

Box Models: residence time • The residence time is the mean time that a given object stays in a reservoir. • Consider a barrel that holds 100 stopwatches. Every minute we zero ten stopwatches and drop them in and every minute we randomly pick out 10 stopwatches and read them. 21

Box Models: residence time • Given the amount of an element in a reservoir and the flux out of the reservoir, we can calculate the residence time of the element in the reservoir – t = Mass/Flux (or other appropriate units: volume, moles…) – This is interpreted as the average time that an atom of the element resides in the reservoir before being removed. – If you calculate a residence time using a flux into the reservoir, you are implicitly making an assumption of steadystate! • Examples: Consider the ocean – Mass of ocean water: 1. 4 x 1021 kg – Total river H 2 O flux: 3. 7 x 1016 kg/yr – Residence time of H 2 O in the oceans: 3. 6 x 104 years 22

Box Models: residence time • Now consider some elements in the ocean: Na – Concentration of Na in ocean water: 1. 08 x 10 -2 – Concentration of Na in river water: 6. 3 x 10 -6 – Mass of Na in oceans: 1. 5 x 1019 kg – Flux of Na into oceans: 2. 3 x 1011 kg/yr – Residence time of Na in ocean from input: 6. 6 x 107 years – Residence time based on Na burial in sediment: 2 x 108 yr • Fe – Concentration of Fe in ocean water: 6 x 10 -11 – Concentration of Fe in river water: 4 x 10 -8 – Mass of Fe in oceans: 8. 4 x 1010 kg – Flux of Fe into oceans: 1. 5 x 109 kg – Residence time of Fe in ocean from input: 57 years – Residence time of Fe from burial in sediment: 0. 7 years • Evidently either the ocean is far from steady state or we are missing some inputs and outputs! 23

Box Models: A motivational example Is mantle convection layered? 40 Ar provides geochemistry’s strongest answer: The mass of 40 Ar in the atmosphere, 66 x 1018 g. We know the K content (250 -285 ppm) and age of the Earth. 40 Ar, produced by decay of 40 K, is retained by the Earth. There should be ~150 x 1018 g. We know the 40 Ar content of the upper mantle from the flux at mid-ocean ridges. 60 -80 x 1018 g are missing, requiring an undegassed, high-K lower mantle. 24

Box Models: Working out Details in a simple case • Consider the following model of the earth’s mantle and the evolution of a radioactive decay system: • Reservoir 1 is the lower mantle, with mass M 1 and concentration of parent element c 1 • Reservoir 2 is the upper mantle, with mass M 2 and concentration of parent element c 2 • The flux of mass from reservoir 1 to 2 it F 12 and we will assume the masses of the two reservoirs are constant, so F 21 = F 12 = F. • The mean life of the radioactive element is t. • For any initial condition (c 1, c 2), we want to be able to calculate the concentrations in the two reservoirs at any time t. 25

A Useful Box Model • We write the equations for conservation of mass of the element in each box: • If M 1 and M 2 are constant, we can take them out of the derivatives and recast these expressions in terms of residence times: Where t 1 = M 1/F and t 2 = M 2/F are the residence time in each reservoir and t is still the mean-life of the isotope. Note the similar role of radioactive mean-life and residence time: mixing and decay act in the same way in this system 26

A Useful Box Model • We can also write the system of first-order differential equations as a matrix problem: – If the concentration vector c is defined by – then the rate of change of c is written – Or more compactly as • It is perfectly straightforward to solve this problem numerically, stepping forward each time-step using the derivative above. However, for the case of constant coefficients a useful analytical solution is available. 27

Solution by eigenmode decomposition • We can quickly find the general solution to this problem, and understand its significance, if we think of this as an eigenvector equation: • Why? Because if c is an eigenvector, then we have the uncoupled system – and we know the solution to this problem, – which is to say, if the initial condition is an eigenvector, then it decays (or grows, or oscillates) without changing the ratio c 1/c 2, or the direction in (c 1, c 2) space of the vector. 28

Solution by eigenmode decomposition • Theorem from linear algelbra: Any non-singular 2 x 2 matrix A has two linearly independent eigenvectors c 1*, c 2* • Furthermore, because any two linearly independent vectors form a basis that spans (c 1, c 2) space, we can write any vector c as a linear combination of the two eigenvectors: • And since the system of equations is linear, by the principle of superposition (not the stratigraphic one), we can write the general solution for any initial vector as a linear combination of the solutions for the two eigenvectors (i. e. , the two eigenmodes): – where the constants a and b are determined by the initial condition co = ac 1* + bc 2*. 29

Solution by eigenmode decomposition • For the problem as formulated a few slides ago, the eigenvectors and eigenvalues are: • The first mode is the decay mode. If we start with the initial condition c 1 = c 2, then mixing between the boxes has no effect on concentration. The abundance in the whole system simply decays with the mean life of the radioactivity. • The second mode is the mixing mode. It shows the time constant on which the difference in concentration between the two reservoirs decays. – We can’t start with an initial condition of pure mixing mode unless we allow negative concentration. – Note, since residence times and mean lives are positive, the mixing mode always decays faster than the decay mode 30

A geochemical application • Now let us take the element to be 238 U, t = 6. 45 x 109 years. • Let’s take box 1 to be the lower mantle and box 2 to be the upper mantle with M 1 = 2 M 2. Also, M 2 ~ 1024 kg. • Let’s imagine that slabs are subducted into the lower mantle, so that the flux between reservoirs is the present-day slab flux: – 43000 km of trench x 100 km thick slabs x 50 km/Ma average rate x 3. 5 x 1012 kg/km 3 = 7. 5 x 1020 kg/Ma – This gives residence time in lower mantle 2. 72 Ga, and residence time in upper mantle 1. 36 Ga. • Let us imagine starting from the initial condition that all the 238 U is in the lower mantle (perhaps the upper mantle was totally stripped of U by early continent formation). So we take c 10 = 1, c 20 = 0. 31

A geochemical application • Here is the solution for the parameters and initial conditions we chose. • The decay of the mixing mode is traced by the ratio of concentrations in the two reservoirs. It decays quickly…the eigenvalue is – 1/(800 Ma) • Once decay mode is gone, c 1 ~ c 2 • Then they both decay together with mean-life of radioactive decay. • So this model predicts, at the present day, equal U concentration in upper and lower mantle. – Could change this conclusion by reducing flux between upper and lower mantle (layered convection…slabs not penetrating). 32 – Intuition may have told you that, but now it is quantified.

A geochemical application: more wrinkles • We know, for example, from He isotopes that the upper mantle and the lower mantle are not equal in composition. Do we conclude from our box model that whole mantle convection is impossible? Well, perhaps not… • Evidently the flux from upper mantle to lower mantle is slabs, which may be different in composition from average upper mantle because of processes at ridges and subduction zones. – So we introduce a fractionation factor a 2, which is the ratio of 238 U concentration in the mass flux to reservoir 1 over the concentration in reservoir 2. • Likewise the flux from lower mantle to upper mantle may be plumes, which might be different in composition from average lower mantle, perhaps because of things going on at the coremantle boundary. – So we introduce factor a 1, the ratio of 238 U concentration in the mass flux 1 ->2 over the concentration in reservoir 1. 33

A geochemical application: more wrinkles • Now the governing equations are • And the eigenmodes become 34

A geochemical application: more wrinkles Now let’s take a 1 = 1 and a 2 = 2 because perhaps the incompatible U concentrated into the basaltic part of the slab is drawn from twice the volume of mantle that will be subducted. Now the eigenvector of the decay mode is • So, over geological time as the mixing mode decays, the system approaches a state where the concentration in the lower mantle is twice the concentration in the upper mantle…even though we have whole mantle convection! – One might try adding He and Pb to this box model, and building a coupled system for parents and daughters. 35

Enough math: where did the atmosphere come from? • Primordial atmosphere formed during accretion by outgassing from the magma ocean and from volatile-rich impacts. It was highly reducing (in equilibrium with H 2, no O 2). – This atmosphere was lost by catastrophic impacts; the Earth did not become closed to Xe loss until ~4. 45 Ga. • In the early Archean, a second atmosphere formed by outgassing from the Earth, now closed to heavy volatile loss. – Arguably some further volatile input with the late veneer – Whether strongly reducing (H 2, CH 4, NH 3) or weakly reducing (H 2 O, CO 2, N 2) depends on the timing of core formation, which oxidized the mantle by withdrawing metallic Fe. – No O 2 so no O 3, so solar UV reached surface. – Earth probably looked like Venus, with enough CO 2 and H 2 O in the atmosphere to make a very hot surface and keep H 2 O in form of steam. • But at 1 A. U. the equilibrium temperature would have settled at ~85°C, just cool enough to condense the ocean 36

Enough math: where did the atmosphere come from? Dauphas, 2003 37

Early atmosphere • A reducing atmosphere or a high solar UV flux favors loss of H 2 from the exosphere. How do we trace how much H a planet has lost? – 1. Stable isotopes of hydrogen, which is directly fractionated by preferential loss of 1 H relative to 2 D • Earth: D/H = 1. 5 x 10 -4 • Mars: D/H = 9 x 10– 4 • Venus: D/H = 0. 022 – 2. Stable isotopes of other noble gases, which are fractionated during hydrodynamic escape driven by a large flux of escaping hydrogen: • All the terrestrial planets have 20 Ne/36 Ar ~ 0. 01 x Solar • Since H is lost as molecular or atomic hydrogen, the effect is irreversible oxidation of the planet – Together with core formation, leads to a weakly oxidizing early atmosphere 38

Evolution of ocean and atmosphere • Once (if) temperature gets low enough to condense H 2 O, an ocean forms, and CO 2 begins to dissolve in the ocean. – Dissolved CO 2 can precipitate as carbonate rocks – The difference between Earth and Venus is that Earth’s CO 2 is stored in rocks, while 90 bars of CO 2 make up the Venus atmosphere. • Free oxygen comes later, and is a unique signature of life. • Rise of oxygen driven by anaerobic cyanobacteria, who used photosynthesis to pollute themselves almost out of existence by filling the atmosphere with O 2, clearing the way for more efficient aerobic organisms to take over. 39