Pseudo-proxy evaluation of Climate Field Reconstruction methods of North Atlantic climate based on an annually resolved marine proxy network

Two statistical methods are tested to reconstruct the inter-annual variations of past sea surface temperatures (SSTs) of the North Atlantic (NA) Ocean over the past millennium, based on annually resolved and absolutely dated marine proxy records of the bivalve mollusk Arctica islandica. The methods are tested in a pseudo-proxy experiment (PPE) set-up using state-of-theart climate models (CMIP5 Earth System Models) and reanalysis data from the COBE2 SST data set. The methods were applied 10 in the virtual reality provided by global climate simulations and reanalysis data to reconstruct the past NA SSTs, using pseudoproxy records that mimic the statistical characteristics and network of Arctica islandica. The multivariate linear regression methods evaluated here are Principal Component Regression and Canonical Correlation Analysis. Differences in the skill of the Climate Field Reconstruction (CFR) are assessed according to different calibration periods and different proxy locations within the NA basin. The choice of the climate model used as surrogate reality in the PPE has a more profound effect on the CFR skill 15 than the calibration period and the statistical reconstruction method. The differences between the two methods are clearer for the MPI-ESM model, due to its higher spatial resolution in the NA basin. The pseudo-proxy results of the CCSM4 model are closer to the pseudo-proxy results based on the reanalysis data set COBE2. The addition of noise in the pseudo-proxies is important for the evaluation of the methods, as more spatial differences in the reconstruction skill are revealed. More profound differences between methods are obtained when the number of proxy records is smaller than five, making the Principal Component 20 Regression a more appropriate method in this case. Despite the differences, the results show that the marine network of Arctica islandica can be used to skilfully reconstruct the spatial patterns of SSTs at the eastern NA basin.


