© Author(s) 2010. CC Attribution 3.0 License.

Abstract. A long standing problem in paleoceanography concerns the reconstruction of water temperature from δ18O carbonate. It is problematic in the case of freshwater influenced environments because the δ18O isotopic composition of the ambient water (related to salinity) needs to be known. In this paper we argue for the use of a nonlinear multi-proxy method called Weight Determination by Manifold Regularization (WDMR) to develop a temperature reconstruction model that is less sensitive to salinity variations. The motivation for using this type of model is twofold: firstly, observed nonlinear relations between specific proxies and water temperature motivate the use of nonlinear models. Secondly, the use of multi-proxy models enables salinity related variations of a given temperature proxy to be explained by salinity-related information carried by a separate proxy. Our findings confirm that Mg/Ca is a powerful paleothermometer and highlight that reconstruction performance based on this proxy is improved significantly by combining its information with the information for other trace elements in multi-proxy models. Although the models presented here are black-box models that do not use any prior knowledge about the proxies, the comparison of model reconstruction performances based on different proxy combinations do yield useful information about proxy characteristics. Using Mg/Ca, Sr/Ca, Ba/Ca and Pb/Ca the WDMR model enables a temperature reconstruction with a root mean squared error of ± 2.19 °C for a salinity range between 15 and 32.


Introduction
To improve our understanding of global change and assess human impact on global warming, reconstructions of past temperatures are essential.Such reconstructions are mostly based on the analysis of trace elements and isotopes in accreting biogenic or abiogenic substrates, called archives.The choice of the parameters (called proxies) to be analysed is based on prior knowledge of their relationship with an environmental variable as derived by observing such relationship in the present-day situation (Kucera et al., 2005).Several natural archives in the terrestrial and the marine environment record environmental information in their trace element and isotope profiles.Bivalve shells, in particular, represent a suitable archive for reconstructing seasonal and long term variations of ambient water conditions and many elemental and isotopic temperature proxies have been proposed and discussed for these archives (e.g., Epstein et al., 1953a, b;Klein et al., 1996a, b;Wanamaker et al., 2006;Freitas et al., 2009).Indeed, bivalves are sensitive to environmental conditions, have a global distribution, and are commonly found in archaeological sites (Pearce and Mann, 2006;Klunder et al., 2008;Butler et al., 2009).
Bivalve shells offer thus the potential for reconstructing environmental conditions for a wide variety of aquatic environments, including fresh water systems (Versteegh, 2009), estuarine and marine environments from tropical (Aubert et al., 2009) to cold polar regions (Tada et al., 2006).
These and other studies reveal that though a given proxy may correlate well with an environmental parameter, the data usually show significant variation around the regression line, Published by Copernicus Publications on behalf of the European Geosciences Union.
reflecting that the process of proxy-incorporation is much more complex than assumed originally (Wanamaker et al., 2007;Gillikin et al., 2005).
Most water temperature reconstructions based on biogenic carbonates are based on δ 18 O records.For instance, for the common blue mussel (Mytilus edulis; the species studied in this paper) it has been shown that temperature reconstructions from shell δ 18 O records can achieve an excellent accuracy of 0.57 • C in Root Mean Squared Error (RMSE) (Wanamaker et al., 2007).However, this paleothermometer equation requires that the δ 18 O value of the ambient water be known.This is obviously not possible for archeological specimens and given that the δ 18 O value of the ambient water strongly depends on salinity (a salinity variation of one can incorrectly be interpreted as a change of 1 • C in water temperature), a proxy or model which is less sensitive to salinity variations may therefore significantly improve paleotemperature reconstructions (Faure, 1986).
Several alternative (salinity-robust) temperature proxies have been proposed (e.g., Mg/Ca-ratios by Klein et al., 1996b; Sr/Ca ratios by Foster et al., 2009).However, proxies mostly appear influenced by several environmental parameters (e.g., Elliot et al., 2009;Foster et al., 2009).Moreover, the fact that these potential temperature proxies are recorded in biogenic material, makes them subject to physiology-related biases such as kinetic effects (Lorrain et al., 2005), methabolic effects (Strasser et al., 2008) and ontogenetic effects (Elliot et al., 2009).It becomes more and more clear that biomineralisation is a complex process, whose adequate study ideally requires the involvement of several disciplines (Weiner and Dove, 2003).
In the present paper we investigate whether more complex, non-linear models are better suited for describing the integrated impact of environmental conditions, physiological state of the organism and a complex suit of biochemical and chemical processes, on proxy incorporation during bivalve shell growth.We propose to combine a suit of proxies, differentially influenced by environmental and biological controls, into a multi-proxy model.Multi-proxy models offer the advantage that variation in the different proxies yields information that is useful to resolve environmental and biological interferences.The proposed multi-proxy model combines information on elemental ratios (in this case Mg/Ca, Sr/Ca, Ba/Ca and Pb/Ca) based on the two general (statistical) assumptions: (i) the proxies are influenced by the same environmental and intrinsic parameters, and therefore combining them may help explaining variation that was not understood before; (ii) the proxies are likely influenced to different degrees by temperature variation and, therefore, using temperature information derived from each of the proxies will yield more robust temperature reconstructions.
The models presented in this paper are not based on a mechanistic understanding of the incorporation mechanisms of the proxies.However, along this paper it becomes clear that the studied proxies do not contribute equally to the final temperature reconstructions.The contribution of each proxy was calculated, from the temperature reconstruction performances of different proxy combinations.

