RBME coral temperature reconstruction: An evaluation, modifications, and recommendations

Available online at www.sciencedirect.com

ScienceDirect Geochimica et Cosmochimica Acta 154 (2015) 66–80 www.elsevier.com/locate/gca

RBME coral temperature reconstruction: An evaluation, modifications, and recommendations Daniel J. Sinclair ⇑ Institute of Marine and Coastal Sciences, Rutgers University, 71 Dudley Road, New Brunswick, NJ 08901, USA Received 22 August 2014; accepted in revised form 8 January 2015; Available online 24 January 2015

Abstract Rayleigh-Based Multi-Element (RBME) coral temperature reconstruction (Gaetani et al., 2011) is an important new method potentially allowing accurate ocean temperatures to be reconstructed from aragonitic corals without the need for in-situ calibration. RBME is based on a mathematical model which assumes that trace element uptake into coral skeleton is mediated by (1) temperature-dependent (inorganic) partition coefficients, and (2) a temperature-dependent change in the proportion of ions in a calcifying fluid that are precipitated as CaCO3. A time-series of Sr, Mg and Ba measurements can therefore be formulated as a set of simultaneous equations allowing temperatures and model fit-parameters to be solved numerically by leastsquares minimization. Here, I critically review the RBME method, focusing on assumptions underlying the mathematical model, implementation of the numerical scheme, and the accuracy of temperature reconstruction. The equations in the RBME are highly nonlinear, and their number scales with the length of the time-series. These ‘highdimensional nonlinear optimization’ problems are characterized by high computational cost, and the difficulty in finding the true ‘global’ minimum rather than false ‘local’ minima. Tests show that several of the RBME equations are highly sensitive to changing parameters thereby making it difficult for minimization algorithms to search broadly in parameter space. The method is therefore exceedingly slow and has a very strong tendency to find local optima that are close to the initial guess (which must be supplied by the user) making it very prone to ‘operator bias’. Here a new formulation of the RBME system is proposed that allows much faster optimization, and a much broader search of parameter space. Testing of the RBME method reveals a number of problems. Firstly, RBME can produce many solutions to a given set of Sr, Mg and Ba data that are nearly identical in their minimization function (i.e., are equally good) but which result in temperatures that are very different. This makes the accuracy of temperature reconstruction by RBME uncertain. Secondly, tests with an artificial Sr, Mg and Ba dataset reveal that even small perturbations, such as analytical noise or vital effects, are enough to destroy the uniqueness of a true ‘global’ solution to the RBME equations. Such sensitivity means that simplifying assumptions made in the mathematical representation of coral biomineralization cannot be safely ignored. Finally, arbitrary choices in the formulation of the equations for the temperature dependence of Rayleigh Fractionation and minimization function can be critical in determining the outcome of the method. I conclude with a list of recommendations for use of the RBME method and reporting of results so that uncertainties in the reconstructed temperatures are accurately represented. Ó 2015 Elsevier Ltd. All rights reserved.

1. INTRODUCTION ⇑ Address: School of Geography Environment and Earth Sci-

ences, Victoria University of Wellington, Cotton Building, Kelburn Parade, Kelburn, Wellington 6012, New Zealand. Tel.: +64 4 463 9755. E-mail address: [email protected]. http://dx.doi.org/10.1016/j.gca.2015.01.006 0016-7037/Ó 2015 Elsevier Ltd. All rights reserved.

Trace elements incorporated into the skeletons of tropical corals are often strongly correlated with the temperature of the water in which the coral was growing. This has led to the widespread use of skeletal trace elements as

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

paleothermometers. Elements that have been reported showing strong seasonal cycles are Sr, Mg, Ba, B, U, Li (Smith et al., 1979; Beck et al., 1992; Min et al., 1995; Mitsuguchi et al., 1996; Sinclair et al., 1998; Marriott et al., 2004). Of these, Sr/Ca ratios appear to be the most reliable (e.g., Quinn and Sampson, 2002), demonstrating consistently strong correlation with temperature with a similar sensitivity between individuals of the same genus. Thus, Sr is the most common choice of paleothermometer. A compilation of Sr/Ca vs SST calibrations from the literature, however, reveals a significant variation in the slope and intercept of the calibration equations (Table 1). This variation is too large to be explained by analytical error, and cannot be fully accounted for by uncertainties in inter-laboratory instrumental calibration or the in-situ temperature records. The temperature dependence of Sr uptake therefore appears to vary from individual to individual within a genus of coral, and the uncertainty this produces in SST reconstruction is large. For example, if there is no a priori information upon which to select a particular equation, the calibrations in Table 1 would produce a temperature record that is uncertain by ±5 °C (24–34 °C) in absolute temperature, and ±0.75 °C (2.1–3.6 °C) in seasonal temperature range (calculated for a ‘typical’ mid-tropics coral with a Sr/Ca seasonal range of 0.00875–

67

0.00910 mmol/mol e.g., Alibert and McCulloch, 1997). This variability is ascribed to ‘vital effects’ which loosely refers to a number of possible physiological and physicochemical processes affecting temperature calibration slopes, offsets between individual calibrations, non-stationary calibration and distortions to the seasonal cycle (de Villiers et al., 1994; Fallon et al., 2003; Lea, 2003). The best practice for circumventing vital effects is to calibrate an individual coral using several years of in-situ temperature data, and to then apply this calibration to the part of the coral that predates the in-situ record. In the case of a fossil coral where the record does not overlap with instrumental data second-best practice is to generate a site-specific calibration using a living coral of the same species from the same location. The former approach requires that temperature calibration remains stationary with time, while the latter approach depends on the assumption that individuals of the same species respond in a similar way within the same location. Unfortunately, vital effects can produce variations in coral Sr/Ca vs SST calibration between individuals of the same species in the same general location (e.g., de Villiers et al., 1994; Felis et al., 2004), and need not remain constant if external stressors affect the physiology of the coral. This uncertainty particularly impacts records being spliced together from overlaps between fossil corals (e.g.,

Table 1 Compilation of Sr/Ca vs SST calibrations for Porites genus corals. Porites species

Location

Constant

Slope

Citation

lobata lobata lobata lobata lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea lutea unknown unknown unknown unknown unknown unknown unknown unknown

New Caledonia New Caledonia Rabaul Hawaii New Caledonia Java Great Barrier Reef Dampier SW Madagascar SW Madagascar SW Madagascar SW Madagascar SW Madagascar SW Madagascar SW Madagascar Fiji Rarotonga Fiji Rarotonga New Caledonia New Caledonia Great Barrier Reef New Caledonia Northern Red Sea Palmyra Island Fanning Island Christmas Island Japan Great Barrier Reef Vanuatu Panama

10.716 10.479 10.64 10.956 10.331 ± 0.055 10.78 10.73 10.68 10.011 10.399 10.323 10.484 10.319 10.348 12.035 10.65 11.12 11.73 11.07 10.12 ± 0.045 10.524 ± 0.043 10.4 ± 0.016 10.383 10.781 ± 0.1181 11.45 10.78 11.18 10.333 ± 0.078 10.47 ± 0.01 10.51 10.94

0.062 0.062 0.061 0.080 0.0504 ± 0.002 0.066 0.064 0.062 0.037 0.053 0.050 0.054 0.049 0.050 0.117 0.053 0.065 0.055 0.063 0.057 0.067 0.0575 ± 0.0005 0.061 0.0597 ± 0.00501 0.088 0.065 0.079 0.051 ± 0.003 0.0615 ± 0.0004 0.062 0.070

Beck et al. (1992) Beck et al. (1992) (revised) Quinn et al. (2006) de Villiers et al. (1994) Stephans et al. (2004) Gagan et al. (1998) Gagan et al. (1998) Gagan et al. (1998) Zinke et al. (2004) Zinke et al. (2004) Zinke et al. (2004) Zinke et al. (2004) Zinke et al. (2004) Zinke et al. (2004) Zinke et al. (2004) Linsley et al. (2000) Linsley et al. (2000) Linsley et al. (2004) Linsley et al. (2004) Quinn and Sampson (2002) Quinn and Sampson (2002) Marshall and McCulloch (2002) Crowley et al. (1999) Felis et al. (2004) Nurhati et al. (2009) Nurhati et al. (2009) Nurhati et al. (2009) Felis et al. (2009) Alibert and McCulloch (1997) Correge et al. (2004) Smith et al. (1979)

68

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

Cobb et al., 2003) as there is no possibility of independently calculating an individual-specific SST-Sr calibration for specimens. Recently, significant effort has been invested in numerical models of coral biomineralization to try to constrain the impact of biology on paleoenvironmental reconstruction. McConnaughey (1986, 1989a, 1989b) provided a conceptual framework for physicochemical biomineralization which became the basis of a stable isotope model by Adkins et al. (2003) that was then adapted to explore trace elements by Sinclair (1999, 2005a), Sinclair and Risk (2006). The latter examined differential transport of trace elements relative to Ca by the ion transport enzyme Ca-ATPase, and the effect of raising the pH and carbonate ion activity of a calcifying fluid on trace element activity. More recently, a number of authors (Cohen et al., 2006; Gaetani and Cohen, 2006; Gagnon et al., 2007, 2012) explored the implications of Rayleigh Fractionation and temperature dependent partition coefficients on trace element coprecipitation into coral skeleton. These studies explore similar Rayleigh Fractionation concepts developed in the study of other carbonate biominerals (e.g., Elderfield et al., 1996). In 2011, Gaetani et al. (2011) published their breakthrough paper describing Rayleigh-Based Multi Element (RBME) paleotemperature reconstruction. This mathematical method builds upon the measured temperature dependence of element partitioning into inorganic aragonite (Gaetani and Cohen, 2006), and utilizes measurement of Sr, Mg and Ba in multiple subsamples of coral to deconvolve two temperature dependent processes: (1) temperature-dependent partitioning of elements into aragonite and (2) temperature-dependent changes in the degree of Rayleigh Fractionation occurring during skeletal formation. Provided sufficient subsamples are analyzed (in this case N P 4), the temperature for each can be calculated. This is potentially a very powerful technique: Because the biological influence on element uptake and precipitation is isolated and removed by the mathematical scheme, the method will potentially work for any aragonitic coral (deep-sea or tropical) without requiring external temperature calibration. What is more, the technique can accommodate cases where changing coral biology results in trace element ‘drift’ through time. Given the potential for RBME to revolutionize coralbased paleotemperature reconstruction, it is essential that this method be given a careful and independent evaluation. In this paper, I critically examine the RBME method. The evaluation focuses on the validity of important assumptions underlying the technique, practical issues associated with numerical solution of the equations in the system, accuracy and uniqueness of solutions, and the effect of uncertainty, noise and natural environmental variation in the measured trace elements on temperature reconstruction. In this paper I present a reformulation of the RBME equations that address some identified problems in the originally-published RBME scheme, and in Electronic Annexes I include an efficient implementation of the RBME method scripted in the programming language of Igor Proe (a timeseries and mathematics package distributed by Wavemetrics Inc.). I also present a modified version of the RBME method