Introduction
Several studies have targeted to reconstruct hemispheric or global average temperature from networks of proxy records (Mann and Jones, 2003;Marcott et al., 2013;Moberg et al., 2005), as well as spatial patterns of past temperature changes at global 25 (Rutherford et al., 2005;Wahl and Ammann, 2007) and regional scale (Ahmed et al., 2013;Luterbacher et al., 2004; during the last millennium. Most of these studies have primarily used terrestrial proxy records including very few marine proxies that could contain information about the vast ocean areas, e.g. over the North Pacific and the North Atlantic Ocean. In the context of marine proxy networks, reconstructions of global SSTs for the past 2000 years derived from individual marine reconstructions were recently published as a global synthesis of SST for the Common Era (CE), .
The performance of CFRs has been tested using PPEs by a number of studies (Rutherford et al., 2003;Von Storch et al., 2004).
Given the importance of the spatial information estimated in CFRs, a growing number of studies has explicitly evaluated the spatial performance of CFRs (Dannenberg and Wise, 2013;Evans et al., 2014;Li and Smerdon, 2012;Smerdon et al., 2008;Wang et al., 2014). However, only a few studies have tested the differences in the spatial performance of CRF methods 80 according to the modelled climate that forms the basis of the PPE (Mann et al., 2007;Smerdon et al., 2016;Smerdon et al., 2011). These studies tested the spatial skill of CFR methods using information from a composite proxy network including mostly terrestrial proxies. As CFRs depend on the characteristics of the proxy network used, such as proxy temporal resolution, growth season, character and level of noise, in our analysis we test CFR methods in the context of the annually resolved marine proxy network of Arctica islandica. The multivariate linear regression techniques tested are Principal Component Regression (PCR) 85 and Canonical Correlation analysis (CCA). PCR has been previously used to reconstruct climate using only marine proxy records (Evans et al., 2002;Marchal et al., 2002), while both methods have been commonly applied in the context of CFRs using annually resolved proxies (Gómez-Navarro et al., 2015;Smerdon et al., 2016;Smerdon et al., 2011). A fundamental assumption of PPEs is that the models can realistically simulate the spatiotemporal characteristics of the observed climate. Therefore, the models that are used in this analysis were chosen based on their ability to simulate the spatiotemporal characteristics of the 90 observed NA SSTs and have been previously evaluated in the context of paleoclimate reconstructions (Pyrina et al., 2017).

Data and Methodology
We performed a PPE using modelled and reanalyzed grid point SSTs co-located with proxy sites of Arctica islandica. The grid point SSTs were taken from the General Circulation Models (GCMs) CCSM4 and MPI-ESM-P, as well as the centennial in-situ observation-based estimate of SSTs, COBE2 (Hirahara et al., 2014). The SSTs are analysed for the summer period (June-95 August) and for the NA region between 60˚ W -30˚ E and 40˚ N -75˚ N, motivated by the growing season of the bivalve shell Arctica islandica (Schöne et al., 2004;Schöne et al., 2005). The pseudo-proxy sample sites are based on five real world proxy sites of Arctica islandica including collection sites in the North Sea (NS: 1˚E, 58.5˚N) (Witbaard et al., 1997), the Irish Sea (IrS: 5˚W, 52.5˚N) (Butler et al., 2009), the coast of Scotland (Sct: 7˚W, 56.5˚N) (Reynolds et al., 2013), the North Icelandic Shelf (IS: 20˚W, 66.5˚N) (Butler et al., 2013) and a location at Ingoya Island (InI: 24˚E, 71.5˚N) (Mette et al., 2015).

Observational and Proxy data
We used the spatially interpolated reanalysis data set COBE2 (Hirahara et al., 2014). The COBE2 data set was developed by the Japanese Meteorological Agency, covers the period 1850-2013 AD and has a spatial resolution of 1˚x1˚ lat x lon. COBE2 data pass first quality control using combined a-priori thresholds and nearby observations and are later gridded using optimal 105 interpolation. Data up to 1941 were bias-adjusted using "bucket correction" (Hirahara et al. 2014). COBE2 combines SST measurements from the release 2.0 of the International Comprehensive Ocean-Atmosphere Data Set, ICOADS (Worley et al., 2005), the Japanese Kobe collection, and readings from ships and buoys. The proxy data of Arctica islandica used in this analysis were downloaded from NOAA's (National Oceanic and Atmospheric Administration) National Centers for Environmental Information (NCEI, https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets) and refer to the 1357-year Arctica 110 islandica reconstructed chronology from Butler et al., 2013. The multi-centennial absolutely dated chronology was reconstructed using annual growth increments in the shell of Arctica islandica and spans the period from 649 AD to 2005 AD.

Models
Two models are employed in this study, the CCSM4 model and the MPI-ESM-P model, which are part of the 5th phase of the Climate Model Intercomparison Project (CMIP5/ http://cmip-pcmdi.llnl.gov/cmip5/). The models' original output was re-115 processed and re-gridded to a regular grid for subsequent comparisons with the COBE2 data set. Therefore, the output was regridded onto a 1˚×1˚ horizontal resolution. The output of the models used in this study includes the combination of the past1000 runs and the historical runs, so that the period used spans from 850 AD to 1999 AD.
The CCSM4 (Gent et al., 2011) uses the atmosphere component Community Atmosphere Model, version 4 (CAM4) (Neale et 120 al., 2013) and the land component Community Land Model, version 4 (CLM4) (Lawrence et al., 2012). Both components share the same horizontal grid (0.9˚ latitude × 1.25˚ longitude). The CCSM4 ocean component model (POP2) is based on the "Parallel Ocean Program", version 2 (Smith et al., 2010). The ocean grid has 320 × 384 points with nominally 1˚ resolution except near the equator, where the latitudinal resolution becomes finer to better simulate ENSO dynamics, as described in (Danabasoglu et al., 2006). CICE4, the CCSM4 sea ice component model is based on version 4 of the Los Alamos National Laboratory 125 "Community Ice Code" sea ice model (Hunke et al., 2008). The atmosphere, land, and sea ice components exchange both state information and fluxes through the coupler for every atmospheric time step. The fluxes between atmosphere and ocean are calculated in the coupler and communicated to the ocean component only once a day. The CCSM4 simulation starts at 850 AD and continues to 1850 AD, where it matches up and is extended as an additional ensemble member of the CCSM4 twentieth century simulations that ends in December 2005 (Landrum et al., 2013

140
In MPI-ESM-P the atmosphere model ECHAM6 (Stevens et al., 2013) was integrated using a horizontal resolution of spectral truncation T63 (1.875), while the ocean/sea-ice model MPIOM (Marsland et al., 2003) features a conformal mapping grid with nominal 1.5° resolution. There is one grid pole over Antarctica and one grid pole over Greenland, which leads to considerably higher resolution in the NA. For land and vegetation the component JSBACH (Reick et al., 2013) is used and for the marine biogeochemistry the HAMOCC5 (Ilyina et al., 2013). The coupling at the interfaces between atmosphere and land 145 processes, and between atmosphere and sea ice occurs at the atmospheric time step, which is also the time step of the land processes, except for the dynamic vegetation, which is updated once a year. The coupling between atmosphere and ocean as well as land and ocean occurs once a day. In the past1000 simulations a prescribed CO 2 forcing is used. For volcanic aerosol optical depth and effective radius the Crowley and Unterman reconstruction (Crowley and Unterman, 2013)  In this work we used three experiments of the MPI-ESM model. The r1 and r2 experiments were initialized with the same ocean state, but they differ in the standard deviation of the assumed lognormal distribution of the volcanic aerosol size, while the simulations r2 and r3 use the same parameter setting but are started from different initial conditions (Jungclaus et al., 2014). For 155 the historical simulations the applied boundary conditions follow the CMIP5 protocol, except for land-cover-changes, where the Pongratz et al. 2008 data set is used.

Methodology
Based on five proxy locations of Arctica islandica we reconstructed the summer SST evolution of the NA region, during the industrial (1850-1999 AD) period of the last millennium, using two CFR methods. In both approaches the goal is to reconstruct  (1950( -1999. The NA fields' inter-annual anomalies were 165 reconstructed using the COBE2 reanalyzed SSTs and the CMIP5 modelled SSTs. To conclude whether the usage of different models or of reanalysis data has an effect on the reconstruction we correlated the original modelled or reanalyzed inter-annual anomalies of the NA field with the reconstructed modelled or reanalyzed inter-annual detrended anomalies, respectively. In addition, we compared the temporal variances at grid-cell scale of the reconstructed and modelled SSTs.

170
The first step of the PPE is the estimation of the NA SST field co-variances for the different calibration periods using Principal Component Analysis (PCA), (Eq. 1). Each eigenvector is associated with a spatial pattern (EOF, Empirical Orthogonal Function) and its temporal evolution (PC, Principal Component). In Eq. (1), x ⃗ t is the field vector of the NA SST anomalies and i the number of eigenvectors. In our analysis we kept the first 10 eigenvectors, as they represent more than 90% of variability. The time, t, depends on the calibration period that we refer to.
For the five sampled SST time series (Proxy j,t with j representing the respective proxy location) we calibrated our regression model (Eq. 2) against the PCs estimated during the calibration period, and predicted the calibration coefficients, α, using PCR.
In our approach we assume that the PCs are linearly related to the pseudo-proxies, so that they represent large-scale climate 180 variations for the SST fields. This relationship is modelled through a disturbance term or error variable ε. The error could be an unobserved random variable that adds noise to the linear relationship between the dependent variable (PC) and the regressors (Proxy SSTs). Based on PCR (Eq. 3) we then predict the principal components PĈ i,t during the reconstruction period, assuming that the calibration coefficients calculated in Eq.
(2), are stationary in time. In this case the time, t, depends on the reconstruction period that we refer to. i=1 Assuming that the dominant patterns of climate variability are similar in recent and past centuries, we predict the x ⃗ t field vector of the NA SST anomalies for the reconstruction period, using the predicted PĈ i,t and the EOF ⃗⃗⃗⃗⃗⃗⃗⃗ i patterns calculated in Eq. (1). This stationarity assumption holds at least for multi-decadal timescales and allows us to deduce back in time the surface temperature patterns (Mann et al., 1998).

Canonical Correlation Analysis
The first step of the PPE using CCA is the eigenvalue decomposition and subsequent truncation of the NA SST field and the proxy SST field (Eq. 1), during the calibration interval. The key feature of this analysis involves concatenating the pseudo-proxy SST time series, or in other words the predictor data, so that we can define both the time and space evolution of the climate system. After the transformation of the data to EOF coordinates we retained five empirical orthogonal functions, for each field,

195
that were subsequently used to calculate the pairs of Canonical Correlation Patterns (CCP NA , CCP pr ) and their time depended Canonical Coefficients (CC NA , CC pr ), as shown in Eq. (4) and Eq. (5) for the NA SST field and the proxy SST field, respectively.
For each estimated CCP the correlation between the respective CCs is maximized and not correlated with the CCs of another pair of patterns. Therefore, the CCs fulfill the condition of orthogonality.
The steps following, are the same as the ones used in the method PCR, but instead of regressing and predicting the PCs we use the CC i,t pr . To reconstruct the x ⃗ t field vector of the NA SST anomalies we used the predicted CĈ i,t pr and the CCP ⃗⃗⃗⃗⃗⃗⃗ i NA patterns calculated in Eq. (4).

205
The PPEs were conducted using idealized pseudo-proxies and noise-contaminated pseudo-proxies. Idealized pseudo-proxies are the raw grid point SSTs, from simulations or reanalysis, co-located with the collection sites of the bivalve shell Arctica islandica, while in the case of the noise-contaminated pseudo-proxies the grid point temperatures are deteriorated by adding statistical noise in order to mimic more realistically the real world proxy records of Arctica islandica. The exact level of non-climatic noise in the real proxies is usually not known and could be strongly dependent on the nature of the proxy indicator (von Storch et al., The dynamics of many physical processes can be approximated by first or second order ordinary linear differential equations, whose discretized versions can be represented by autoregressive processes (Von Storch and Zwiers, 2001). An autoregressive processes of order k=0, where k is the time lag, is white noise (Z t ). An autoregressive processes of order k=1, or AR(1),

215
represents a discretized first order linear differential equation and can be written as: where α 1 is the damping coefficient and Z t represents a random variable uncorrelated in time. The Yule-Walker equations can be used to derive the first k+1 elements of the autocorrelation function and for an AR(1) process they give α 1 =ρ 1 , where ρ k is the autocorrelation for lag k. For a positive damping coefficient, an AR(1) process is unable to oscillate. Its 'spectral peak' is located 220 at frequency ω=0 and therefore the variation X t behaves as red noise. Furthermore, the stationarity condition for an AR (1) process implies that |α 1 |<1. Therefore, to approximate the variation of noise in Arctica islandica the autocorrelation function of an Arctica islandica chronology located at the Icelandic shelf was calculated and found equal to 0.4 at lag 1 year. Assuming that the relative amount of noise is constant for all the noise-contaminated pseudo-proxies, red noise pseudo-proxies were constructed by fitting the parameters of the AR(1) process to the simulated grid-point temperatures.

Ideal pseudo-proxies
In this section the results are based on five ideal pseudo-proxies co-located to Arctica islandica sites that are "sampled" from the SST output fields of three realizations of the MPI-ESM model, the CCSM4 model and the COBE2 data set. The correlation between the reconstructed and the original SST-anomaly evolution of the NA field is calculated for the industrial reconstruction  5S. Moreover, the ratio of the estimated Standard Deviation (SD) according to the reconstructed SST-anomaly evolution of the NA field and the estimated SD according to the original SST-anomaly evolution was calculated (SD reconstruction /SD original ) and shown in Fig. 1 and Fig. 2 for the recent and LIA calibration periods, respectively.

240
Looking at the results of the same model simulation and method but for different calibration periods we can see that the spatial skill of the reconstruction does not significantly change according to the calibration period of the regression models. Therefore, the calibration coefficients can be considered stationary. Additionally, the SST patterns shown by the correlation maps do not profoundly change when a different method is applied, indicating that with both methods the SST evolution of the eastern NA basin can be reconstructed using only five proxy locations of Arctica islandica. More pronounced differences are found between 245 the pseudo-reconstructions conducted with each model and also with the pseudo-reconstruction using the reanalysis. For instance, when the industrial period is reconstructed using the recent calibration period ( Fig. 1) according to CCSM4 and COBE2, the SSTs between Iceland and Norway can be skilfully reconstructed with a correlation coefficient approximately equal to r≈+0.8, but according to MPI-ESM the skill drops for all realizations (Appendix, Fig. 3S). Another interesting finding is that the SSTs along the south-east coast of Greenland can be reconstructed with a good skill using the CCSM4 model (r≈+0.9), while 250 the reconstruction skill drops for the COBE2 data (r≈+0.7) and the MPI-ESM model (r≈+0.5). The correlation maps provided from the CCSM4 model and the COBE2 data are more homogenous than the ones given by the MPI-ESM-P model. The areas that exhibit correlation values higher than r≈+0.6 are more widespread in CCSM4 and COBE2 compared to MPI-ESM-P. The standard deviation ratios indicate that variance is mostly preserved in areas where field correlations are high. According to both models and the COBE2 reanalysis data, PCA results to a loss in variance in most of the NA basin, with the largest decrease 255 occurring over the West Atlantic basin. Generally, the CCA derived CFRs are also afflicted by a loss in variance, but enhanced variance is found over the West and South NA basin.

Noise contaminated pseudo-proxies
The results of this section are illustrated as in section 3.1, but the calculations were performed with noise-contaminated pseudoproxies. The methods calibrated during the recent period and LIA are shown in Fig. 3 and Fig. 4, respectively, for the CCSM4 meaning. According to the noise contaminated experiment, the standard deviation ratios indicate that the CFRs conducted using either the PCA or the CCA method are afflicted with a loss in variance over the whole NA basin.

Test of proxy locations
In the following chapter we test the contribution of the different Arctica islandica proxy sites on the reconstruction skill. We  NP=3). The reconstruction is performed with both statistical methods, using ideal proxies (Appendix, Fig. 6S) and noise contaminated proxies (Appendix, Fig. 7S). In Fig. 5 the results using three Arctica islandica sites are shown for both ideal proxies and noise contaminated proxies.

285
Regarding the ideal PPE (Appendix, Fig. 6S), only the results of the r2 realization of the MPI-ESM-P model indicate that the PCA method leads to a greater reconstruction skill than the respective results according to the CCA method. That applies in both cases, for two and three proxy locations. Moreover, anti-correlation between the reconstructed and original SST evolution occurs with the CCA method, mostly when only two proxy sites are used. Additionally, the differences amongst models and amongst 290 models and data are more profound for the CCA method than for the PCA, because the number of CCPs that can be used for the CFR is limited by the number of proxies used (two CCPs for NP=2 and three CCPs for NP=3). Comparing the results of Fig. 5, regarding the reconstruction using three Arctica sites, to the corresponding results of Fig. 1 that are based on the reconstruction using five Arctica sites, we see that in most cases the sites in the Irish Sea and the off the coast of Scotland do not only contribute to the reconstruction of the SSTs on their surrounding waters but they additionally increase the reconstruction skill over the east

Stationarity of the calibration coefficients
Both CFR methods showed that the spatial skill of the reconstruction does not profoundly change according to the period during which the regression models were calibrated, even when the calibration period chosen relates to the recent period, a period dominated by anthropogenic forcing. This result is robust amongst models and amongst models and data and could be an describing an unforced climate. In contrast to those previous results, these authors found that the simulated SST modes, during such long time intervals, describe regional SST variability patterns in the NA consistent with those simulated and observed over 320 the last 160 years.
Another reason that could potentially explain the similarities of the reconstruction skill amongst different calibration periods could be that the relationship of the SSTs of the regions where the proxies are located and the modes of variability that influence the NA basin remain unchanged. Lohmann et al., 2005 found that during the period 1900-1998 AD the Arctic Oscillation-related 325 temperature teleconnections, show weak decadal variations in some regions of the NA. The Arctic Oscillation signature in climate variables was detectable also during the spring season, which is of practical relevance as the climate information obtained from most terrestrial and marine proxy archives is more linked to the growing season rather than to winter (Cook et al., 2002). The NA regions that Lohmann et al., 2005 identified with stable teleconnections include almost all the proxy locations used in our study (see in Lohmann et al., 2005; Fig. 4), but their analysis regards only the recent period and decadal variations.

330
Regarding the North Atlantic Oscillation (NAO), variable teleconnections have been detected in coupled ocean-atmosphere model simulations (Raible et al., 2001;Zorita and González-Rouco, 2002). Moreover, model inaccuracies in the representation of the mean NA climate, NAO and other modes of climate variability that influence the NA SSTs could lead to a poor representation of the NA regions. However, as CFRs are covariance-based approaches, a poor representation would lead to a generally low reconstruction skill in the NA basin, which is not consistent with the results of this study. Therefore, a more 335 plausible reason that explains the similarities of the reconstruction skill amongst different calibration periods relates to the minimization of possibly existing non-stationarities by using proxy sites from multiple regions (Batehup et al., 2015).

Model and method dependent spatial skill
As presented above, more profound spatial differences in the reconstruction skill are found for the different model simulations, rather than for the different periods in the last millennium that were used to calibrate the regression models. The pseudo-proxy 340 results of the CCSM4 model are found to be closer to the pseudo-proxy results of the COBE2 data set, than the ones calculated based on the MPI-ESM-P model. In the study of Pyrina et al. 2017, and in

365
The differences of the reconstruction skill according to CCA or PCR can be better identified in the case of the MPI-ESM model, because of the spatially more heterogeneous correlations between the original and the reconstructed SST anomalies. For the reconstruction of the NA SSTs using the PCR method, we take advantage of the dominant patterns of SST variability of the NA basin, while with the usage of the CCA method we exploit only those patterns with time histories related to the time histories of the proxy network. In the case of five proxy locations, and even worse of three and two, we see that this method is problematic 370 because we lose information due to the limited number of CCP pr that can be derived from our small sized proxy network.
Moreover, the differences in the methods can be better seen when we use the noise contaminated pseudo-proxies, as in this case the loss of information is greater, because the noise contaminated pseudo-proxies do not explain 100% of the SST signal.
Regarding the ideal experiment, PCA generally leads to an underestimation of variance over the whole NA basin, while it contaminated experiment both methods lead to a loss in variance. The decrease in variance can be expected for regression-based CFR methods that blend signal and error variances as a characteristic of formulation (e.g. von Storch et al., 2004). The source of a loss in variance is likely associated with the manner in which the eigenvalue spectra are truncated by the two methods 380 (Smerdon et al., 2010).

Conclusions
From the PPEs derived herein we have demonstrated that a small sized proxy network of the marine bivalve mollusk Arctica islandica (five proxy locations) can produce skilful spatial SST reconstructions for the eastern NA basin. The tested CFR methods can alter the spatial skill of the reconstruction especially when noise contaminated pseudo-proxies are used. Therefore,

385
it is important to assess CFR methods with noise-contaminated proxies, rather than with ideal ones. Moreover, the CFRs conducted by both methods suffer variance losses. Even though the CCA method is problematic when a significantly low number of proxies is used (two and three proxies), the spatial skill of the reconstruction using CCA and five proxy locations is similar to the results calculated with the PCR method. The calibration period does not significantly affect the reconstruction, while the most profound changes in the spatial skill of the reconstruction are caused by the model simulation used.

Author Contributions
The analysis was performed by M. Pyrina with the consultation of S. Wagner and E. Zorita. M. Pyrina prepared the manuscript with contributions of co-authors.  1 and 4) and SD ratio (columns 2 and 4) between the reconstructed and the original SST-anomaly 600 evolution of the NA field during the industrial era. The results are given for the recent calibration period (1950( -1999 and for the two different reconstruction methods (1 st and 2 nd column CCA, 3 rd and 4 th column PCA) for the models CCSM4 (first row) and MPI-ESM-P r1 (second row) and the reanalysis data COBE2 (third row).

605
The results are given for the LIA calibration period (1650-1699 AD) and for the two different reconstruction methods (1 st and 2 nd column CCA, 3 rd and 4 th column PCA) for the models CCSM4 (first row) and MPI-ESM-P r1 (second row).
Figure 3: Same as in Fig. 1, but for the noise-contaminated pseudo-proxies.

610
Clim . Past Discuss., doi:10.5194/cp-2017-61, 2017 Manuscript under review for journal Clim. Past Discussion started: 12 April 2017 c Author(s) 2017. CC-BY 3.0 License. Figure 4: Same as in Fig. 2, but for the noise-contaminated pseudo-proxies. Figure 5: Correlation coefficient between the reconstructed and the original SST-anomaly evolution of the NA field during the industrial era for the ideal (columns 1 and 3) and the noise contaminated experiment (columns 2 and 4) when the number of proxy locations used is NP=3. The results are given for the recent calibration period (1950( -1999, for two different reconstruction methods (1 st and 2 nd column CCA, 3 rd and 615 4 th column PCA) and for the models CCSM4 (first row) and MPI-ESM-P r1 (second row), and the reanalysis data COBE2 (third row). Figure 1S: Correlation coefficient between the reconstructed and the original SST-anomaly evolution of the NA field for the ideal (column 1 and 3) and the noise contaminated experiment (column 2 and 4). The results are given for the MCA calibration period (1000-1049 AD) and 620 for the two different reconstruction methods (1 st and 2 nd column CCA, 3rd and 4th column PCA) for the models CCSM4 (1st row) and three realizations of the MPI-ESM-P model (2 nd , 3 rd and 4 th row). Figure 2S: As in Fig. 1S, but the results are given for the LIA calibration period (1650-1699 AD).

APPENDIX figure captions
625 Figure 3S: As in Fig. 1S, but the results are given for the recent calibration period (1950( -1999. The last row contains the pseudo-proxy results of the COBE2 reanalysis data. Figure 4S: As in Fig. 3S, but the results are given for the industrial calibration period (1850-1999 AD).
630 Figure 5S: As in Fig. 3S, but the results are given for the pre-industrial calibration period (850-1849 AD). Figure 6S: Correlation coefficient between the reconstructed and the original SST-anomaly evolution of the NA field during the industrial era, when the number of proxy locations used is NP=2 (1 st and 3 rd column) and NP=3 (2 nd and 4 th column). The results are given for the recent calibration period (1950( -1999 and for the two different reconstruction methods (1 st and 2 nd column CCA, 3 rd and 4 th column PCA) for 635 the model CCSM4 (1 st row), three realizations of the MPI-ESM-P model (2 nd , 3 rd and 4 th row) and the COBE2 reanalysis data (5 th row). Figure 7S: As in Fig. 5S, but the results are given for the noise contaminated pseudo-proxy experiment. Figure 8S: SD ratio between the reconstructed and the original SST-anomaly evolution of the NA field during the industrial period, for the 640 ideal (column 1 and 2) and the noise contaminated (column 2 and 4) pseudo-proxy experiment. The results are given for the recent calibration and for the two different reconstruction methods (1 st and 2 n column CCA, 3 rd and 4 th column PCA) for the models CCSM4 (1 st row) and three realizations of the MPI-ESM-P model (2 nd , 3 rd and 4 th row), as well as the COBE2 data (5 th row). Clim. Past Discuss., doi:10.5194/cp-2017Discuss., doi:10.5194/cp- -61, 2017 Manuscript under review for journal Clim. Past Discussion started: 12 April 2017 c Author(s) 2017. CC-BY 3.0 License. Calibration period: 1950-1999AD Reconstruction period: 1850-1999 Figure 8S: SD ratio between the reconstructed and the original SST-anomaly evolution of the NA field during the industrial period, for the ideal (column 1 and 2) and the noise contaminated (column 2 and 4) pseudo-proxy experiment. The results are given for the recent calibration and for the two different reconstruction methods (1 st and 2 nd column CCA, 3rd and 4th column PCA) for the models CCSM4 (1 st row) and three realizations of the MPI-ESM-P model (2 nd , 3 rd and 4 th row), as well as the COBE2 data (5 th row).