Why multi-proxy models?
As mentioned in the introduction two reasons can be invoked for promoting the use of multi-proxy models.
The first and most important reason (i, above) is synthesized by the set of Eqs.(1), representing a linear multiple regression model with a limited number of parameters.
These equations express how a number of environmental parameters (e.g., temperature, salinity, chlorophyll concentration) all contribute to the final trace element signature of the archive.Solving this set of equations for the environmental parameters involves a new set of equations in which all environmental parameters can be described by multi-proxy equations, implying that all proxies add some information to the final paleo-temperature equation.For example: by combining an element that is mainly influenced by salinity with another element influenced by both, temperature and salinity, it is possible to construct a model that is more robust across a range of salinities.
Though Eqs. (1) as shown include only environmental parameters (Temp., Sal., Chlorophyll) it is clear that other, organism-related parameters may be included as well.For example shell growth, spawning events and metabolic activity can be included.Such multi-proxy equations would resolve part of the "vital effect" commonly invoked to explain a chemical response that is not understood.
Note that solving a non-linear model with a large number of parameters is much more complex, but the idea behind it would be the same.Although it is algebraically possible to reverse such multiple regression equations when there are as many proxies as environmental parameters, this would induce large errors on the estimated parameters.Therefore the multi-proxy models obtained in this paper are considered as black-box models that cannot be reversed to obtain a mechanistic understanding in the proxy incorporation.
(ii) A second drive for using multi-proxy models is rather intuitive.Assuming that different proxies each carry some temperature information, it seems reasonable that a model based on the information of several proxies will yield more robust and accurate reconstructions, though this requires proper weighing of each proxy.The weight given to a proxy depends on the quality of the proxy environmental relationship in the calibration or training set and less importance is given to proxies that show a less clear or noisy relationship with the environmental condition (temperature).Noise may result from the large influence of an additional environmental or biological condition or from measurement uncertainty.This means that proxies that have a large load of environmental information have the largest influence on the final reconstruction, even though other proxies are used to explain or confirm parts of the signal.
Despite these clear advantages, applications of multiproxy models are scarce in bivalve sclerochronology literature.Some steps in this direction are made by Klein et al. (1996b) and Schöne et al. (2006), though these authors rather use a secondary proxy to confirm a signal that is revealed by a primary proxy.Gentry et al. (2009) and Bice et al. (2006) discuss two approaches in which the influence of salinity on δ 18 O carbonate is eliminated by formulating an initial guess of the δ 18 O water using information from a secondary proxy.However, to the best of our knowledge multi-proxy models in which a given environmental parameter is described by a combination of several proxies have not been published yet, one exception being the work of Freitas et al. (2006) who demonstrate that a linear multiple regression analysis using Sr and Mg, significantly improves temperature estimates.