utilizing just Mg/Ca and Sr/Ca measurements, which eliminates the dependence of the method on Ba/Ca – an element affected by significant natural environmental variation. Finally, I provide a set of recommendations for application of the RBME method and reporting of results that addresses some of the issues discussed in this paper. 2. BACKGROUND In the following sections, I present a brief overview summary of the RBME method as outlined and implemented in Gaetani et al. (2011). In Section 3 I present my analysis of this methodology. 2.1. Coral biology and Rayleigh Fractionation Corals are assumed to precipitate CaCO3 from a small membrane-isolated pocket of extracellular calcifying fluid (ECF) (McConnaughey, 1986, 1989a, 1989b). This fluid is derived from seawater, and its ion content is modified by the coral to raise the supersaturation state with respect to CaCO3. RBME makes no explicit assumption about how this is achieved, although use of Ca-ATPase to exchange 2H+ ions in the ECF for one Ca2+ ion from outside the ECF is speculated to be an essential step in the mechanism (Gaetani and Cohen, 2006). The elevated saturation state of the ECF causes CaCO3 to precipitate spontaneously, and this precipitation continues until some proportion of the Ca in the ECF has precipitated, whereupon the coral replaces the ECF with fresh seawater and another ‘batch’ of ion-modification-then-precipitation begins. When aragonite precipitates from solution, trace elements either preferentially incorporate into the crystal phase (e.g., Sr and Ba) or the solution phase (e.g., Mg) relative to Ca. In a closed system, this results in Rayleigh Fractionation where Mg/Ca in solution increases, while the Sr/Ca and Ba/ Ca decrease as precipitation progresses. Late-precipitated aragonite therefore tends to be richer in Mg and depleted in Sr and Ba than early-precipitated aragonite. If the proportion of Ca precipitated from the ECF during each batch of replacement/modification/precipitation increases, then the average Ba and Sr content of new coral skeleton will decrease, and the Mg content will increase. 2.2. RBME method – equations The RBME method is based on two temperature dependent processes that change trace element incorporation into skeletal aragonite: (1) The temperature-dependence of the inorganic partition coefficients for Ca, Sr, Mg and Ba (Kinsman and Holland, 1969; Dietzel et al., 2004; Gaetani and Cohen, 2006). (2) An assumed temperature-dependent change in the proportion of Ca2+ precipitated from the ECF during each batch of precipitation. Together these two processes can be encapsulated in a system of nonlinear simultaneous equations which describe

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

the relationship between temperature and each trace element measured within a subsample. In the following summary, Mg is used as an example, but the form of the equations is essentially the same for Sr and Ba. The subscript ‘i’ denotes parameters that are different for each individual subsample of coral (e.g., i = 1–100 for a trace element record with 100 data points). The subscript ‘o’ denotes parameters applicable at the beginning of each batch of replacement/modification/precipitation (before precipitation has commenced). Partitioning of Mg between solution and the growing crystal is given by the Nernst (single element) partition coefficient: DMg;i ¼

½Mgaragonite;i ½MgECF;i

decreases with temperature (Gaetani and Cohen, 2006) so increasing Rayleigh Fractionation must be responsible. This is shown in a graph of FL vs T using the parameters for Eq. (5) published in Gaetani et al. (2011) supporting material (Fig. 1). The physical mechanism responsible for this increase in the proportion of CaCO3 precipitated during each cycle of replacement/modification/precipitation is unknown, but Gaetani et al. speculate that it is the increase in CaCO3 precipitation rate with temperature (the implication being that the rate of ECF replenishment remains constant with time). Eqs. (1)–(5) can essentially be combined into one equation (although in practice this is complex and is not necessary for the implementation of a numerical solution):

ð1Þ 1.002

The temperature sensitivity of the partition coefficients were measured by Gaetani and Cohen (2006) using the protocols of Kinsman and Holland (1969). Data are modeled as exponential equations (or a linear equation in the case of Sr): K K1þ T2 i

ð2Þ

and K

DCa;i ¼ e

K3þ T4

ð3Þ

i

where: Ti is the temperature (in Kelvin) represented by subsample ‘i’, and K1–K4 are (known) empirical fit coefficients published in Gaetani and Cohen (2006). The effect of Rayleigh Fractionation on trace element partitioning is described by the Rayleigh Equation:   D     1  FLi Mg;i Mg Mg   ¼ ð4Þ Ca i Ca 0;i 1  FLDCa;i i

Mg

Where: Ca 0;i is the Mg/Ca ratio in the initial ECF of each batch, after the coral has modified its composition, but before precipitation has begun. FLi is the “mass fraction of initial solution remaining after precipitation has ended” (Gaetani et al., 2011). This essentially defines how much of the Ca2+ and CO2 3 in solution have precipitated by the time the coral replenishes the ECF. A more detailed discussion of this parameter is given below in Section 3.2. FLi itself is assumed to be temperature dependent with a polynomial form: FLi ¼ C 1 þ

C2 þ C 3 T i þ C 4 T 2i Ti

ð5Þ

where: C1–C4 are (unknown) parameters. Eq. (5) is completely arbitrary. The justification given in Gaetani et al. (2011) is: “We found that the algorithm does a better job of avoiding local minima for Eq. (7) when FL is assumed to be a polynomial function of T” (p. 1922, emphasis added). In essence, Eq. (5) is a ‘black box’ for which there is no theoretical basis. The only thing known for sure about FL is that it must decrease (implying greater proportion of CaCO3 precipitated) with temperature because Mg/Ca increases with temperature in corals. This cannot be driven KdMg because it

Published Alternative 1 Alternative 2 Alternative 3 Alternative 4

1.000 Value of Function F

DMg;i ¼ e

69

0.998

0.996

0.994 20

22

24

26

28

30

Temperature (°C)

Fig. 1. A selection of alternative Parameter Sets for the FL equation. FL is defined by 4 parameters (C1–C4) which are modified during RBME fitting. The parameters generating these 5 curves are presented in Table 3. FL is sensibly constrained to between 1.000 and 0.996. Outside these values, the RBME equations become poorly behaved (see Section 3.2). The FL function for the example parameter set published in Gaetani et al. (2011) is graphed in red (solid line). Alternative 1 was a viable solution for the Madrepora dataset (see Section 4.2) found during RBME fitting. These parameters are within 0.06% of the published parameters, but already the behavior of FL is significantly different from the published function. Sensitivity testing shows that randomly perturbing values for C1–C4 by more than 0.1% results in FL behavior that can fall significantly outside the tolerance range of 1.000–0.996. Alternatives 2 and 3 were found by randomly perturbing Alternative 1 and selecting the few FL curves (out of many thousands of trials) that fell within the tolerance range. Note that the perturbations of C1–C4 are several orders of magnitude larger than the 0.1% sensitivity described above, so these alternatives represent distant, isolated ‘islands’ of viable parameters that would be difficult for a search heuristic to find when starting from the published results. Alternative 4 represents an even larger perturbation and was found by drawing an arbitrary curve passing through the tolerance region for FL, then fitting the FL equation to estimate C1–C4. This demonstrates that there exist potentially many arbitrary, isolated sets of C1–C4 that produce viable function behavior for the FL equation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

70

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

    ! Mg Mg ¼ f T i ; C1; C2; C3; C4; Ca i Ca o;i

ð6Þ

Analogous equations can be defined for Sr/Ca and Ba/Ca (noting that KdSr is modeled with a linear, rather than exponential temperature dependence in Gaetani and Cohen, 2006):     ! Sr Sr ð7Þ ¼ f T i ; C1 ; C2 ; C3 ; C4 ; Ca i Ca o;i !     Ba Ba ¼ f T i ; C1 ; C2 ; C3 ; C4 ; ð8Þ Ca i Ca o;i If N = the number of data points in a coral time series, then Eqs. (6)–(8) represent a system of 3N equations in 4N + 4 unknowns. Since the number of equations can never exceed the number of unknowns, this system of equations can never be solved. To solve these equations, one must make the simplifying    Sr  Ba assumption (Assumption 1) that Mg ; and Ca are Ca o Ca o o the same for all i subsamples. Assumption 1 essentially means is that for each batch of replacement/modification/ precipitation, the coral modifies the initial seawater consistently each time so that each new batch of ECF starts out with precisely the same TE/Ca ratio (but not necessarily absolute concentration of elements). Assumption 1 reduces the system to 3N equations in N + 7 unknowns. By the time N = 4, there are 12 equations in 11 unknowns and each of the unknown constants, and each of the temperatures, can be calculated. Measuring more subsamples allows the estimates of the unknowns to be refined using least-squares minimization routines. It should be noted that several other assumptions are implicit in the mathematical treatment of biomineralization built into RBME. Firstly, since trace elements in the ECF are sourced ultimately from seawater, this implies (Assumption 2) that trace elements do not vary in the natural environment. Secondly (Assumption 3) temperature is assumed to be the only process influencing FL such that there is a unique value of FL (hence TE/Ca) associated with any specific temperature. The merits of Assumptions 1–3 are discussed in detail in Sections 5.1–5.3.

Where: subscripts ‘model’ and ‘meas’ refer to the modeled and measured (observed) TE/Ca values for each data point in the trace element time series, and SFSr, SFMg and SFBa represent scaling factors. In the scheme published by Gaetani et al. (2011) these scaling factors are, respectively,  Sr    Ba ; Mg and Ca . Ca meas;i Ca meas;i meas;i The heuristic favored by Gaetani et al. (2011) is the Nelder–Mead Simplex algorithm (Lagarias et al., 1998) which systematically tests N + 1 trial solutions and modifies the poorest fit each iteration to narrow in on a best solution. 3. EVALUATING THE IMPLEMENTATION OF RBME 3.1. Numerical optimization Numerical minimization of Eq. (9) can be challenging. The equations involved are highly nonlinear (containing exponential and polynomial terms), and the dimensionality of the problem scales with the number of observations made. For a moderate-length trace-element record (e.g., 10 years of monthly-resolution data) Eq. (9) has 127dimensions! Much has been written about high-dimensionality nonlinear optimization (Weise, 2013). The primary concerns are the cost (in computational time) to find a solution, and the presence of local optima which represent essentially ‘false positives’. Different search heuristics perform differently for these two criteria, and commonly there is a trade-off between the two – algorithms that are efficient (fast to converge on a solution) are often prone to finding local optima, while algorithms that are better at avoiding local optima are often very slow. Early tests using the equations and numerical scheme published in Gaetani et al. (2011) revealed that the RBME method suffers from both problems. Computation takes an exceedingly long time (sometimes many days to converge on a single solution) and the solution found was not the global minimum. One of the big problems is that the equation set proposed for RBME is ‘poorly behaved’ – exhibiting extreme sensitivity to input parameters that makes aggressive optimization impossible, while increasing the number of local minima within the solution space.

2.3. RBME method – numerical solution

3.2. Poorly behaved equations

Eqs. (1)–(5), (and their analogous counterparts for Sr and Ba), cannot be solved explicitly and must instead be solved numerically. This is achieved by starting with initial guesses    Sr  Ba for all unknowns (Ti, C1–C4, Mg ; and Ca ). From Ca o Ca o o these initial guesses, a ‘model’ value for each Mg  Sr  Ba ; and Ca are calculated from Eqs. (6)–(8). A Ca i Ca i i search heuristic is then used to minimize the least-squares function:

At the most basic level, numerical optimization algorithms work by generating new – hopefully better – solutions to a function based on information provided by previous solutions (Weise, 2013). This works best when the solution space is smoothly sloping so that previous solutions provide information about the gradient and thus point in the direction of maximum or minimum value. However, if the solution space is chaotic either because of noise, the presence of multiple local optimums, extreme gradients, or regions where the function returns nonnumerical results (e.g., division by zero), then previous solutions may offer no meaningful information upon which to base the next generation of new solutions, causing search algorithms to fail. One way around this is to reduce the distance between previous and new solutions until a scale is

0   !2 !2  Sr  Mg   Sr X  Ca  Mg Ca model;i Ca model;i Ca meas;i meas;i @ Min ¼ þ SF Sr SF Mg i !2 1 Ba Ba  Ca meas;i Ca model;i A ð9Þ þ SF Ba

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

reached where the solution space changes smoothly. However, this can dramatically increase the computation time, reducing the region of solution space that can be searched and increasing the likelihood of finding local optima rather than the global optimum. The equations presented in Gaetani et al. (2011) demonstrate extreme sensitivity to several of the input parameters, and this can be traced to the parameter ‘FL’. FLi ¼ C 1 þ

C2 þ C 3 T i þ C 4 T 2i Ti

ð5Þ

FL quantifies the degree of calcification and is defined as “the mass fraction of the initial fluid remaining at the end of the experiment” (‘experiment’ in this context being each batch of replenishment/modification/precipitation by the coral). This definition is common for melt crystallization, but is less intuitive for precipitation from aqueous solution where the vast majority of fluid mass is water. The upper bound on FL is 1.000 (since values greater than unity imply that CaCO3 is dissolving not precipitating). The lower bound is set by the mass of DIC in solution that can precipitate as CaCO3. In a calcifying fluid derived from seawater with a pH raised to 9.0 (e.g., Al-Horani et al., 2003), the DIC could be up to 20,000 lmol/kg (estimated using PHREEQ C for solution in equilibrium with atmospheric CO2 – see Sinclair and Risk, 2006). If this all precipitated as CaCO3, 1 kg of seawater would lose approximately 0.002 kg. Thus the mass fraction of fluid remaining after calcification has finished would be 0.998. FL is therefore only meaningful for a very narrow range of values (0.998–1.000) but there is nothing in the numerical system that constrains FL to this tolerance range. This means that optimization algorithms can very easily stray into values where FL that has no physical meaning, and/ or that result in erratic function behavior. C1–C4 in Eq. (5) are unknowns that are modified by the search algorithm. Using typical values (Gaetani et al., 2011 supplementary material), the 4 terms of Eq. (5) can be calculated (Table 2): Note that the final value of FL is a delicate balance between positive and negative terms that are 4000– 14,000 times larger than the tolerated range of FL. Sensitivity analysis shows that changing C1 or C3 by just 0.04& (or C2 or C4 by 0.1&) will be enough to shift FL outside the tolerances. Outside of the tolerance range, the equations start to behave poorly. This derives from Eq. (4) (and the analogs for other elements). Taking Sr/Ca as an example:   D     1  FLi Sr;i Sr Sr   ¼ ð10Þ Ca i Ca 0;i 1  FLDCa;i i

Typical values for DSr and DCa are 2022 and 1729, respectively (Gaetani et al., 2011 supplementary material). For  Sr   Sr  ¼ Ca to better than 1 part values of FL below 0.996, Ca i 0 in 1000. This results in large areas of the solution space that have essentially a flat gradient (giving the algorithm no information that can be used to generate better fits to the Sr/Ca data). However, this flat solution region changes suddenly into extreme sensitivity to FL: between FL = 1.000  Sr and FL = 1.050, Ca increases by a factor of over i

71

Table 2 Calculation of FL using typical parameters and temperature (Ti) = 298.15 K. Constant

Valuea

Term

Value

C1 C2 C3 C4

27.7873886 2614.50620 0.09137820 0.0001037778

C1 C2  Ti C3  Ti C4  Ti2 FLb Tolerancec

27.787389 8.769097 27.244410 9.225164 0.999045 0.998–1.000

a From example calculation presented in Gaetani et al. (2011), supplementary material, p. A1. Note  in Gaetani et al. the values for terms C3 and C2 have been accidentally inverted. b FL = sum of the terms. c See text for calculation of tolerance range.