Why nonlinear multi-proxy models?
Considering that physiological processes are nonlinearly influenced by environmental conditions, as is the case for instance for temperature, plankton blooms (Cloern et al., 1995), optimal feeding temperature (Yukihira et al., 2000), the occurrence of nonlinear relationships between proxies and environmental conditions does not come as a surprise.Figure 2 shows an example of a substantially nonlinear relationship between bivalve shell proxies and water temperature (Vander Putten et al., 2000), highlighting a direct but complex influence of temperature on trace element uptake.However, such relationships have been traditionally described using linear equations (Klein et al., 1996b;Wanamaker et al., 2008), though some recent publications describe or advocate the use of inverse exponentials (Clarke et al., 2009), exponentials (Freitas et al., 2005) and even dynamical (Klunder et al., 2008) relationships.
Nonlinear relationships between proxies and environmental conditions are difficult to describe in a single mathematical equation but they can be modeled by several modern multivariate statistical techniques (Izenman, 2008).Most scientists are familiar with the classical linear multiple regression and dimensionality reduction methods, such as Principle Component Analysis (PCA), Cluster Analysis, etc.However, these methods are developed to detect linear relationships and are not applicable to datasets that behave substantially nonlinear.To detect nonlinear relationships in a multi-dimensional space, recently developed multivariate statistical tools are needed (Izenman, 2008).The best known nonlinear multivariate statistical techniques in paleoclimatology are Artificial Neural Networks which are being used for reconstructing ENSO events from coral records (Juillet-Leclerc et al., 2006) and in dendrochronology to reconstruct precipitation rates (Woodhouse, 1999) and temperature (Guiot et al., 2005).However, other techniques such as Support Vector Machines and Manifold Learning can be used for the same purpose (Bauwens et al., 2010).
Different nonlinear multivariate statistical techniques are thus available to analyze multidimensional datasets, but the choice of a specific technique will depend on characteristics of the dataset such as number of data, intrinsic variance, smoothness, periodicity.As a consequence, each dataset has its own "best method".In a previous paper we compared three nonlinear multiple regression methods: two of the three nonlinear regression methods explored in that paper reduce the multi-proxy problem into a single dimensional problem by observing that the proxies lie on a onedimensional manifold.One of the two is based on intuition and tailored for temperature reconstruction using bivalve shells.The other is a new system identification approach, Weight Determination by Manifold Learning (WDMR), and based on manifold learning.The third approach,Support Vector Regression (SVR), does not rely on an assumption of a manifold in the proxy space; it rather increases the dimensionality of the problem by creating "new proxies" from nonlinear combinations of the original proxy data.In Bauwens et al. (2010) it is concluded that manifold based methods are the most powerful tools for reconstructing paleo-environmental conditions based on proxy records in shells of short-lived bivalves, suggesting that the proxyenvironmental relationships are straightforward and no extra information is gained by using a more complex SVR model (Bauwens et al., 2010).
In the present paper we use the manifold based method called Weight Determination Manifold Regularization (WDMR) (Ohlsson et al., 2008(Ohlsson et al., , 2009a, b) , b) to build a salinity-robust model for reconstructing temperature using shells of the common blue mussel Mytilus edulis.

Raw data
The trace element datasets used in this paper were originally published by Vander Putten et al. (2000) and Gillikin et al. (2006a, b).Both datasets consist of spatially well resolved measurements of Mg/Ca, Sr/Ca, Ba/Ca and Pb/Ca ratios along the shell's main growth axis for approximately two years old M. edulis specimens.For both studies laser ablation craters (from LA-ICP-MS analyses) were produced in the calcitic layer of the shell.The ablation craters were approximately 50 µm in diameter and were spaced every 250 µm.For each shell of 45 to 65 ablations were performed over the shell section that grew during the period of monitoring.All specimens were sampled in the Scheldt Estuary (The Netherlands, Belgium); the exact geographical position of the four study sites is shown in Fig. 1.The reader interested in more details about these data sets is referred to the papers by Vander Putten et al. (2000) and Gillikin et al. (2006a, b).
The Gillikin et al. dataset consists of proxy profiles for a single shell sampled at the Knokke site and monitored from February to September 2002.Since the blue mussel stops growing when temperature drops below 8 • C (usually in autumn; Gillikin et al., 2009), the analyzed February to September period closely corresponds to a complete growth season.The Vander Putten et al. data set concerns seven blue mussel shells from Terneuzen, four shells from Ossenisse and four shells from Breskens (Figs. 1 and 3).These data cover the period from April to June 1996, and do not cover the full growth season of the mussel, though it includes the spring period when shell accretion is fastest and variations in trace element concentrations largest.The total dataset covers a salinity range from 15 to 31 and a temperature range from 6.8 • C to 18.6 • C for 1996 and from 8.7 • C to 19.3 • C for 2002.

Linking proxy data to environmental information
The proxies were measured along the largest growth axis (i.e., along a distance scale) starting at the margin of the shell moving towards the umbo.Since temperature measurements are obtained on a time scale, linking proxy data to environmental information is not straightforward.For both data sets the link between spatial and temporal scales was established using the anchor point-method (Paillard et al., 1996), implying that between anchor points, growth is assumed linear.The anchor points for the Vander Putten et al. shells were T 0 (marking on the shell), T final (date of collection) and recognizable patterns in trace-elemental chemistry, such as a conspicuous Ba-peak associated with the spring bloom.The anchor points for the Gillikin dataset were obtained from pattern similarities between the δ 18 O profile of shell carbonate and the water temperature profile monitored at the study site.The assumption of subsequent linear growth events, however, is an approximation since shell growth is variable (Schöne et al., 2005).Other methods to reconstruct the shell growth, as reviewed in de Brauwere et al. ( 2008), could not be applied to the datasets used in the present study, since these methods are designed for periodic signals and are not applicable to records covering only a single season, as is the case here.

Normalized data
Proxy signals in different specimens from the same species sampled at the same location are often similar but seldom identical.Since environmental variability is unlikely over the small spatial scale of a mussel bank, the variation can be seen as an intrinsic and unexplained variation that we shall call "noise".Besides noise, site-and year-specific variation can occur.By normalizing the data the reconstructed environmental parameter will become dependent on the overall shape of the proxy record.Normalization was done by dividing the data by the standard deviation and subtracting the mean.This offers the advantage of the data becoming less sensitive to site and year specific variability as well as concentration shifts (see Fig. 3; and also Stecher et al., 1996;Gillikin et al., 2008) since these effects will be filtered out.The disadvantages, however, are that (1) some potential useful information may be lost and (2) that temperature reconstructions are not possible from individual measurements since the normalization of data is only possible for more than one data point.Note also that the normalization of data is preferentially done on data of a whole season to avoid over or under estimations of the reconstructed temperatures.

Training and validation data
The data were divided into two parts: a training dataset consisting of 6 shells from the Terneuzen site in the Scheldt Estuary and a validation dataset consisting of shells from all 4 sites along the Scheldt Estuary, i.e., one shell from Knokke sampled in 2002, four shells from Breskens, one shell from Terneuzen and four shells from Ossenisse, all sampled in 1996 (Figs. 1 and 3).The fact that the Knokke specimen  3 The methods

Linear multiple regression
Linear multiple regression is the most commonly used multivariate method to describe the linear relationship between two or more explanatory variables (here proxies) and a response variable (here temperature).This is done by fitting a linear equation to observed data.An equation similar to the equations in Eq. ( 1) (Temp = α 1 .Mg + α 2 .Sr + α 3 .Ba + ... + C α ) describes how temperature co-varies with the proxies.A limited number of parameters α 1 , α 2 ,...α n define the slope of the regression line and a coefficient C α defines the offset.
The main advantages of linear multiple regressions are that appropriate toolboxes are available on all statistical software packages and that models have a limited number of parameters and model outputs which renders interpretation easier.A large disadvantage, however, is that linear models are not able to fit nonlinear relationships which are likely to occur in biogenic archives.

Weight Determination by Manifold Regularization (WDMR)
The mathematical details of the method called Weight Determination by Manifold Regularization (WDMR) are beyond the scope of the present paper and the interested reader is referred to (Ohlsson et al., 2008(Ohlsson et al., , 2009a, b), b).
Interested users can also download a Matlab WDMR toolbox that is added as supplementary material to this paper, although we recommend contacting the corresponding author to ensure correct use of the WDMR toolbox.
In the following we briefly describe the concept of this approach.Manifold learning is an umbrella term for methods for describing low-dimensional structures in data.A manifold is defined as a low-dimensional structure which underlies a collection of high dimensional data, for example a curve in the space of Mg, Sr, Ba and Pb concentrations.An algorithm that builds on concepts from manifold learning is the nonlinear semi-supervised regression method called Weight Determination by Manifold Regularization (WDMR) (Ohlsson et al., 2008).WDMR, like a manifold learning algorithm, finds descriptions of manifolds but unlike most manifold learning methods WDMR can utilize a training set for the description.If the temperature associated with a specific measurement of Mg, Sr, Ba and Pb in the training set is known, that information can be used in WDMR to impose a one-dimensional description of the curve imitating the temperature.In the case that proxy composition are controlled solely by water temperature, concentrations of Mg, Sr, Ba and Pb would be restricted to a one-dimensional curve in the four-dimensional measurement space with each position on the curve having a temperature value associated with it.As a result the curve can be parameterized by the water temperature.The computed WDMR-model can then be used to estimate the water temperature for any other dataset of Mg, Sr, Ba and Pb.As in all real-world problems there is noise associated with all measurements.And more importantly, it is unlikely that the concentrations of Mg, Sr, Ba and Pb will only depend on water temperature.Rather, they will depend also on other conditions such as salinity, food availability, shell growth or metabolism) and therefore the data will scatter around a one-dimensional curve in the Mg, Sr, Ba and Pb space.
The assumption of a one-dimensional manifold is therefore only an approximation, but the performance of the computed models shows that this approximation is appropriate (or at least a better approximation than assuming linearity).

Method
To investigate the benefit of using nonlinear methods rather than linear methods we compared the reconstruction performance of models generated using WDMR with models obtained by classical linear multiple regression.Six shells from Terneuzen were used to train both the linear model and the WDMR model.The linear multiple regression analysis where don on not normalized data, since these analysis are traditionaly done on raw data.The model performances were calculated for the four validation sets consisting of shells from the 4 study sites, including one additional shell from the training site (see Fig. 1).
To calculate the model performance the Root Mean Squared Error (RMSE) between measured and reconstructed temperatures for each data point was used.The reconstructed temperatures for the nonlinear WDMR model and the linear multiple regression model were compared and the differences between their RMSE used to verify whether some proxy combinations benefited more than others from the nonlinear model.