1,000,000 meaning tiny changes in FL made by the algorithm could result in large and unpredictable changes in Sr/Ca. The only way that an algorithm can optimize equations this sensitive is to begin with an initial guess that is close to the ‘correct’ solution, and drastically restrict incremental changes in C1–C4 each iteration. However, by doing this, the algorithm will be very slow to find optima that are even slightly different from the initial guess (a matter made worse by the low efficiency of the search algorithm preferred by Gaetani et al.), and may not be able to find optima that are separated by greater distances in solution space. This means that it is hard to conduct a broad ‘survey’ of the solution space, and the numerical system is strongly inclined to find optima that are close to the initial guess supplied to the algorithm. This would be fine if realistic values for FL were defined by only a very narrow range of values for C1–C4. Unfortunately, there are multiple, widely separated sets of C1–C4 that produce FL within the region of 1.0. The terms in C1 and C3 act as a matched pair (likewise for C2 and C4): large in value but opposite in sign. For any value of C1, a value of C3 can be determined that satisfies this balance. Thus there are multiple realistic solutions to the FL equation (local optimums) that are distantly separated in parameter space (e.g., see Fig. 1 and Table 3). While some of these might be connected by a continuous gradient (a narrow ‘ridgeline’ in multidimensional parameter space), the need for an algorithm to make tiny incremental steps to avoid catastrophic changes in FL ensures that finding these different solutions can be prohibitively costly in computational terms (i.e., slow). 3.3. Modification of the RBME method The poor behavior of the RBME method can be partially addressed by modifying the defining equations and the numerical method used to optimize them. I make the following changes to the RBME method: (1) I replace ‘FL’ (the mass fraction of solution remaining after precipitation has ended) with the parameter ‘P’ (defined as the proportion of Ca remaining in

72

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

Table 3 A selection of alternative parameter sets for function FL graphed in Fig. 1. Alternative 3b

Alternative 4c

27.7874 2614.5062 0.0913782 0.000103778

27.8000 2613.2900 0.0914148 0.000103716

28.5065 2706.7100 0.0928070 0.000103921

28.5770 2694.1000 0.0932010 0.000103990

29.7900 3038.5000 0.0900390 0.000092749

A viable solution found during routine fitting. See caption for Fig. 1 for details. A random perturbation of Alternative 1. See caption for Fig. 1 for details. An arbitrary parameter set found by fitting the FL equation through an arbitrarily drawn curve. See caption for Fig. 1 for details.

solution after precipitation has ended). In the context of precipitation from aqueous solution, P is more intuitive than FL. However, the primary advantage of this substitution is that P is sensibly constrained between 1 and 0, which makes it possible to design a temperature-dependence equation that is naturally bounded by these values (see point 5). (2) I replace the Nernst single-element partition coefficients (DSr, DMg, DBa, DCa) with the ‘Distribution Coefficients’ (KdSr, KdMg, KdBa) defined as:     TE TE ð11Þ KdTE ¼ DTE =DCa ¼ Ca Crystal Ca Solution

 Pi ¼

1þe

C 1 T i eC 2

 !1 ð13Þ

This equation represents a ve sloping function decreasing from 1 to 0 via a sigmoid-shape (see Fig. 3). Parameter C1 corresponds to the mid-point of the sigmoid (the T where P = 0.5). Parameter C2 controls how steep the sigmoid transition is. The equation can approximate a ve sloping straight line, as well as decreasing and increasing curvature, and

(a + b/T)

1.6

The formal derivation of this expression is presented in Electronic Annex 1. (4) The temperature dependence of the three KdTE parameters was calculated from the original Ca, Sr, Mg and Ba data presented in Gaetani and Cohen (2006). Kd was calculated for each precipitation experiment, and graphed as a function of temperature. Exponential equations of the form K2

KdTE ¼ eK 1 þ T were then fitted to this data (see Fig. 2). (5) The temperature dependence of P is given by:

Kd = e

1.4

a = -12.9 ± 0.6 b = 1860 ± 180

-3

Kd is the more commonly used coefficient in cases where there is kinetic/thermodynamic competition between minor ions (Sr, Mg, Ba) for ion sites in a major phase (CaCO3). Thus, Kd has historically been used to describe trace element uptake when a mineral precipitates from aqueous solution (as is the case for corals). This substitution has one major advantage. Since DSr, DBa and DCa are of similar magnitudes (Gaetani et al., 2011), KdSr and KdBa are close to 1.0. This means that equations where KdSr and KdBa appear as exponential terms (see point 3) do not display the extreme sensitivity that characterizes Eq. (4) (see Section 3.2). (3) Having made substitutions (1) and (2), the Rayleigh Fractionation equation describing the TE/Ca of the batch of aragonite precipitated can then be reformulated as (using Sr as an example):       1  P KdSr;i i Sr Sr ð12Þ ¼ Ca i Ca 0;i ð1  P i Þ

KdMg (x10 )

c

Alternative 2b

1.2 1.0 0.8

Mg

0.6 1.2 a = -1.86 ± 0.37 b = 600 ± 116

KdSr

a b

Alternative 1a

1.1 1.0 0.9

Sr

2.5 KdBa

C1 C2 C3 C4

Published

a = -8.9 ± 1.0 b = 2890 ± 300

2.0 1.5 1.0

Ba

0.5 300 320 340 Temperature (K)

Fig. 2. Calculating the temperature dependence of Kd from the original temperature partitioning experiments presented by Gaetani and Cohen (2006). Partitioning data for Ca, Mg, Sr and Ba were converted into Mg/Ca, Sr/Ca and Ba/Ca data for crystals and solutions, and the Kds calculated according to Eq. (11). An exponential equation of the form Kd = e(a+b/T) was then fitted to the data. The fit parameters and their 95% confidence intervals are presented in the figure.

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

0.8

1.0000 0.9998

0.6

0.9996 0.4

0.9994 FL Function using C1-C4 from Gaetani et al. (2011)

0.9992 0.9990

0.2

Value of function for P

Value of function for FL

1.0002

P Function using C1 = 20.2 °C and C2 = 1.24

0.9988

0.0 10

15

20 Temperature (°C)

25

30

Fig. 3. The original function (FL) and the new function proposed here (P) for the temperature dependent change in amount of Ca precipitated during each cycle of precipitation. See text for definition of FL and P. Parameters for generating FL (C1–C4) are those published in Gaetani et al. (2011) – see Table 1. Note that the proposed equation for P can reproduce a similar shape of temperature dependency to the original function ‘FL’ proposed by Gaetani et al. (2011).

thus can reproduce similar T-dependence shapes to the equation for FL used by Gaetani et al. (2011) (Fig. 3) As with FL, the form of this equation is arbitrary; however, it has been designed to have the following important property: C1 and C2 can take any value (from 1 to +1), and P will always be bounded by +1 and 0. This means that an optimization algorithm can use much larger increments for the unknown parameters (C1 and C2) without risk of producing unrealistic values of P. Combined with the better behavior of the redefined Rayleigh equation (see step 3), this allows the use of far more aggressive optimization algorithms, and a faster (hence wider) search of solution space for minima.Although the sigmoid form of the P equation (where P approaches 1.0 and 0.0 asymptotically) is arbitrary, a case can be made that this behavior is plausible. Coral calcification drops off as temperatures declines (Marshall and Clode, 2004). This would equate to a P value asymptotically approaching 1.0 (zero Ca precipitation) as T drops. Conversely, P decreases (more Ca precipitation) as temperatures increases, but as Ca is removed from solution the supersaturation state of the fluid decreases and calcification rate drops off. Thus full calcification (P = 0) would be expected to be approached asymptotically. This is, however, a simplification and there are multiple factors that might control the shape of the P vs T curve. (6) The scaling factors (SFSr, SFMg and SFBa) used in the minimization function (Eq. 9) have been changed. Gaetani et al. scale the offset between observed and modeled TE/Ca to the observed TE/Ca. Although not obvious, this subtly biases the optimization algorithm toward solutions that offer the best fit to the element with the largest seasonal amplitude (in this case Mg/Ca). This can be rationalized as follows: In a typical tropical coral (e.g., Sinclair et al., 1998), the seasonal amplitude of Sr/Ca spans 7.9–

73

8.4 mmol/mol – an amplitude of 6% relative to the average Sr/Ca. In contrast, Mg/Ca ranges from 3 to 4.3 mmol/mol – an amplitude of 36% relative to the average Mg/Ca. A modeled solution that erroneously results in a doubling of the seasonal amplitude for Sr will produce values that are at most 3% deviated from the observed values. A modeled solution that erroneously results in a doubling of the seasonal amplitude for Mg will produce values that are up to 18% deviated from the observed values. Thus, P the ‘cost’ (as (meas  obs)2) of doubling the Mg amplitude will be greater than the cost of doubling the Sr amplitude, and the algorithm would favor solutions that best fit the Mg/Ca seasonality. To counter this bias toward Mg/Ca, I have scaled the offset between observed and modeled TE/Ca to the seasonal amplitude observed for each element (defined as maximum  minimum). The effect of this on the solutions that the algorithm finds is subtle, since the best solution is still one where all seasonal amplitudes are correct. However, in cases where there is noise in the element profiles, my minimization scheme tends to produce better fits for all 3 elements, at the expense of slightly poorer fits to the Mg/Ca profile. (7) I use the ‘Quasi-Newton’ optimization algorithm (Dennis and Schnabel, 1996) which is natively implemented in Igor Proe instead of the slower Nelder– Mead Simplex algorithm. The use of this more aggressive algorithm is made possible in part by the modifications listed in points 1–5, which result in more robust behavior of the numerical scheme.

4. TESTING THE RBME METHOD 4.1. Numerical optimization strategy The Quasi-Newton and Nelder–Mead Simplex algorithms are both categorized as ‘gradient-climbers’. The Nelder–Mead Simplex is better able to ‘step over’ local minima to find a global minimum; however, it is regarded as an obsolete search heuristic because of its inefficiency and tendency to get stuck in non-convergent oscillation (McKinnon, 1999). The Quasi-Newton method, in contrast, is extremely fast but tends to find local minima, necessitating a strategy of multiple restarts (up to several-thousand perturbed initial guesses). Tests on a laptop computer (2.2 GHz Intel i7 processor) showed that the Nelder–Mead simplex took several days to find a single solution, and needed occasional perturbation to break out of non-convergent states. The Quasi-Newton method was able to find a single optimum within 10 s of seconds, and – with multiple perturbed initial guesses – an excellent solution within a few hours. This method has an additional advantage: the redesigned equations (Section 3.3) are very tolerant of large perturbations of the initial guess, allowing a broad search of the solution space. Occasionally ‘rogue’ solutions (good minima based on significantly different fit parameters) would be found. By iteratively using these as an initial guess for a new cycle

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