Results
The RMSE are smaller for the WDMR than for the linear approaches .The nonlinear WDMR model results in a better reconstruction of the seasonal temperature pattern for the Knokke site, as shown in Fig. 4. Also for the three other sites and for most proxy combinations the nonlinear WDMR model performs better than the linear multiple regression models in RMSE-sense (Fig. 5).This is true, in particular for the temperature reconstructions at Terneuzen and Knokke where only the Sr-only and the combined SrPb proxies do better with a linear model.The reconstruction performance of the nonlinear WDMR model is up to 1.5 • C better than the one for the linear model.Furthermore, the performance of nonlinear models is increasing when more proxies are included.This result confirms that relationships between a proxy and the controlling environmental condition can indeed be nonlinear.However, the weaker reconstruction performances of the nonlinear model for the Breskens and Ossenisse sites indicate that the nonlinear model over-fits the training data for some proxy combinations, such that in these cases linear models result in better reconstruction (Fig. 6).This is in particular true for Ba at the Breskens site, showing a distinct site-specific behavior which results in the linear model performing better.

Discussion
For most proxy combinations Fig. 4 clearly shows that the nonlinear WDMR model results in more accurate temperature reconstructions than the linear multiple regressions.The reconstruction performance of the nonlinear model is up to 1.5 • C better than for the linear model.However, Fig. 5 also shows that some proxy combinations do not benefit from the nonlinear model.A linear model is less sensitive to model errors related to over-fitting.The temperature reconstructions from the Breskens shell, for instance, are improved when using a linear model based on proxy combinations containing Ba information.The site specificity of Ba that can be observed in Figs. 3 and 5 is discussed later in Sect.5.3.2.Several relationships between proxy and environmental control factors reported in literature, behave linearly (Wanamaker et al., 2008;Carre et al., 2006) and in these cases a linear model with a lower number of parameters is still preferable.However, this should not be a reason for not using nonlinear methods since nonlinear methods since a linear model is, generally spoken, a special case of a nonlinear one and therefore nonlinear methods can fit linear data.

Method
To investigate the benefit of a multi-proxy approach using the WDMR method and to examine the contribution of the different proxies, different models were constructed based on a limited number of proxies.In total 15 combinations of proxies were investigated.The RMSE values obtained on the validation data were used to quantify the model performances.For the nonlinear models seven unique contribution factors were defined in order to quantify the contribution of each proxy.Every contribution factor quantifies how much a specific proxy contributes to a specific model; in other words the contribution factor informs on how much the RMSE decreases by including the information of the investigated proxy into a specific model.For example one of the seven contribution factors for Mg is 1.62.This means that the RMSE of a MgSr model was 1.62 lower than the RMSE of a Sr-only model.Negative contribution factors, on the other hand, reflect that including a specific proxy in the model has a negative influence the on model reconstruction performance.All contribution factors are defined as the difference between the RMSE between two model configurations (i.e.models run with different combinations of elemental ratios) (Table 1).This enables evaluating model performance change due to inclusion of additional proxies.All trace element combinations were tested for their robustness to salinity by using the different models to reconstruct the temperature based on the validation shells from the four sites along the estuarine salinity gradient.

Results
Figure 6 demonstrates that the four-proxy model generated with the WDMR method is relatively insensitive to changes in salinity, since the model is able to reconstruct the temperature for all study sites along the estuarine salinity gradient, without systematic errors due to differences in salinity.The overall trend of reconstructed temperatures is very similar to the measured temperature, but the reconstructed temperature profiles show more variability.
Though the best reconstruction is obtained for the validation shell from the same study site and collected at the same time as the training shells (RMSE = ±1.29 • C), the temperature The reconstruction performance of models trained for different proxy combinations is shown in Fig. 7.In general RMSE decreases with an increasing number of proxies.This trend is also observed in Table 1 where it is demonstrated that the use of an additional proxy in a multi-proxy model greatly improves reconstruction performance since most contribution factors (i.e.RMSE with proxy -RMSE without proxy) are positive.The benefit of using a multi-proxy model is thus significant, although it is clear that not all proxies contribute equally to the final reconstruction and the fourproxy model is not necessarily the best model.
Table 1 shows that on average all proxies contribute positively to the final reconstruction.Mg can improve the RMSE of a temperature reconstruction with 0.72, on average.Ba improves the RMSE of a temperature reconstruction with 0.42.Pb and Sr, however, show lower contribution factors of 0.20 and 0.04, respectively.The average contribution factors shown in Table 1 thus suggest that Mg and Ba contribute the most to an accurate temperature reconstruction.Ba, however, shows several negative contribution factors for the Breskens site, revealing site specific effects.However, information stored in the Sr-signature of the shell almost completely compensates for these site specific effects.This can clearly be seen by comparing the performance of the MgBa-model with the one of the MgSrBa-model in Fig. 7, with the latter yielding fairly accurate and salinity robust temperature reconstructions.Adding Pb to this MgSrBamodel does slightly improve the reconstruction, although by not more than 0.2 • C.