of optimization, it was possible to survey different ‘islands’ of solutions, further broadening the search. Running several searches concurrently on different threads in a multicore processor allowed a full spectrum of results could be generated within 2 days. While the Quasi-Newton method is faster and able to search widely in parameter space, the Nelder–Mead Simplex was sometimes better at refining a new solution (especially for noisy datasets – see Section 4.3). In cases like this, a hybrid approach was sometimes most effective, using the Quasi-Newton method to search widely, and the Nelder– Mead Simplex to refine a solution.

minimization function values, but which return different temperature profiles, and (2) it was possible to find excellent solutions where temperature averaged 25 °C, despite the fact that this is a long way from the known temperature in which the coral was growing. A test run using the full original RBME method also converged on a 25 °C solution (Fig. EA 1), indicating that this is not a consequence of modifications made to the RBME method, but rather the initial guess supplied to the RBME algorithm. A direct comparison of solutions generated here with that published by Gaetani et al. (2011) is not strictly valid because different optimization criteria were used in the algorithms (see modification 7 above). However, when measured by the new minimization function, a significant number of the solutions found were significantly better than the solution published by Gaetani et al. (2011) (Fig. 5). When measured by the minimization function originally used by Gaetani et al. many of the solutions (including those from the 25 °C initial guess) were very close to the Gaetani et al. solution (Fig. 4). One of the best solutions found (Fig. 4b) is indistinguishable from Gaetani et al.’s when evaluated using their original optimization criterion, and significantly better when evaluated using the new optimization criterion (due to a notably improved fit to Sr/Ca). This overall best solution returns a temperature profile that averages 7 °C (the average T of the initial guess), but with a significantly larger seasonal amplitude (3.4 °C) than both Gaetani et al.’s solution, and the in-situ temperature record, which vary by 1.1 °C and 0.7 °C, respectively.

4.2. Degenerate solutions and the influence of the initial guess Initial testing of the modified RBME system was undertaken using the dataset for Madrepora published in Gaetani et al. (2011). Due to an oversight, a starting temperature profile of a flat 25 °C was used as an initial guess. A suite of simulations was run (using both the original published RBME method and my modified version) until it was realized that Madrepora is a deep water coral, and that the specimen analyzed by Gaetani et al. grew in a near-constant 7 °C water. A second set of tests was therefore undertaken with a starting temperature of 7 °C. Results from these initial tests are presented in Figs. 4 and 5. The erroneous 25 °C initial guess proved instructive, returning a suite of RBME temperature solutions averaging near 25 °C. Two observations can be taken from these tests: (1) there exist a wide suite of solutions with similar

New RBME, Best 7°C Solution Using New Minimization Criterion

-3

10.6 10.4 10.2 10.0 Mg/Ca

-3

Mg/Ca (mol/mol x10 )

9.8 3.8 3.6 3.4 3.2

Observation Gaetani et al. (2011)

3.0 Ba/Ca

-6

11.0

Gaetani Solution New Min Function = 3.53 Old Min Function = 0.179 New Solution New Min Function = 2.73 Old Min Function = 0.188

10.5 10.0 9.5 9.0 0

20

40

60

Point in Time Series

80

100

Temperature (°C)

4 10.8

Sr/Ca

-3

Sr/Ca

Sr/Ca (mol/mol x10 )

10.8

6

Sr/Ca (mol/mol x10 )

4

8

10.6 10.4 10.2 10.0 9.8 3.8

Mg/Ca

-3

6

10

New RBME, Best 25°C Solution Using New Minimization Criterion

Mg/Ca (mol/mol x10 )

8

B

Temperature

3.6 3.4 3.2

Observation Gaetani et al. (2011)

3.0 11.0

Gaetani Solution New Min Function = 3.53 Old Min Function = 0.179

Ba/Ca

-6

10

New RBME, Best 7°C Solution Using Old Minimization Criterion

12 Temperature (°C)

A

Temperature

Ba/Ca (mol/mol x10 )

-6

Ba/Ca (mol/mol x10 )

-3

Mg/Ca (mol/mol x10 )

-3

Sr/Ca (mol/mol x10 )

Temperature (°C)

12

Ba/Ca (mol/mol x10 )

74

New Solution New Min Function = 2.85 Old Min Function = 0.179

10.5 10.0 9.5 9.0 0

20

40

60

Point in Time Series

80

100

30

C

Temperature

25 20 15 8.0 7.0 10.8

Sr/Ca

10.6 10.4 10.2 10.0 9.8 3.8

Mg/Ca

3.6 3.4 3.2

Observation Gaetani et al. (2011)

3.0 11.0

Gaetani Solution New Min Function = 3.53 Old Min Function = 0.179

Ba/Ca

New Solution New Min Function = 2.98 Old Min Function = 0.195

10.5 10.0 9.5 9.0 0

20

40

60

80

100

Point in Time Series

Fig. 4. 3 New solutions to the Madrepora dataset found using the Modified RBME system. These solutions all result in an excellent match between the observed and modeled trace element profiles, yet produce temperatures that are not realistic, and which do not match the profile published by Gaetani et al. (2011). Note that the values for both the original minimization criterion and the new minimization criterion are presented in the legend. (A) The best fit found for an initial guess of 7 °C according to the new minimization criterion. (B) The best fit found for an initial guess of 7 °C according to the original minimization criterion published in Gaetani et al. (2011). Note that the new solution is – within uncertainty – just as good as the originally published solution when evaluated by the original minimization criterion while being significantly better than the originally published solution when evaluated by the new minimization criterion. (C) The best fit found for an initial guess of 25 °C. Note that this solution is very nearly identical to the one published by Gaetani et al. (2011), except that it fits the Sr data significantly better.

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

75

Fig. 5. Temperature solutions that have similar or better minimization functions (using the new criterion) than the solution published by Gaetani et al. (2011). Note that these represent trials beginning with average 25 °C and 7 °C initial guesses only, so solutions at other average temperatures may also exist. Note also that the solution published by Gaetani et al. (2011) was not found – likely due to a combination of the alternative minimization criterion and the new formulation of the RBME equations.

This degeneracy of solutions implies that there are too many degrees of freedom within the RBME equations (even my modified scheme which contains 2 fewer degrees of freedom than the original) to allow it to produce an unequivocal ‘best’ solution. For example, a suite of multiple solutions may be possible where a steeper slope of the P (or FL) vs T function is balanced by a decreased seasonal T cycle, resulting in an equivalent DP. None of the solutions found to the Madrepora dataset match the one published by Gaetani et al. (2011). One might argue that the change from a polynomial function for FL (Eq. 5) to a sigmoid function for P (Eq. 13) prevents the modified RBME system finding the ‘correct’ solution. However, this raises a critical point: both function forms are arbitrary. If correct functioning of RBME is dependent on a polynomial rather than a sigmoid equation, then this must be backed up with theory, as these equations are just two of many possible functions that could have been chosen. 4.3. Artificial test case and effect of noise Corals are imperfect recorders of climate (e.g., Sinclair et al., 1998, 2011) and trace element profiles will therefore include some variation that is not driven by temperature. Analytical methodology will also introduce noise to any measured trace element profiles. To test the effect of noise on the reconstruction of temperature, I generated a synthetic solution to the RBME equations. This was done by starting with an artificial temperature profile generated from a sinusoid function (see Table 4), and scaled to a temperature range appropriate for a tropical coral (average 25 °C, range 23–28 °C). Typical parameters were then used in the RBME equations (Table 4) to generate profiles for Sr/Ca, Mg/Ca and Ba/Ca. The original temperature profile and equation parameters thus represent a ‘perfect’ solution

to the RBME system: minimization should find this solution and return a value of 0.0 for the minimization function. Starting with an initial guess that was close to the ‘perfect’ solution, the RBME system converged on the correct solution. However, perturbing the initial guess by as little as 1 part in 1000 caused the RBME system to converge on a different solution – where the temperature averaged (coincidentally) 7 °C instead of 25 °C. The fit between model and ‘measured’ TE/Ca profiles for this was excellent (Fig. 6), although not perfect. This behavior proved to be very robust: despite thousands of trials with both initial guesses based on perturbations of the perfect solution, and initial guesses a long way from the perfect solution, the system almost always converged on a solution with temperature averaging around 7 °C. This worrying result confirms that very good (but incorrect) local minima exist to the RBME equations, and indicates that if a clear global minimum does exist, the algorithms may not be able to find it. To test the effect that random (non T-derived) variation in the trace element profiles would have, I added Gaussian noise to the trace element profiles and again searched for solutions. The amplitude of the noise was calculated by spectral analysis (Sinclair et al., 2005, 2011) of a monthlyresolution tropical coral dataset (a Christmas Island coral analyzed for Ba/Ca by ICP-MS and Mg/Ca + Sr/Ca by ICP-OES). The noise produces misfits between the (smooth) ‘perfect solution’ and (noisy) data: calculating the minimization function for the perfect solution thus returns a value of 1.20, instead of 0.0. Optimization this solution improves the fit by matching some of noise variations, returning a minimization function value of 0.656 (Fig. 7). This theoretically represents the best possible fit to the noisy data. As before, perturbing the initial guess by too much results in a 7 °C solution with a minimization function

76

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

Table 4 Parameter values used to generate artificial data. Parametera

Value

P equation parameter C1: Temperature (K) at P = 0.5 P equation parameter C2: Rate at which sigmoid function changes (Sr/Ca)o (mol/mol) (Mg/Ca)o (mol/mol) (Ba/Ca)o (mol/mol) Temperature (at point x) (100 data points long)

289.295 1.60017 0.00954024 1.20667 7.77495E06 SST = 25.5 + 2.5sin(0.025px)sin(0.06px)

a

Temperature

Min Function = 0.021 Avg. T = 7.1 °C

Sr/Ca

9.90 9.95 10.00 10.05

-3

Mg/Ca (mol/mol x10 )

3.9 3.8 3.7 3.6 3.5 8.4

24 22

Min Function = 0.656 Avg T = 25.4 °C

20 7.8 7.2 6.6 9.80

Min Function = 0.656 Avg T = 23.3 °C

Sr/Ca

8.8 9.0 9.2 9.4

Min Function = 0.690 Avg T = 7.1 °C

'Observed' RBME Modeled

9.85 9.90 9.95 10.00 10.05

Model TE/Ca data are from the 23.3 °C solution

10.10

Mg/Ca

4.0 3.9 3.8 3.7 3.6 3.5 8.4

Ba/Ca

8.6

26

4.1

Mg/Ca

4.0

-6

-3

Mg/Ca (mol/mol x10 )

-6

'Observed' RBME Modeled

9.85

10.10 4.1

Ba/Ca (mol/mol x10 )

Temperature (°C)

24

-3

26

22 7.15 7.10 7.05 7.00 9.80

Temperature

28

Sr/Ca (mol/mol x10 )

28

Ba/Ca (mol/mol x10 )

-3

Sr/Ca (mol/mol x10 )

Temperature (°C)

Parameters come from a 25 °C solution to the Madrepora dataset.

Ba/Ca

8.6 8.8 9.0 9.2 9.4

0

20

40

60

80

Point in Time Series

Fig. 6. A solution to the Artificial Dataset. An artificial SST and TE/Ca dataset (red dashed lines) was generated from the data in Table 2. Solutions to the RBME equation were then generated using a hybrid of the Quasi-Newton and Nelder–Mead Simplex methods (blue solid lines). (A) The solution presented here converged on a temperature averaging 7.1 °C, providing an excellent, but not quite perfect, fit to the artificial TE profiles (B– D). Note the axis break and change in scale for Temperature axis (the modeled temperature is offset a long way from the ‘correct’ temperature and spans a significantly reduced seasonal scale). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of 0.690 (Fig. 7). By eye, this solution is indistinguishable from theoretical best solution. Several other solutions were also found with minimization function values between 0.690 and 0.700. These spanned temperatures averaging between 20 °C and 27 °C, and differing in amplitude by

0

20

40

60

80

Point in Time Series

Fig. 7. Three solutions to the Noisy Artificial Dataset. An artificial SST and TE/Ca dataset (red dashed lines) was generated from the data in Table 2, and the TE/Ca data were then perturbed with a Gaussian noise with an amplitude calculated from spectral analysis of a real coral dataset (see text). Solutions to the RBME equation were then generated using a hybrid of the Quasi-Newton and Nelder–Mead Simplex methods (blue solid lines). The 3 SST solutions (graphed in A) represent 3 of the best solutions found. The values of the minimization function for these solutions are added as annotations. Note the axis break and change in scale for Temperature axis (the amplitude of the 7.1 °C solution is magnified). The modeled TE/Ca datasets (blue solid lines in B–D) come from the 23.3 °C solution, but are essentially indistinguishable from the results for the 25.4 °C and 7.1 °C solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

up to 1.5 °C. The overall best solution found (Fig. 7) averaged 23 °C (instead of 25 °C) and had a minimization function value of 0.656, which is identical to the theoreti-

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

cal best solution. These results make it clear that that adding noise to the trace element profiles can degrade a perfect solution to the point where it may be no better than local minima, thus destroying the uniqueness of a single correct solution. 4.4. Summary The results of the RBME tests lead to the following conclusions: (1) There are a (potentially very large) number of solutions to the RBME equations (fits between modeled and observed Sr/Ca, Mg/Ca and Ba/Ca profiles) that are all very nearly as good as each other. (2) The average temperature of the initial guess can have a significant influence over the average temperature of the optimized solution. (3) The RBME method is able to find excellent solutions that nonetheless return ‘incorrect’ temperatures. (4) The true global minimum may be difficult for the search algorithm to find. (5) Noise can degrade the ‘true’ temperature solution to the point where it is no better (and can be worse than) other local minima. This leads to the important caveat: (6) The RBME method can return a solution that ‘confirms’ a users initial temperature guess, but unless the solution space is systematically tested, the user cannot be certain that this is the best or only temperature profile that could produce the observed trace element profiles.