Discussion
Using the WDMR-method to construct paleo-thermometer models yields accurate temperature reconstructions for shells from Terneuzen where the training set was sampled.This reconstruction shows that it is possible to reconstruct the temperature based on Mg, Sr, Ba and Pb.The reconstruction performance is slightly poorer for shells from the other sites suggesting that the model is sensitive to site-specific variations.However, considering the salinity range from 32 (Knokke) to 15 (Ossenisse), the reconstruction performance (RMSE lower than 2.19 • C) for shells from a different site (and salinity) than the training set, is promising.Compared to other approaches for reconstructing water temperature based on the blue mussel archive (Epstein et al., 1953b;Wanamaker et al., 2006;Klein et al., 1996b) the performance of the method proposed here is of similar standard, if not better.
The multi-proxy model presented in this paper is built on four proxies of which two (Ba and Pb) were previously not considered to have potential as paleo-thermometers.It is thus probable that this method will provide even better reconstructions when trained on a set of well known temperature sensitive proxies or when combined with another paleothermometry method (e.g., δ 18 O, Epstein et al., 1953b).Nevertheless, the use of nonlinear methods in general allows discovering less obvious (nonlinear) relationships between proxies and temperature.Consequently, it is possible that the use of modern nonlinear multivariate statistics (among which the WDMR method) will help to find new proxies with hidden paleothermometer potential.The use of nonlinear models in general will probably open new research paths in paleoclimatology.
Figure 7 clearly shows that models based on a combination of proxies perform better than single proxy models.But it is also clear that not all proxy combination perform as well.Table 1 gives an objective overview of the contributions of Mg/Ca, Ba/Ca, Sr/Ca and Pb/Ca to paleotemperature models.It thus appears that Mg, already known as a temperature proxy (Klein et al., 1996b;Wanamaker et al., 2006), shows the highest contribution to the temperature reconstruction.More surprising is that Ba and Pb, which have not been proposed as temperature proxies, seem to contribute more to the temperature reconstruction than Sr which has been suggested as paleothermometer (Wanamaker et al., 2008).

Magnesium
Our results confirm the paleothermometry capacity of the Mg/Ca ratio as reported for several bivalve species by others (e.g., Klein et al., 1996b;Wanamaker et al., 2006).
However, Fig. 2 clearly shows that the Mg-temperature relationship is not linear.The Mg-temperature relationship seems to reflect that Mg incorporation in M. edulis is driven by a physiological response to temperature, with a maximal Mg incorporation around 16 • C. Except for the work of Vander Putten et al. (2000), which is based on the same dataset, a similar Mg/Ca-temperature relationship showing maximal Mg uptake at an intermediate temperature has not been reported in literature.Most published papers propose linear Mg-temperature relations for bivalves (e.g.Richardson et al., 2004;Pearce and Mann, 2006;Klein et al., 1996b, a).Freitas et al. (2006) observe an exponential Mg-temperature relationship for different bivalve species.That relationship is similar to the abiogenic Mg/Ca-temperature relationship reported by Oomori et al. (1987) and the temperature dependent Mg-incorporation in foraminifera reported by Barker et al. (2005).
On the other hand, it has been shown that Mg/Ca ratios in shells are influenced by growth rate (Ford et al., 2008) and by metabolic activity (Strasser et al., 2008).Moreover, Mg is shown to be incorporated for at least 33% in shell organic matrix (Foster et al., 2009;Takesue et al., 2008).Such biological controls on Mg may explain why combining Mg with other proxies results in better reconstructions (see Fig. 7).For instance if Mg incorporation in the shell depends indeed on physiology it is reasonable to assume that Mg incorporation will also be influenced by other vital factors, since the animal's physiological condition will be influenced by metabolic activity, growth rate, food availability and/or ontogenetic stage.Therefore, Sr (being a potential proxy for metabolic activity and growth rate), Ba (being a potential proxy for food availability) and Pb (also being influenced by ontogenetic stage) may explain some of the variation in the Mg/Ca profile of a shell.