5. DISCUSSION 5.1. Examining Assumption 1: constancy of (TE/Ca)0,i One of the critical assumptions (Section 2.2) that must be made for the RBME method to work is that the initial TE/Ca composition of the calcifying fluid (after modification by the coral, before precipitation commences) is the same for each cycle of replenishment/enrichment/precipitation. Gaetani et al. (2011) justify this assumption based on two pieces of evidence. Firstly, even assuming that CaATPase is 100% specific for Ca2+ (an unresolved point of debate in the literature – e.g., Ip and Krishnaveni, 1991; Ip and Lim, 1991; Ferrier-Page`s et al., 2002; Allison et al., 2011) Ca pumping will result in at most an 8% decrease in the TE/Ca ratio of the calcifying fluid relative to seawater – a result that broadly accords with modeling by Sinclair (2005a) and Sinclair and Risk (2006). Variation in the amount of pumping will therefore have only a small effect on the TE/Ca ratio of the calcifying fluid. Secondly, Gaetani et al. (2011) calculate that the change in proportion of Ca precipitated throughout the seasonal cycle can be entirely accounted for by the known (inorganic) temperature-dependence of the rate of CaCO3 precipitation. In

77

other words, the coral could be maintaining a seasonally constant rate of ECF replenishment and ATPase ion pumping (hence initial TE/Ca ratio) and the faster CaCO3 precipitation in summer would be enough to reproduce the calculated DP. The second point does not constitute proof that the coral is not modifying the ion pumping rate, and several lines of evidence suggest that corals do indeed exert control over the rate of ion pumping to their ECF under different conditions (e.g., Cohen et al., 2009; McCulloch et al., 2012; Gagnon et al., 2013; Venn et al., 2013). While the amplitude of this effect may be small, the RBME solutions presented here demonstrate that seasonal differences as small as 3% in the TE/Ca cycles characterize solutions that reproduce entirely incorrect temperature profiles (e.g., Fig. 4). Thus, even a small violation of Assumption 1 may potentially result in large uncertainties in the reconstructed temperature. 5.2. Examining Assumption 2: Ba/Ca and RB2E Assumption 1 carries a second implicit assumption: in order for Mg/Ca Ba/Ca and Sr/Ca to remain constant in the ECF, they must also be constant in ambient seawater from which they are derived. Magnesium is conservative in seawater, and the assumption of constancy is justified. Sr is a nearly-conservative element, although slight surface-water depletion (Alibert et al., 2003) by scavenging means that varying oceanic conditions (e.g., upwelling) may result in a seasonally-varying Sr/Ca that could bias the RBME calculation. Although the assumption of constant Sr may be justified, the sensitivity shown by the RBME system (Section 5.1) means that caution should be taken, especially in regions affected by upwelling. Unlike Sr and Mg, however, Ba is not close to being conservative in the ocean. Ba is enriched in riverine discharge and submarine groundwaters (Moore, 1997) and is strongly depleted in surface waters (Bacon and Edmond, 1972; Chan et al., 1976). It is also a highly bioactive element, being taken up and released by marine organisms and particles. There are therefore very few situations where Ba could be trusted to remain constant throughout a seasonal cycle. Moreover, corals have demonstrated anomalous skeletal variations that may be related to endogenous processes (Sinclair, 2005b). The combination of these factors means that Ba is a very poor choice for RBME temperature reconstruction. Tests of the RBME system with a perturbed Ba signal clearly demonstrate that non-constant Ba in the ECF can produce very anomalous T reconstruction. RBME temperature reconstruction, however, does not require 3 elements. If the equations for Ba are simply dropped from the system, then the modified scheme described in Section 3.3 becomes one where N subsamples gives 2N equations in N + 4 unknowns. This system becomes resolvable by the time N P 6. Although requiring a few more observations, there is no particular reason why RBME requires the use of Ba in addition to Sr and Mg. In Electronic Annex 2 I present alternative code for a two-element variation of the RBME system (dubbed the RB2E

78

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

system). Tests show that this scheme performs on par with the RBME system, returning similar temperature solutions to the 3 element calculation in cases where Ba/Ca in ambient water can be assumed constant (e.g., the Madrepora dataset and artificial datasets). 5.3. Other Assumptions The entire basis of RBME temperature reconstruction is built on the assumption (Assumption 3) that temperature is the only parameter affecting the precipitation efficiency of the coral such that there is a single value of FL (or P), and hence TE/Ca, for a given temperature. In essence, this encapsulates all of the necessary simplifications to the biomineralization physiology that enable it to be numerically soluble, and is common to all methods of temperature calibration (not just RBME). There are a number of external parameters in addition to temperature which could independently affect coral calcification physiology, including pH, nutrient/energy supply, light, disease, etc. A full review of these is beyond the scope of this paper, although there is clear evidence that deep sea corals growing in very constant temperature environments can display a wide range of Sr/Ca and Mg/ Ca values (e.g., Gagnon et al., 2007, 2012). Given the demonstrated sensitivity of the RBME method to subtle changes in seasonal TE/Ca variations, even second-order effects by these other environmental parameters could potentially affect the accuracy of temperature reconstruction. While Rayleigh Fractionation may remain the dominant process influencing trace-element variation, the RBME method should be applied with caution if there is evidence that other parameters might be influencing coral health. 6. CONCLUSIONS AND RECOMMENDATIONS RBME is a potentially very powerful system for paleotemperature reconstruction. Under well constrained conditions, this method offers the potential for eliminating some of the uncertainty caused by the possibility that individual corals exhibit individual-specific Sr/Ca paleothermometry calibrations. However, as I have demonstrated, the system suffers from the same problems that are common to most nonlinear high-dimensionality optimization techniques: namely sensitivity and bias-toward initial guesses, sensitivity to noise, degeneracy of solutions, computational overhead and the impossibility of determining if a global optimum has been identified. (Necessary) simplifications and arbitrary choices in the biomineralization equations also add another layer of uncertainty. RBME, therefore, is not a golden fix for all the problems of coral paleotemperature reconstruction, but should be regarded as one of a number of tools and techniques that could be employed. Because of these issues, output of RBME modeling needs to be presented with due diligence given to the uncertainties inherent in the method. To that end, I recommend the following strategies when undertaking RBME paleotemperature reconstruction:

(1) Downsample coral data to the lowest resolution that captures the full T cycle (e.g., a minimum of 6 samples/seasonal cycle). Reducing the number of data points reduces the dimensionality and therefore computational overhead associated with searching for optimums. (2) Use aggressive optimization routines using the betterbehaved equations presented here. Together with 1) this allows the user to run many trials (see points 3–7). (3) Use multiple different initial guesses for the temperature so that bias toward the initial guess is minimized and the widest possible range of solution space is searched. I would, for example, recommend bracketing the ‘expected’ temperature (by several °C beyond the maximum known constraints) including: a. Different average temperatures (e.g., 0.5 °C increments). b. Different amplitudes of T cycle (e.g., 1.0 °C increments). (4) Use different suites of the other fit parameters in the initial guess, especially if you are using the FL equation rather than the P equation (as this tends to produce more discrete local optimums). (5) Use different variations of the RBME system (e.g., Gaetani et al. original equations, Gaetani et al. original optimization criterion, RB2E etc). (6) Use different search heuristics or a hybrid strategy. For example, I found an effective combination was to do a wide search of solution space using the Quasi-Newton method (using multiple random perturbations from an initial guess), and to refine the best looking candidates using the Nelder–Mead Simplex (using a random perturbation each time it got cycle-locked). Other search heuristics might also be useable: initial tests with Simulated Annealing proved to be extremely slow, while Genetic Algorithms may have similar behavior to the Nelder– Mead Simplex, given their similarity (Weise, 2013). (7) Having run multiple trials, users should offer an estimate of the uncertainty of a preferred temperature solution by reporting (graphing) a range of solutions (I suggest all solutions with optimums within 20% of the best solution). (8) For any published solution, the user should report the full data associated with this solution (observed coral data, initial temperature guess, final fit parameters, details of the numerical scheme used, and – importantly – the value of the optimization function) so that others may test the result using different strategies. The RBME method is not perfect, but it does represent a significant and important new direction in coral paleotemperature reconstruction. I hope as our knowledge of biomineralization processes advances, and experiments yield better quantitative constraints on the impact of physiology on skeletal chemistry, that the RBME method will continue to evolve and improve.

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80 ACKNOWLEDGEMENTS I am very grateful for the support of Professor Alan Saul (Georgia Regents University, Department of Ophthalmology) who wrote the Igor Pro code for the Nelder–Mead Simplex algorithm, and who provided important feedback and discussion on the code and manuscript. I would also like to acknowledge helpful discussions with members of the Igor Pro email list, especially John Hellstrom (University of Melbourne), which helped a lot in the development of the code presented here. Finally, I wish to acknowledge the invaluable support of Professor Robert M. Sherrell who provided financial and logistical support during my employment at the Institute of Marine and Coastal Sciences, Rutgers University. This research was supported by NSF grants 0752544, 0962260 and 1103478. I am very grateful for insightful reviews by A. Gagnon and an anonymous reviewer that significantly improved this manuscript.

APPENDIX A. SUPPLEMENTARY DATA Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.gca.2015.01.006. REFERENCES Adkins J. F., Boyle E. A., Curry W. B. and Lutringer A. (2003) Stable isotopes in deep-sea corals and a new mechanism for “vital effects”. Geochim. Cosmochim. Acta 6, 1129–1143. Al-Horani F. A., Al-Moghrabi S. M. and de Beer D. (2003) The mechanism of calcification and its relation to photosynthesis and respiration in the scleractinian coral Galaxea fascicularis. Mar. Biol. 142, 419–426. Alibert C. and McCulloch M. T. (1997) Strontium/calcium ratios in modern Porites corals from the Great Barrier Reef as a proxy for sea surface temperature: calibration of the thermometer and monitoring of ENSO. Paleoceanography 12, 345–363. Alibert C., Kinsley L. P. J., Fallon S. J., McCulloch M. T., Berkelmans R. and McAllister F. (2003) Source of trace element variability in Great Barrier Reef corals affected by the Burdekin flood plumes. Geochim. Cosmochim. Acta 67, 231– 246. Allison N., Cohen I., Finch A. A. and Erez J. (2011) Controls on Sr/Ca and Mg/Ca in scleractinian corals: the effects of CaATPase and transcellular Ca channels on skeletal chemistry. Geochim. Cosmochim. Acta 75, 6350–6360. Bacon M. P. and Edmond J. M. (1972) Barium at GEOSECS III in the Southwest Pacific. Earth Planet. Sci. Lett. 16, 66–74. Beck J. W., Edwards R. L., Ito E., Taylor F. W., Recy J., Rougerie F., Joannot P. and Henin C. (1992) Sea-surface temperature from coral skeletal strontium/calcium ratios. Science 257, 644– 647. Chan L. H., Edmond J. M., Stallard R. F., Broecker W. S., Chung Y. C., Weiss R. F. and Ku T. L. (1976) Radium and barium at GEOSECS stations in the Atlantic and Pacific. Earth Planet. Sci. Lett. 32, 258–267. Cobb K. M., Charles C. D., Cheng H. and Edwards R. L. (2003) El Nino/Southern Oscillation and tropical Pacific climate during the last millennium. Nature 424, 271–276. Cohen A. L., Gaetani G. A., Lunda¨lv T., Corliss B. H. and George R. Y. (2006) Compositional variability in a cold-water scleractinian, Lophelia pertusa: new insights into “vital effects”. Geochem. Geophys. Geosyst. 7, Q12004.

79

Cohen A. L., McCorkle D. C., de Putron S., Gaetani G. A. and Rose K. A. (2009) Morphological and compositional changes in the skeletons of new coral recruits reared in acidified seawater: insights into the biomineralization response to ocean acidification. Geochem. Geophys. Geosyst. 10, Q07005. Correge T., Gagan M. K., Beck J. W., Burr G. S., Cabioch G. and Le Cornec F. (2004) Interdecadal variation in the extent of South Pacific tropical waters during the Younger Dryas event. Nature 428, 927–929. Crowley T. J., Quinn T. M. and Hyde W. T. (1999) Validation of coral temperature calibrations. Paleoceanography 14, 605–615. de Villiers S., Shen G. T. and Nelson B. K. (1994) The Sr/Catemperature relationship in coralline aragonite: influence of variability in (Sr/Ca)seawater and skeletal growth parameters. Geochim. Cosmochim. Acta 58, 197–208. Dennis, Jr., J. E. and Schnabel R. B. (1996) Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice Hall, Englewood Cliffs. Dietzel M., Gussone N. and Eisenhauer A. (2004) Co-precipitation of Sr2+ and Ba2+ with aragonite by membrane diffusion of CO2 between 10 and 50 °C. Chem. Geol. 203, 139–151. Elderfield H., Bertram C. J. and Erez J. (1996) A biomineralization model for the incorporation of trace elements into foraminiferal calcium carbonate. Earth Planet. Sci. Lett. 142, 409–423. Fallon S. J., McCulloch M. T. and Alibert C. (2003) Examining water temperature proxies in Porites corals from the Great Barrier Reef: a cross shelf comparison. Coral Reefs 22, 389–404. Felis T., Lohmann G., Kuhnert H., Lorenz S. J., Scholz D., Patzold J., Al-Rousan S. A. and Al-Moghrabi S. M. (2004) Increased seasonality in Middle East temperatures during the last interglacial period. Nature 429, 164–168. Felis T., Suzuki A., Kuhnert H., Dima M., Lohmann G. and Kawahata H. (2009) Subtropical coral reveals abrupt earlytwentieth-century freshening in the western North Pacific Ocean. Geology 37, 527–530. Ferrier-Page`s C., Boisson F., Allemand D. and Tambutte´ E. (2002) Kinetics of strontium uptake in the scleractinian coral Stylophora pistillata. Mar. Ecol. Prog. Ser. 245, 93–100. Gaetani G. A. and Cohen A. L. (2006) Element partitioning during precipitation of aragonite from seawater: a framework for understanding paleoproxies. Geochim. Cosmochim. Acta 70, 4617–4634. Gaetani G. A., Cohen A. L., Wang Z. and Crusius J. (2011) Rayleigh-based, multi-element coral thermometry: a biomineralization approach to developing climate proxies. Geochim. Cosmochim. Acta 75, 1920–1932. Gagan M. K., Ayliffe L. K., Hopley D., Cali J. A., Mortimer G. E., Chappell J., McCulloch M. T. and Head M. J. (1998) Temperature and surface-ocean water balance of the midHolocene tropical western pacific. Science 279, 1014–1018. Gagnon A. C., Adkins J. F., Fernandez D. P. and Robinson L. F. (2007) Sr/Ca and Mg/Ca vital effects correlated with skeletal architecture in a scleractinian deep-sea coral and the role of Rayleigh fractionation. Earth Planet. Sci. Lett. 261, 280–295. Gagnon A. C., Adkins J. F. and Erez J. (2012) Seawater transport during coral biomineralization. Earth Planet. Sci. Lett. 329– 330, 150–161. Gagnon A. C., Adkins J. F., Erez J., Eiler J. M. and Guan Y. (2013) Sr/Ca sensitivity to aragonite saturation state in cultured subsamples from a single colony of coral: mechanism of biomineralization during ocean acidification. Geochim. Cosmochim. Acta 105, 240–254. Ip Y. K. and Krishnaveni P. (1991) The incorporation of strontium (90Sr2+) into the skeleton of the hermatypic coral Galaxea fascicularis. J. Exp. Zool. 258, 273–276.

80

D.J. Sinclair / Geochimica et Cosmochimica Acta 154 (2015) 66–80

Ip Y. K. and Lim A. L. L. (1991) Are calcium and strontium transported by the same mechanism in the hermatypic coral Galaxea fascicularis?. J. Exp. Biol. 159, 507–513. Kinsman D. J. J. and Holland H. D. (1969) The co-precipitation of cations with CaCO3—IV. The co-precipitation of Sr2+ with aragonite between 16° and 96 °C. Geochim. Cosmochim. Acta 33, 1–17. Lagarias J. C., Reeds J. A., Wright M. H. and Wright P. E. (1998) Convergence properties of the Nelder–Mead Simplex Method in low dimensions. SIAM J. Optim. 9, 112–147. Lea D. W. (2003) Elemental and isotopic proxies of marine temperatures. In The Oceans and Marine Geochemistry (ed. H. Elderfield). Elsevier-Pergamon, Oxford, pp. 365–390. Linsley B. K., Wellington G. M. and Schrag D. P. (2000) Temperature variability in the subtropical South Pacific from 1726 to 1997 A.D. Science 290, 1145–1148. Linsley B. K., Wellington G. M., Schrag D. P., Ren L., Salinger M. J. and Tudhope A. W. (2004) Geochemical evidence from corals for changes in the amplitude and spatial pattern of South Pacific interdecadal climate variability over the last 300 years. Clim. Dyn. 22, 1–11. Marriott C. S., Henderson G. M., Belshaw N. S. and Tudhope A. W. (2004) Temperature dependence of d7Li, d44Ca and Li/Ca during growth of calcium carbonate. Earth Planet. Sci. Lett. 222, 615–624. Marshall A. T. and Clode P. (2004) Calcification rate and the effect of temperature in a zooxanthellate and an azooxanthellate scleractinian reef coral. Coral Reefs 23, 218–224. Marshall J. F. and McCulloch M. T. (2002) An assessment of the Sr/Ca ratio in shallow water hermatypic corals as a proxy for sea surface temperature. Geochim. Cosmochim. Acta 66, 3263– 3280. McConnaughey T. (1986) Oxygen and Carbon Isotope Disequilibria in Galapagos Corals: Isotopic Thermometry and Calcification Physiology. University of Washington, Seattle, WA, United States, p. 340. McConnaughey T. (1989a) 13C and 18O isotopic disequilibrium in biological carbonates: I. Patterns. Geochim. Cosmochim. Acta 53, 151–162. McConnaughey T. (1989b) 13C and 18O isotopic disequilibrium in biological carbonates: II. In vitro simulation of kinetic isotope effects. Geochim. Cosmochim. Acta 53, 163–171. McCulloch M., Trotter J., Montagna P., Falter J., Dunbar R., Freiwald A., Fo¨rsterra G., Lo´pez Correa M., Maier C., Ru¨ggeberg A. and Taviani M. (2012) Resilience of cold-water scleractinian corals to ocean acidification: boron isotopic systematics of pH and saturation state up-regulation. Geochim. Cosmochim. Acta 87, 21–34. McKinnon K. I. M. (1999) Convergence of the Nelder–Mead Simplex Method to a nonstationary point. SIAM J. Optim. 9, 148–158. Min G. R., Edwards R. L., Taylor F. W., Recy J., Gallup C. D. and Beck J. W. (1995) Annual cycles of U/Ca in coral skeletons and U/Ca thermometry. Geochim. Cosmochim. Acta 59, 2025– 2042. Mitsuguchi T., Matsumoto E., Abe O., Uchida T. and Isdale P. J. (1996) Mg/Ca thermometry in coral skeletons. Science 274, 961–963. Moore W. S. (1997) High fluxes of radium and barium from the mouth of the Ganges–Brahmaputra River during low river

discharge suggest a large groundwater source. Earth Planet. Sci. Lett. 150, 141–150. Nurhati I. S., Cobb K. M., Charles C. D. and Dunbar R. B. (2009) Late 20th century warming and freshening in the central tropical Pacific. Geophys. Res. Lett. 36, L21606. Quinn T. M. and Sampson D. E. (2002) A multiproxy approach to reconstructing sea surface conditions using coral skeleton geochemistry. Paleoceanography 17, 1062. Quinn T. M., Taylor F. W. and Crowley T. J. (2006) Coral-based climate variability in the Western Pacific Warm Pool since 1867. J. Geophys. Res.: Oceans 111, C11006. Sinclair D. J. (1999) High Spatial-Resolution Analysis of Trace Elements in Corals Using Laser Ablation ICP-MS, Research School of Earth Sciences. Australian National University, Canberra, p. 386. Sinclair D. J. (2005a) Correlated trace element ‘vital effects’ in tropical corals: a new tool for probing biomineralization chemistry. Geochim. Cosmochim. Acta 69, 3265–3284. Sinclair D. J. (2005b) Non river-flood barium signals in the skeletons of corals from coastal Queensland, Australia. Earth Planet. Sci. Lett. 237, 354–369. Sinclair D. J. and Risk M. J. (2006) A numerical model of traceelement coprecipitation in a physicochemical calcification system: application to coral biomineralization and traceelement ‘vital effects’. Geochim. Cosmochim. Acta 70, 3855– 3868. Sinclair D. J., Kinsley L. P. J. and McCulloch M. T. (1998) High resolution analysis of trace elements in corals by laser-ablation ICP-MS. Geochim. Cosmochim. Acta 62, 1889–1901. Sinclair D. J., Sherwood O. A., Risk M. J., Hillaire-Marcel C., Tubrett M., Sylvester P., McCulloch M. T. and Kinsley L. P. J. (2005) Testing the reproducibility of Mg/Ca profiles in the deep-water coral Primnoa resedaeformis: putting the proxy through its paces. In Cold Water Corals and Ecosystems – Selected Papers from the 2nd International Symposium on Deep Sea Corals (eds. A. Freiwald and J. M. Roberts). SpringerVerlag, Berlin Heidelberg, pp. 1039–1060. Sinclair D. J., Allard G., Williams B., Ross S. and Risk M. J. (2011) Reproducibility of trace element profiles in a specimen of the deep-water bamboo coral Keratoisis sp. Geochim. Cosmochim. Acta 75, 5101–5121. Smith S. V., Buddemeier R. W., Redalje R. C. and Houck J. E. (1979) Strontium–calcium thermometry in coral skeletons. Science 204, 404–406. Stephans C. L., Quinn T. M., Taylor F. W. and Corre`ge T. (2004) Assessing the reproducibility of coral-based climate records. Geophys. Res. Lett. 31, L18210. Venn A. A., Tambutte´ E., Holcomb M., Laurent J., Allemand D. and Tambutte´ S. (2013) Impact of seawater acidification on pH at the tissue–skeleton interface and calcification in reef corals. Proc. Natl. Acad. Sci. USA 110, 1634–1639. Weise T. (2013) Global Optimization Algorithms Theory and Application. University of Kassel, p. 820. Zinke J., Dullo W. C., Heiss G. A. and Eisenhauer A. (2004) ENSO and Indian Ocean subtropical dipole variability is recorded in a coral record off southwest Madagascar for the period 1659 to 1995. Earth Planet. Sci. Lett. 228, 177–194. Associate editor: Yair Rosenthal