Barium
Except for the specimens from Breskens, the nonlinear Bamodel results in fairly good SST reconstructions, indicating that Ba uptake in the shell of M. edulis is partly driven by temperature.It is probable that the Ba-temperature relationship is indirect and rather reflects temperature driven plankton blooming or water mixing events (Lazareth et al., 2003;Barats et al., 2009).These indirect relationships can be informative but one should be aware of the model errors that could be created, possibly biasing the temperature reconstruction.Indeed bloom events are quite complex and are influenced by many environmental parameters such as river discharge, wind speed, insulation etc (Cloern et al., 1995).The failure of the Ba-model at the Breskens site is probably due to such model errors.Indeed Fig. 3 shows that Breskens is the only site where a second Ba-peak is observed, although the temperature profiles at the three study sites monitored in 1996 are very similar.The model trained on shells of Terneuzen incorrectly couples the barium peak to temperature increase since all training shells independently showed a Ba peak coinciding with temperature increase in spring.As a result the model provides a similar interpretation for the second Ba peak observed for the Breskens shells although the origin of this second Ba peak is probably different.To avoid this kind of model errors it is not recommended to use Ba/Ca ratios as stand alone temperature proxy.
However, this does not mean than Ba/Ca ratios can not add information into a multi-proxy model.Several studies report that phytoplankton bloom events can influence the metabolism of the filter feeding bivalve, thereby inducing variation in shell growth rate (Versteegh, 2009;Schöne et al., 2006;Gillikin et al., 2008).Therefore, it can be expected that the combination of Ba and Sr (a potential proxy for shell growth and metabolism) in a multi-proxy model will contribute to resolving variations in other proxies which are due to shell growth.

Strontium
Contrasting with the studies that report a relation between Sr/Ca ratios and water temperature in calcitic bivalve shells (Carre et al., 2006;Freitas et al., 2005;Wanamaker et al., 2008)

Lead
The Pb/Ca-model does not result in accurate temperature reconstructions (Fig. 4) and when adding Pb to a multiproxy model, low or negative impacts are observed.This means that Pb uptake is poorly influenced by temperature and also that the variation in Pb/Ca ratios do not contain much information helping to explain variation in other proxies.Nevertheless, when Pb/Ca is added to the MgBa model, contribution factors increase (Table 1).This suggests that Pb/Ca and Ba/Ca are influenced by a common parameter.A common forcing for Pb/Ca and Ba/Ca, however, has not been reported in literature.However, it has been shown that Pb incorporation in a bivalve shell is influenced by temperature (Strasser et. al., 2008) and Pb/Ca profiles sometimes show ontogenetic trends (Dick et al., 2007).These facts may explain the positive contribution factors of Pb in the multiproxy models.However, Pb/Ca ratios in shells have been shown to be strongly influenced by anthropogenic activities (Gillikin et al., 2005;Richardson, 2001) rather than natural climate related changes.So, we do not recommend including Pb in a multi-proxy model.

Year to year and site specific variations
Several studies reveal that trace element profiles in shells may vary significantly between successive years (Barats et al., 2009) and between different study sites (Gillikin et al., 2006a).Our study as well reveals year to year and site specific variations (see Fig. 3).However, the accurate temperature reconstruction based on the shell from Knokke sampled during a different year, at a different site relative to the training site suggests that the models are relatively robust against year to year and site specific variations in trace element composition.Moreover, even though the distance between Terneuzen and Knokke is not more than 40 km the two sites strongly differ in environmental conditions: the Terneuzen site is an estuarine environment with a lower salinity compared to the Knokke site which is a more saline costal environment, therefore the model is probably fit for application to a wider environment than studied here.
However, we did observe differences in site specificity for different proxy combinations (e.g. the Ba/Ca problem that is observed for the Breskens site) and therefore the site specificity of every proxy has to be investigated.This needs to be done independently and in combination with other proxies before a model based on a specific proxy combination can be extrapolated to a broader environment.

Species specificity
The model presented in this paper is trained on M. edulis shells.Although we do not expect this model to be directly applicable to other species because Mg/Ca (the main player in the temperature reconstruction) is assumed to be driven by a physiological temperature response that is probably species specific, some preliminary tests suggest that the models may be extrapolated to other bivalves with calcitic shells.Moreover it should be possible to generate a specific WDMR model for other substrates such also corals, trees, sediments, etc.
Thus WDRM method could be used to develop nonlinear models to reconstruct the paleoenvironment for all different types of natural climate achieves.

Building new models using the WDMR method
As mentioned before, the model presented in this paper is species specific, implying that a different model needs to be constructed for other species.Furthermore, we believe that the WDMR method could also be used to build a stronger model for M. edulis shells.The model presented in this paper is based on trace elements of which some have never been linked to temperature before (i.e.Ba/Ca and Pb/Ca).Although Ba/Ca has clearly been shown to improve temperature reconstructions, a multi-proxy model that uses even more proxies with paleothermometry capacity would significantly improve the temperature reconstructions.Therefore, we encourage the construction of a WDMR model using high resolution measurements of Li/Ca ratios (Thebault et al., 2009), deuterium (Carroll et al., 2006) and oxygen isotopes (Epstein et al., 1953b).On the other hand we also encourage exploring other elemental and isotopic measurements using the WDMR-method since this method is able to detect less straightforward relationships between a potential proxy and its environment.The WDMR toolbox for Matlab is added as supplementary material to this paper.

The benefit of nonlinear methods
In this paper we observed that some proxy environmental relationships are substantially nonlinear and we demonstrate that using a nonlinear model to describe a proxy data set can improve temperature reconstruction performance with more than 1 • C compared to classical multiple regression techniques.

The benefit of combining proxies
Furthermore, we demonstrate that combining different proxies results in better temperature reconstructions.However, it is clear that not all proxies contribute equally to the final result.Our tests confirm that the Mg/Ca ratio in bivalve shells is a successful paleothermometer.We suggest that the Mg biomineralization is driven by a physiological response to changing temperature, which is possibly perturbed by metabolic activity and variable growth rate of the bivalve.The Combination of Mg, Ba and Sr into a multi-proxy model was successful because Ba and Sr reduce interfering effects due to metabolism and growth rate variation, thereby reducing the variance of the temperature prediction based on Mg.

The robustness of the WDMR method
The nonlinear multi-proxy model obtained by the WDMR is able to reconstruct temperature with a RMSE of less than 2.19 • C for a salinity ranging from 32 to 15.In comparison with other paleothermometry methods the performance using WDMR is good, if not better.This stresses that there is indeed a significant underlying low-dimensional structure in the proxy space.Although WDMR is a complex and sophisticated method, its success and robustness relies on its capability to nonlinearly combine proxy measurements into a multi-proxy model.One of the main messages of this contribution is therefore to encourage other researchers to combine their proxy measurements in one nonlinear multiproxy model, since this will allow identifying new proxies with paleothermometer potential.

Fig. 1 .
Fig. 1.Geographical position of the study sites in the Scheldt Estuary.Boxed: six shells from Terneuzen were used for training the models.Circled: the shells used for validation.

Fig. 2 .
Fig. 2. Left: Ba/Ca, Sr/Ca, Mg/Ca and Pb/Ca ratios plotted against water temperature (Vander Putten et al., 2000).Right: Ba/Ca, Sr/Ca and Mg/Ca concentrations plotted against each other.The shown curve indicates how the concentrations are changing with water temperature (in colour).

Fig. 3 .
Fig. 3.Chemical signature along the growth axis of the shells used to train (first column) and to validate (columns 2 to 5) the models.The trace element/Ca ratios are in mmol/mol.

Fig. 4 .Fig. 5 .
Fig. 4. Detailed visualization of the temperature reconstructions for the shell from the Knokke study site for all proxy combinations, the x-axis corresponds with sample number along the along the shell's growth axis.

Fig. 6 .
Fig. 6.Measured water temperature (dashed line) and reconstructed water temperature (solid line) obtained by the four-proxy WDMR model trained on Terneuzen data set and tested on a single validation shell from Terneuzen, Breskens, Ossenisse and Knokke.

Fig. 7 .
Fig. 7. Model reconstruction performance expressed as RMSE obtained on the validation data for the four sites, as based on the Terneuzen training dataset using different proxy combinations.

Table 1 .
Seven unique contribution factors are defined per proxy.Every contribution factor is defined by the difference between the RMSE of a model based on a proxy combination whith the investigated proxy and the RMSE of a model based on a proxy combination without the investigated proxy.The average contribution of each proxy is given per study site and for the total validation set.Negative contribution factors are marked in red and mean that the corresponding proxy does not contribute to a better reconstruction.
our results indicate that Sr/Ca ratios do not carry www.geosci-model-dev.net/3/653/2010/Geosci.Model Dev., 3, 653-667, 2010 much temperature information.The Sr/Ca-model computed in this paper does not result in satisfactory temperature reconstructions, neither for shells from the Terneuzen training site nor for the other sites.Moreover, when Sr is added to a multi-proxy model often a negative impact is seen, indicating that Sr uptake is poorly influenced by temperature and also that the variations in Sr/Ca ratios do not contain significant information that assists in resolving the variation in other proxies.Nevertheless, Sr/Ca seems to have a positive influence on the site specificity of Ba/Ca, suggesting that Ba/Ca ratios and Sr/Ca ratios are influenced by a common environmental factor.Lazareth et al. (2003) also observed some Sr/Ca maxima to coincidence with Ba/Ca-peaks.It is possible that the incorporation of both elements is influenced by shell growth rate as suggested at least for Sr/Ca by Carré et al. (2005), Gillikin et al. (2005) and Foster et al. (2009).Therefore, considering that Mg/Ca is a potential temperature proxy, even though it appears affected by variable shell growth rate and metabolic activity (Takesue et al., 2008), the combination of Mg/Ca with Sr/Ca and Ba/Ca can help explain a considerable fraction of the Mg/Ca signal noise.This is indeed observed in our dataset where the RMSE of the MgSrBa model is a significantly lower than the RMSE of the Mg/Ca model (RMSE (MgSrBa−model) -RMSE (Mg−model) = 1.28; 0.30; 0.40 and 0.34 for Terneuzen, Breskens, Ossenisse and Knokke, respectively).