Atlantic Multidecadal Variability from the Last Millennium Reanalysis

The Last Millennial Reanalysis (LMR) employs a data assimilation approach to reconstruct climate fields from annually-resolved proxy data over years 0-2000CE. We use the LMR to examine Atlantic Multidecadal Variability (AMV) over the last two millennia, and find several robust thermodynamic features associated with a positive Atlantic Multidecadal Oscillation (AMO) index that reveal a dynamically-consistent pattern of variability: the Atlantic and most continents warm; sea ice thins over the Arctic and retreats over the Greenland-Iceland-Norwegian Seas; and equatorial precipitation shifts northward. 5 The latter is consistent with anomalous southward energy transport mediated by the atmosphere. Net downward shortwave radiation increases at both the top-of-atmosphere and surface, indicating a decrease in planetary albedo, likely due to a decrease in low clouds. Heat is absorbed by the climate system and the oceans warm. Wavelet analysis of the AMV time series shows a reddening of the frequency spectrum at the 50-to-100-year time scale, but no evidence of a distinct multidecadal or centennial spectral peak. This latter result is insensitive to both choice of prior model and calibration dataset used in the data assimilation 10 algorithm, suggesting that the lack of a distinct multidecadal spectral peak is a robust result.

as it is made clear later, they refer to an an ensemble approach to do global past climate reanalyses.I would suggest to append "approach" to the current title.We have revised the title of the manuscript to "Insights into Atlantic Multidecadal Variability using the Last Millennium Reanalysis Framework."9. Does the abstract provide a concise and complete summary?From the abstract only, and even in the paper, it is unclear until page 3 in the manuscript whether the authors have done any new assimilation in this study, or used a previous reanalysis product (LMR) as documented in Hakim et al. (2016).Later it is understood that the use the same LMR approach as in Hakim et al.(2016), but the author re-do the assimilation step.This should be made more clear in the abstract.10.Is the overall presentation well structured and clear?Yes 11.Is the language fluent and precise?Yes 12. Are mathematical formulae, symbols, abbreviations, and units correctly defined and used?Yes, in general.
13. Should any parts of the paper (text, formulae, figures, tables) be clarified, reduced, combined, or eliminated?See specific comments below.
14. Are the number and quality of references appropriate?Basic references to the data assimilation method are needed, as it is the base for the reanalysis on which the study is built upon.
15.Is the amount and quality of supplementary material appropriate?The authors indicate the use multitaper spectra to evaluate if amplitudes in the wavelet analysis are higher than red noise at 95% confidence level.I would suggest to include this analysis as supplementary material (if not in the manuscript itself).We have added a Supplemental Information section that includes the multitaper spectra of the different reconstructed AMO indices described in the manuscript.

General comments:
Section 2 (Methods) is poor.The authors should be more explanatory about some aspects in the assimilation approach.The assimilation approach itself is unclear.For example, they do not mention the word "ensemble" at all in the manuscript.Thus it would appear they use a standard Kalman filter, which is unfeasible given the size of the problem.One need to resort to Hakim et al.(2016) too often to understand the basics of their assimilation approach.Thank you for this critical comment on the methods.We didn't want to burden this paper with a repeat of details presented elsewhere, but we realize now that the description of the method was too brief and as a result unclear.In response to this, and a comment from reviewer #2, we have added new material (p4, lines 3-7 and 29-30; p5 lines 3-33) to describe the data assimilation approach in words, and to map out some details of the calculation that are described subsequently.We have also expanded our description to include information on the ensemble approach to the prior, and the ensemble square-root solution method.We provide sufficient detail that a reader familiar with ensemble filters should be able to reproduce the results, and references for more details for readers less familiar with this approach.Hakim et al.(2016) use a square root Kalman Filter (Whitaker and Hamill, 2002), an ensemble approach, at the core of their assimilation method.If this paper also uses the same approach, the authors should at least include at this reference, indicate the ensemble size, and explain the ensemble generation method for the background fields.Lines p4.l23-24indicate they used CCSM4 from CMIP5 as background.How exactly they transformed a single simulation in a background ensemble?Although all these could potentially assumed similar to Hakim el al.(2016) it should be described here as well, as this is strongly related to the results.All of these issues are now clarified in the revised text as described above.Also the explanation of the forward model H would benefit from stating more clearly what is H in general and how it was explicitly obtained here (1st paragraph in Section 2.4 in Hakim et al (2016) explains this, but at least the should summarise it here).Also as Hakim et al. (2016) mention, temperature alone is not a good predictor of tree-ring widths everywhere.This, for example, will result in the linear model regression for the tree-ring width with temperature as the only predictor variable having very high regression residuals.There is not really anything terrible in this in the DA context, as by construction the R matrix originates not only from measurement errors (i.e.errors in the proxy data), but also from errors in the forward operator H, as is then this case.But please note down, as this is not evident for readers not so familiar with assimilation.As part of the revisions described below, we have added to the description of the observation model, and note the additional sources of error that contribute to matrix R (see p4, lines 30-31 and p5, lines 1-2).
I am also wondering about the general use of the linear model as H for all proxies.Have other relationships, apart from the linear one with temperature, been explored for H in the LMR?This could be e.g. a nonlinear proxy∼T relationship, or the possible use of other state variables (moisture) or geographical coordinates as predictor variables (e.g.latitude -see related note by the Editor-).We note here that recent (e.g., Dee et al. 2016; doi : 10.1002/2016MS000677) and ongoing research has been devoted to exploring the impact of bivariate and more complex observation models for trees, but that is beyond the scope of the present paper.
A note is the background state for the LMR as described in Hakim et al.(2016) is biased by design (as temporal trend are removed from the background state), and it is up to the proxy data to bring the re-analysed fields to somewhere in between the two sources of information (the background and the observations).Hence the resulting reanalysis will give, also by design, a "smoothed" version of the variability indicated by the proxy observations.The authors should comment on this, as this would be indeed a factor leading to reduced amplitudes in the higher frequencies of the resulting time series, which can also explain that the multidecadal variability is not so evident in this reanalysis.We have added a paragraph discussing this exact point.In particular, we note that while the reviewer's interpretation applies to a single proxy, the reconstruction at each location in fact depends on many proxies having unequal weighting.As a result, the temporal variability in the reconstructed time series at a point can be larger or smaller than that of a single record at that point.Validation of the results provides a way to know if the weighting, in aggregate, is appropriate.We have added several sentences to discuss how this method has been validated against instrumental data and independent (not assimilated) proxies.Please see p5, lines 3-12 and 19-33.Also, if I have understood it right, the regression in Fig. 2 (and related results) is based on the single member for CCSM while it is based on the mean of the ensemble for LMR.Right?Please, state more clearly.This is correct.We have clarified this in the Methods section.
A figure showing the mean and the variance of the ensemble background state (i.e.prior ccsm4) for SST and surface temperature over the world would be helpful (for comparison with Figure 2).We have included a new figure in the Supplemental Information section showing the mean and variance of the surface temperature and SST in the CCSM4.
On the other hand, I find the remaining Sections in the manuscript quite clear and adequate.Specific comments: p4.l5.The explanation for the calculation of the B matrix is certainly wrong "(computed as the sample mean of x_b)".Correct.Yes, thank you.This has been deleted.p4.l32.Clarify why you chose these parameters (specifically the 31 coefficients) for the Lanczos filter.Add references as needed.The Lanczos filter has superior spectral properties for the current application, and the coefficients are simple to compute.We have added a reference to Luchan (1979) for further information regarding the Lanczos filter and computing the filter weights (see p6, line 4).p5.l6.As the wavelet transform is a series of bandpass filters, it is better to give also here the range for the scales used in the transform, apart from spacing between scales (even if it can be understood form Figures later on).Also, indicate units.The range of scales (2 to 1000 yrs, with only 2 to 150 yrs shown) is now given in the Methods section (p6, line 9) and the units (in units of variance of the AMO index, K 2 ) are now included in the captions to Figures 9 and 10.
Figure 1.y-axis seems wrong in Fig. 1a.If not, explain units.We had (incorrectly) not removed the mean temperature from the area-averaged north Atlantic surface temperatures when plotting the AMO index in Figure 1a.We have corrected this in the revised Figure 1.p6.12.As mentioned, the weaker amplitude in the LMR with respect to the other re-analyses (Fig. 1c) is likely related to my previous general note about the background state generation, related to temporally de-trending the background.The authors should could comment on this, and in other related discussion throughout the paper.We have included this caveat in the revised manuscript (p5, lines 29-33): We note that because the prior is the same for every year, all temporal variability is determined by the proxies.As a result, the weighting in equation ( 1) involves a temporally invariant prior and temporally variable proxies so that, for any given proxy, reduced temporal variability may be expected.However, since the time series at any point in the reconstruction depends on the prior and many proxies, all having different errors, the variability at a proxy location may in fact be larger or smaller than that of a given proxy at that point.

Overview
In this paper, Singh and coauthors examine the characteristics of North Atlantic Mul-tidecadal SST variability in the Last Millennium Reanalysis (LMR).The LMR is a sim-ulation carried out with CCSM4 using data assimilation from proxies over the years 0-2000CE, and is mainly documented in a previous study.They find some known changes associated with the Atlantic Multidecadal Variability (AMV), e.g.: warming of the northern hemisphere, decrease in sea ice and northward shift in the ITCZ, and some unknown changes, in particular that there is a difference in the ocean heat trans-port between CCSM4 and the LMR.In addition, they carry out a wavelet analysis and show that there is no distinct peak at multidecadal timescales in the AMV index in the LMR, consistent with the notion that the AMV may be interpreted as a red noise process (which is in contrast to many modeling studies argue that there is multi-decadal peak in the AMV driven by oscillations in the AMOC).I think the results are interesting and fit in the broader discussion of the mechanisms and timescales of the AMV.The paper is well written, however I recommend some changes to improve the clarity in the presentation, which I elaborate in more detail below.

General comments
1.It is quite hard to understand for a general audience how this data assimilation works, how it compares with other studies, and in particular why there are so many differences with CCSM4 prior.When discussing ocean and heat transports I don't understand why LMR is so different from CCSM4.In general, I find it hard to follow because there is no comparison with mechanisms operating in the real world and would be easier to understand if it was compared with direct observations (even if we know that observations over the 20th century are affected by anthropogenic forcings).Thank you for this comment.We have rewritten portions of the method section to make the material more accessible for non specialists, and also to provide more details for readers that may wish to repeat these calculations (see p4 and p5 in the revised manuscript).As for why the LMR is different than the CCSM4, it derives from the information obtained from the proxies.As we describe in the revisions, all temporal variability comes from the proxies (not from the prior), which is one reason they are different.Another is that the proxies also change the spatial covariance information (i.e., the prior and posterior state-covariance matrices differ as a result of the information in the proxies).
Comprehensive analysis of the differences in dynamics is beyond the scope of this paper, but an area of current research.
2. In the paper AMV and AMO are used interchangeably, making it quite confusing to understand what the authors refer to.AMV is typically used to refer to internal+forced North Atlantic variability while AMO is used to refer to internal variability only (see Booth 2015).
We agree that AMV would imply both forced and unforced variability.We have removed most references to AMV, unless the reference is to both forced and unforced variability.
3. The AMV index is defined as the area average SST index from 120W to 0 and from 0 to 60N, and some studies are referenced including Clement et al. 2015. Clement et al. used 80W-0,0-60N as well as other studies which use 70/75W.In other words, 120W seems too far west, including a lot of not ocean (land) areas.The authors should check if their results are consistent when using the most commonly utilized index (80W-0,0-60N).The 120W is an unfortunate misprint in the manuscript.We computed the AMO index identically to that described in Clement et al (2015), with the western border of the bounding region at 80W, not 120W.We apologize for the confusion.
4. A lot of times the authors generally refer to similarities amongst figures in a very qualitative and general way.They have to be more quantitative and compute pattern correlations between each figure in order to draw any solid conclusion.
We have now included correlations to support our claims of similarity between different figures and time series as follows: (1) That the LMR reconstruction of the AMO index over the modern period (1850 to 2000) closely follows the AMO index derived from other reanalyses (GIStemp, ERA-20C, and 20CR).See p7, lines 13-14.(2) That the spatial pattern of the surface temperature anomaly associated with a positive AMO index (as computed from the LMR) has remained relatively stationary over the last 2000 years.See p7, lines 24-25.(3) That the spatial pattern of the sea ice area anomaly associated with a positive AMO index (as computed from the LMR) has remained relatively stationary over the last 2000 years.See p11, lines 3-4.(4) That the spatial pattern of the precipitation anomaly associated with a positive AMO index in the CCSM4 prior is similar to that reconstructed in the LMR over all time periods.See p11, lines 17-18.
Specific comments P1L20: Kushnir 1994 analysis was based on the hypotheses of the earlier Bjerknes 1964 paper We have now referenced both Kushnir (1994) and Bjerknes (1964).See p1, line 20.
P2L30: please cite also Murphy et al. 2017 andBellucci et al. 2017 which are more recent papers addressing the role of external radiative forcing in driving the AMV We have added references to Murphy et al (2017) and Bellucci et al (2017).See p3, lines 6-7.P2L10: "Several dynamical studies have suggested that zonal and meridional oscillations in the AMOC on multidecadal time scales may drive changes in north Atlantic SSTs (Dijkstra et al., 2006(Dijkstra et al., , 2008))."Can you elaborate what you mean by "zonal oscillations in the AMOC"?I am not aware of any zonal variations, usually the AMOC is seen as a meridional source of heat transport changes.
The description of the how zonal and meridional changes in the overturning drive a multidecadal oscillation in an ocean-only model may be found in Te Raa & Dijkstra (2002).We have included this reference and re-worded this sentence as follows (see p2, lines 13-16): Several dynamical studies, using both ocean-only and fully-coupled models, have suggested that zonal and meridional variations in the AMOC on multidecadal time scales may drive changes in north Atlantic SSTs (for a description of the dynamical mechanism, see Te Raa & Dijkstra 2002; see also Dijkstra et al 2006Dijkstra et al , 2008)).
P2L20: Tandon and Kushner 2015 also find that there's a lot of inter-model spread in the lead-lag correlation between the AMV and AMOC indices We agree.We have added this caveat to the description of the Tandon & Kushner ( 2015) study (see p2, lines 23-24).
P6L25: The fact that the anomalies are smaller after low-pass filtering is expected as a consequence of the low-pass filter.
We don't think that this must be the case.In general, correlations will be dominated by high-frequency variability.If the high-frequency variability is uncorrelated with an index, but low-frequency is strongly correlated, it is possible to obtain higher correlations following filtering.We agree that in this case, the anomalies are, indeed, smaller following the low-pass filtering.

Introduction
Modeling and observational studies have shown that north Atlantic sea surface temperatures (SSTs) co-vary with precipitation over the African Sahel (Zhang and Delworth, 2006), Eurasia (R. Lu and Ding, 2006;Zhang and Delworth, 2006), and North America (Enfield et al., 2001); hurricane development and intensity over the Atlantic (Zhang and Delworth, 2006); drought over the North American interior (McCabe et al., 2004;Nigam et al., 2011); summer temperatures over Europe and the Americas (Sutton and Hodson, 2005); sea ice thickness and extent over the Arctic (Miles et al., 2014); climate variability over the Pacific (Dong et al., 2006); and marine primary productivity (Henson et al., 2009).Improved prediction of global climate perturbations likely depends on better understanding of north Atlantic variability.Kushnir (1994), following earlier hypotheses of Bjerknes (1964), shows that north Atlantic sea surface temperatures (SSTs) appear to vary at both interannual and interdecadal time scales, and that variability at these different time scales is substantively different.While wind and surface pressure appear to covary with north Atlantic SSTs at short time scales, longer time scale variability appears to be uncorrelated with these fields.Kushnir (1994) concludes that these different time scales constitute 1 different phenomena, and suggests that the longer time scale portion of the variability could be driven by ocean-atmosphere coupling.This longer time scale variability is known today as Atlantic Multidecadal Variability (AMV), with the unforced, internally-driven component referred to as the Atlantic Multidecadal Oscillation (AMO).
Studies of the observational and paleoproxy records provide evidence of multidecadal variability centered over the Atlantic.
The Central England Temperature record, the longest observational temperature time series, suggests a 65-year time scale of variability over the last 350 years (Tung and Zhou, 2013).Evidence of multidecadal Atlantic variability (with time scales anywhere between 20 and 80 years) has also been found in proxy tree ring records (Delworth and Mann, 2000;Gray et al., 2004), annually-resolved ice cores (Chylek et al., 2011), and coral isotope records (Hetzinger et al., 2008).Lengthy proxy records extending over the last 8000 years also show multidecadal spectral power, though this power is not stationary over space or time (Knudsen et al., 2011).
Many studies have explored ocean-atmosphere interactions as driving factors for AMV, with changes in SSTs in the North Atlantic controlled by natural (unforced) fluctuations in the Atlantic meridional overturning circulation (AMOC) (see, e.g., Delworth et al., 1993;Polyakov et al., 2005;Zhang et al., 2007).Several dynamical studies, using both ocean-only and fullycoupled models, have suggested that zonal and meridional oscillations in the AMOC on multidecadal time scales may drive changes in north Atlantic SSTs (for a description of the dynamical mechanism, see Raa and Dijkstra, 2002;Dijkstra et al., 2006Dijkstra et al., , 2008)).While long-term observational records of the AMOC state are unavailable, observational evidence of sea surface height appear to support the idea that north Atlantic SSTs co-vary with changes in sea surface height along the eastern seaboard of the United States, which is consistent with changes in the AMOC (McCarthy et al., 2015).
Nevertheless, there remains disagreement regarding the role of AMOC changes in driving north Atlantic SST variability over multidecadal time scales.Tandon and Kushner (2015) show that the relationship between AMV and the AMOC may not be stationary between the preindustrial and industrial periods: while AMOC fluctuations appeared to lead changes in north Atlantic SSTs over the preindustrial period in models, this lead-lag relationship reversed during the industrial era with strong anthropogenic forcing of the climate system (though there was substantial variations in the lead-lag relationship between models).Recently, Clement et al. (2015) show that stochastic coupling between the atmosphere and oceanic mixed layer is sufficient to recover the AMV found in most climate models, calling into question the role of the oceanic MOC as a driver of the AMO.Indeed, only some state-of-the-art global climate models (GCMs) simulate multidecadal variability in Atlantic SSTs (Clement et al., 2015), while one older GCM has been shown to display spurious multidecadal variability due to poor parameterization of ocean mixing processes (see Danabasoglu, 2008).Hakkinen et al. (2011) suggest that changes in the strength of the north Atlantic subpolar and subtropical gyres may be more closely linked to the AMO than variations in the AMOC; such multidecadal shifts in the gyres alter atmospheric blocking patterns which may, in turn, perturb SSTs over the entire Atlantic basin.
The role of aerosols and other forcings as external drivers of north Atlantic SST variability has also been debated.Aerosol release by volcanic eruptions has been suggested to act as an external driver of AMV in the preindustrial period (Ottera et al., 2010;Knudsen et al., 2014), contrary to other studies that suggest that AMV is internally-driven by ocean-atmosphere interactions.Over the industrial period, the role of anthropogenic aerosols in driving AMV also remains an open question.Booth et al. (2012) argue that aerosol direct and indirect effects over the modern era (post-1850) have strongly impacted multidecadal SST variability over the north Atlantic, resulting in an evolution of north Atlantic SSTs that mimic an internallydriven oscillation.This argument, however, has been challenged by Zhang et al. (2013), who suggest that the GCM used in Booth et al. (2012) incorrectly modeled aerosol effects, and show that improved representation of these effects reveals that 20th century surface temperature variations can not be explained in full by changes in anthropogenic aerosol emissions.
Nevertheless, more recent work by Murphy et al. (2017) and Bellucci et al. (2017) further challenges the premise that north Atlantic SST variability over the industrial period is unaffected by anthropogenic forcings.
Given major disagreements between climate models in AMO representation, it's clear that study of AMO dynamics would benefit from using observational data from the real Earth system.However, detailed study of the AMO has been severely limited by the brevity of the instrumental record, which is less than two centuries long.Furthermore, these observations are from a time period characterized by intense anthropogenic forcing of the climate system through greenhouse gases, aerosols, land use changes, ozone-depleting substances, and others.While the instrumental record is short and likely biased towards a forced response, the proxy record does not suffer from this shortcoming.Indeed, annually-resolved proxies for temperature from tree rings, corals, ice cores, and others, exist at annual resolution over the last several millennia.Although proxy records exist over much longer periods than observations, an objective procedure for assimilating these proxies into coherent spatial fields remains an ongoing challenge.
The Last Millennial Reanalysis (LMR) provides a robust framework for such objective assimilation of proxies for the reconstruction of spatial fields (Hakim et al., 2016).In this study, we use climate field reconstructions made possible by the LMR to analyze the climate dynamics of the AMO.By assimilating temperature-dependent information from proxy records, we consider what thermodynamic features characterize the AMO over the last 2000 years, and whether these features have changed with time, particularly between the preindustrial (yrs 0 to 1850) and modern (yrs 1850 to 2000) eras.Since the data assimilation procedure used in the LMR relies on the use of a prior dataset (as described in Steiger et al., 2014;Hakim et al., 2016), we will also note what additional information assimilation of proxies yields beyond that found in this prior dataset.
Finally, we will use the LMR to consider what, if any, robust time scales characterize the AMO over the last two millennia, and how these time scales computed from assimilation of the observations compare to those found in previous studies.
The remainder of this paper is organized as follows.We describe the LMR proxy assimilation and climate field reconstruction procedure, along with details of our regression and time scale analysis methodologies in §2.In §3.1, we consider the AMO index over the last 2000 years, as inferred from the LMR, and in §3.2, we describe the basic climate state that corresponds to a positive AMO index.We highlight the energetics associated with a positive AMO index in §3.3, particularly top-of-atmosphere and surface fluxes, the implied meridional energy transports, and flux imbalances.In §3.4,we consider time scales of the AMO using wavelet analysis.We conclude by discussing our findings, their implications, and caveats of our approach in §4.
The climate field reconstruction methodology used in the LMR is described in detail by Steiger et al. (2014) and is validated for temperature reconstructions over the modern era (yrs 1850 to 2000 CE) in Hakim et al. (2016).LMR uses a novel ensemble data assimilation procedure, where information from a prior expectation of the climate, derived from a climate model, is weighted against information in proxy records.Weights are determined from the relative error in these two estimates of the climate, as defined by the update equation in the Kalman filter, which is optimal if the errors are Gaussian distributed.We proceed with an overview explanation of the Kalman update equation, before describing the solution method for this equation.
In brief, the data assimilation procedure updates the state of the prior, x p to some new state x a using information from the proxies y that is weighted using a Kalman gain matrix K: Here, H is the observation model that converts the prior state to its proxy equivalent, H(x p ).The innovation, y H(x p ), contains new information from the observations that is not already present in the prior.The Kalman gain matrix K, which weights the innovation and maps from proxy space to physical space, is where B is the error covariance matrix for the prior data, H is the linearization of the observation model H about the prior mean, and R is the error covariance matrix for the proxies (described further below).The numerator of matrix K, BH T , is the covariance expectation between the prior and the prior-estimated observations, and acts to spread information from proxy locales into the physical space.Solution of equations ( 1) and ( 2) is fully described in Steiger et al. (2014) and Hakim et al. (2016), and a summary of the essential elements follow.First, we describe our method for estimating the proxies from the prior data (i.e., H), and the associated errors (i.e., R).Second, we described our method for creating the prior and the associated errors (i.e., x p and B) using an ensemble sampling technique.This ensemble technique relates directly to the solution method for (1) and ( 2), which completes our summary.
The observation model H, a forward model which maps from physical space to the proxy space, is constructed using a calibration dataset and assumes a linear relationship between the proxy state and temperature over the calibration period: where T 0 is the annual-mean 2-m air temperature anomaly from the calibration dataset, 0 is the intercept, 1 is the slope, and ✏ is a Gaussian random variable with zero mean and a variance of 2 .The calibration temperature T 0 for a given proxy record is chosen from the gridpoint closest to the proxy location at concurrent time periods, and parameters 0 , 1 , and are computed using a linear least-squares fit against the NOAA Merged Land-Ocean Surface Temperature Analysis (MLOST; Smith et al., 2008) during the time period 1900-2000.The diagonal elements of R, the error covariance matrix for the proxies, are the variance of the regression residuals, 2 .In principle, other sources of error contribute to R, such as instrumental (laboratory) error, and representativeness error (the climate-model grid cell resolves spatial averages rather than the points that are observed), but the errors in (3) are large by comparison and so we approximate R accordingly.
In operational weather data assimilation, prior data is determined by a short-term forecast that is initialized from an analysis derived at an earlier time.Because the forecast applies to a short time interval (1-12 hours), it provides an accurate estimate of the observations at the future time; that is, it is a well-informed prior.In contrast, climate model forecasts on proxy timescales (yearly for the proxies used here) contain little forecast skill, and are expensive to compute.Consequently, model simulations on annual timescales are nearly agnostic of the initial state, and thus are statistically nearly the same as forecasts drawn randomly from the climate of the model.This fact motivates an "offline" approach to the prior and data assimilation method, where prior data is defined using an ensemble sampling strategy.Here we draw an ensemble of size 100 randomly from an existing climate model simulation.Specifically, we use the Community Climate System Model version 4 (CCSM4) Last Millennium run (from Phase 5 of the Climate Model Intercomparison Project) for the prior, and use the same 100-member sample from this simulation as the prior for each year in the reconstruction.The mean value of this sample serves as x p in equation ( 1).
We use an ensemble square-root solution method to equations ( 1) and ( 2), including serial observation processing (Whitaker and Hamill, 2002), so that proxies are assimilated one at a time for each year of the reconstruction.In this approach, the matrix inverse in equation ( 2) becomes a trivial scalar inverse, and the covariance matrix B, which can be enormous, is never explicitly formed; rather, only the ensemble-estimated covariance between the proxy and the reconstructed field at each point is needed (see Steiger et al. (2014) for further details).Proxy estimates H(x p are computed in advance, and updated each year along with the reconstructed fields with the assimilation of each proxy. Using this method, we reconstruct annual-mean fields between 1 to 2000 CE.The proxies used are a total of 465 timeseries from the PAGES2K dataset (Consortium, 2013;Hakim et al., 2016), which includes lake sediment calcium to phosphorus ratios, ice core water isotopes (both oxygen and hydrogen), coral oxygen isotopes, tree rings (width and density), and speleothem water isotopes.As described in Hakim et al. (2016), the LMR reconstructions of temperature and 500 hPa geopotential height have been validated over the instrumental period.Moreover, results have also been validated against independent proxies (not assimilated) before and during the instrumental period and to insure that the error in the ensemble mean is statistically equivalent to the predicted error in the ensemble variance as expected from theory.Though the reanalysis results obtained when using different prior and observation-model calibration datasets are quantitatively distinct (see Figure 12 in Hakim et al., 2016), they all yield qualitatively similar results.Therefore, much of the ensuing analysis, apart from calculation of time scales (see below), is performed using this single reanalysis, hereafter referred to as MLOST-CCSM4.
We note that because the prior is the same for every year, all temporal variability is determined by the proxies.As a result, the weighting in equation (1) involves a temporally invariant prior and temporally variable proxies so that, for any given proxy, reduced temporal variability may be expected.However, since the time series at any point in the reconstruction depends on the prior and many proxies, all having different errors, the variability at a proxy location may in fact be larger or smaller than that of a given proxy at that point.
Following Clement et al. (2015) and others, we quantify AMV using the AMO index, which is computed annually as the average of the area-weighted SSTs from 60N to the equator and 80W to the prime meridian.Links between climate variables and the AMO are defined by regression onto the AMO index and computing the value of the variable at two standard deviations of the AMO index.In order to isolate multidecadal (and longer) time scale variability, both the AMO index and variable fields are filtered using a 20-year low-pass Lanczos filter with 31 filter weights prior to computing the regressions (unless otherwise noted; see Duchon, 1979, for a full description of the Lanczos filter).The AMO index and climate variable fields are detrended before regression.Regression results from the LMR are computed for three time periods to highlight potential non-stationarity in the AMO: years 0 to 1850 (the preindustrial era), 1200 to 1850 (the part of the preindustrial era when proxy coverage is substantial), and 1850 to 2000 (the industrial era).These results are compared to similar results from the CCSM4 prior.
AMO index time scale analysis was performed with wavelets using the methodology described by Torrence and Compo (1998).The Paul wavelet is used with a scale spacing of 0.25 and a total of 52 scales spanning 2 to 1000 years.Significance is computed at p < 0.05, using a one-step autoregression (AR(1); i.e. red noise) process as the null hypothesis; significance 10 results were found to be insensitive to choice of the mother wavelet.Wavelet analysis was performed on six different AMO index reconstructions using 2 different priors and 3 different calibration datasets (see Table 1).

The AMO index in the LMR
Over the last 2000 years, the AMO index reconstructed from the LMR resembles the global mean surface temperature (GMST) timeseries (Figures 1a and b, respectively; r = 0.97 at lag 0).Like the GMST, the AMO has the least variability over the first 1000 years of the reanalysis, possibly due to the lack of globally-distributed proxy records over this early period of the reanalysis.Over the latter portion of the reanalysis (year 900 CE forward), a cooling trend is evident in both the AMO and GMST, which is punctuated by a warmer period centered at year 1500 CE.A warming trend is evident in both the AMO and GMST timeseries following year 1900 CE.
In Figure 1c, we show the LMR reconstruction of the AMO index over years 1900 to 2000 CE, along with three different reanalysis datasets: NASA's GIStemp, ERA-20C, and the 20CR.Comparison with these other reanalyses suggests that the LMR MLOST-CCSM4 reconstruction reasonably re-creates the the north Atlantic SST record over the last century: as expected, the LMR and the three reanalyses all show an upward trend in the AMO index over this time period, and all reveal intermittent periods of rapid and slower temperature rise (r = 0.64, 0.67, 0.73 between the LMR MLOST-CCSM4 AMO index reconstruction and the GIStemp, ERA-20C, and 20CR AMO indices, respectively), suggesting a multi-decadal oscillation.The LMR displays slightly weaker amplitude in this multi-decadal timescale than the other reanalyses.We return to this question of timescales in §3.4.

Thermodynamics of the AMO in the LMR
Globally, the temperature field associated with a positive AMO index is characterized by warming over the continents and especially strong warming over the far northern reaches of the Atlantic and Arctic (Figure 2).Warming is greater over the northern hemisphere (NH) than the southern hemisphere (SH) over all time periods, and the magnitude of warming is strongest over all regions in the CCSM4 prior compared to that in the LMR.Over both the CCSM4 prior and the LMR, there is warming over the Arctic, though the location of maximum warming is different: in the CCSM4 prior, warming is greatest over the Greenland-Iceland-Norwegian (GIN) seas, while warming is greatest in the LMR over the Barents sea.In the LMR, variability is relatively stationary over all time periods (r = 0.92, 0.73 between the regression over years 0 to 1850, and years 1200 to 1850 and years 1850 to 2000, respectively), though differences in warming over the high northern latitudes are apparent; in the pre-industrial period (pre-1850), warming over northern Eurasia is most prominent, while warming over boreal North America and central Asia dominates in the industrial period.
Comparison of the 20-year low-pass filtered regression of the AMO index on temperature (Figure 2a to d) with the unfiltered regression (Figure 2e to h) reveals that the temperature response associated with the AMO is stronger at the shorter time scales that dominate the unfiltered regression.
Over all time periods in the LMR (years 0 to 1850, years 1200 to 1850, and years 1850 to 2000), warming is not uniform over the north Atlantic basin: maxima in warming are evident at 50N and 5N (Figure 3a to d).Over shorter timescales than those explored in the remainder of this study, however, the two-lobed pattern becomes a horseshoe with the tropical and midlatitude

9
warming maxima linked by warming over the eastern Atlantic (compare Figure 3a-d to Figure 3e-h).This distinct horseshoe pattern of warming over the midlatitudes and tropics that characterizes north Atlantic SSTs variability at shorter time scales has been noted by others, and is the leading mode of variability found with empirical orthogonal functional decomposition analysis over this region (see, e.g., Hartmann, 1994).Because our focus is on longer timescale variability, we will utilize 20-year low-pass filtering for the rest of our analysis (except for time-series analysis of the AMO index, §3.4).hemisphere or a decrease in reflected SW due to fewer clouds in the warmer hemisphere.The outgoing LW flux increases nearly globally, commensurate with increased temperatures at the effective radiating level associated with warming during the positive phase of the AMO.This increase in outgoing LW radiation dominates in the SH, while the decrease in outgoing SW radiation dominates in the NH, implying a net southward meridional energy transport in the LMR (see Figure 8a and accompanying text).At the equator, there is an increase in the outgoing SW radiation, which is balanced by a decrease in the outgoing LW radiation; both are consistent with greater cloudiness in the deep tropics and strengthening of the equatorial precipitation noted earlier (Figure 5).
In the CCSM4 prior, the TOA radiative fluxes suggest a very different picture.While the outgoing SW flux decreases (Figure 6c) and the outgoing LW flux increases (Figure 6d) as in the LMR, the sum of these, the net upward TOA flux (Figure 6a), is positive in the high latitudes in both hemispheres and negative in the lower and mid-latitudes; as a result, the implied meridional total energy transport in the CCSM4 prior is polewards in both hemisphere.
At the surface, both the LMR and the CCSM4 prior have net fluxes associated with a positive AMO index that are downwards in the tropics, upwards in the subtropics and midlatitudes, and downwards in the NH high latitudes (Figure 7a).This pattern is most pronounced in the CCSM4 prior, but is also evident over all time periods in the LMR.Globally, we find that this spatial pattern results from decreases in both the net (upward) surface SW radiation (i.e. an increase in SW coming down at the surface; Figure 7b) and the net (upward) surface LW radiation (i.e. an increase in the LW coming down at the surface; Figure 7c); these increased surface SW and LW radiative fluxes are mostly balanced by an increase in the surface latent heat flux (Figure 7e), which renders the net surface flux small over most latitudes.At the equator, an increase in the net (upward) surface SW radiation and a decrease in the net (upward) surface LW radiation is consistent with greater cloud cover.Sensible heat flux anomalies at the surface are small everywhere (Figure 7d).These results from the LMR suggest that the positive phase of the AMO is associated with greater atmospheric water vapor, likely a result of Clausius-Clapeyron scaling of atmospheric moisture with temperature.These findings, together with the TOA fluxes described earlier, also suggest that the positive phase of the AMO is associated with fewer (reflective) low clouds in the extra-tropics and more clouds near the equator; the latter is in agreement with positive precipitation anomalies in the deep tropics (recall Figure 5).
As alluded to earlier, differences in the net TOA radiative fluxes between the LMR and the CCSM4 prior imply very different total energy transport anomalies.The total implied meridional energy transport, T ET , at latitude 0 can be computed using the zonally-averaged net TOA fluxes R T OA ( ) as where and ✓ are the latitude and longitude coordinates, respectively, r earth is the radius of Earth, and R storage ( ) is the heat storage tendency of the climate system at latitude (Peixoto and Oort, 1992).We find that the T ET is anomalously southward north of 20S and northward south of 20S for the LMR over later time periods (1200 to 1850 and 1850 to 2000), and is particularly large between 10S and 20N (Figure 8a).In the CCSM4 prior, on the other hand, the T ET is anomalously southward south of 10N and anomalously northward north of 10N.
The atmospheric energy transport, AET , at latitude 0 is computed using the difference between R T OA ( ) and net surface flux R Sf c ( ) as The AET anomaly associated with a positive AMO is southwards at most latitudes in both the LMR (over all time periods) and in the CCSM4 prior (Figure 8b).The maximum southward AET is further northward in the LMR compared to the CCSM4 prior (close to the equator in the former, while in the southern midlatitudes in the latter).The increase in southward crossequatorial energy transport in the LMR is consistent with the northward shift of the ITCZ noted previously (recall Figure 5).The increase in southward cross-equatorial energy transport is also consistent with an increase in the inter-hemispheric temperature difference associated with the AMO which preferentially warms the NH over the SH (recall Figure 2; also see Friedman et al., 2013), which arises from a strengthening of the Hadley circulation in the cooler hemisphere.On the other hand, the increase seen in the CCSM4 prior occurs most prominently in the midlatitudes, suggesting a strengthening of poleward energy transport by midlatitude eddies in the SH rather than an adjustment of the tropical circulation.
The oceanic energy transport OET ( 0 ) is computed as a residual from the T ET and AET : From Figure 8c, it is clear that, aside from the subtropical SH, the LMR and CCSM4 prior differ significantly in anomalous OET associated with the positive phase of the AMO.In the LMR (over 1200 to 1850 and 1850 to 2000), energy transport by the ocean is anomalously southward in the NH and northward in the SH, with a significant southward cross-equatorial component; much of the anomalous transport is confined to the subtropics.These anomalies are larger over the 1850-2000 time period, and are consistent with a decrease in oceanic energy transport by the wind-driven subtropical gyres and subtropical cells when the AMO is in its positive phase.A large southward cross-equatorial component to the OET anomaly suggests a decrease in northward cross-equatorial energy transport by the AMOC.In the CCSM4 prior, on the other hand, the OET response is anomalously positive at all latitudes, suggesting that an increase in northward energy transport by the AMOC corresponds to a positive AMO index.
In both the LMR and CCSM4 prior, we find that the average net TOA flux over the last 2000 years is positive, indicating that energy is being removed from the Earth system; when the AMO index is positive, however, the net TOA flux decreases such that energy is being anomalously added to the Earth system relative to the mean state (Table 2).The extent to which this occurs is similar between the CCSM4 prior and the LMR over all time periods, suggesting that this aspect of the AMO has been stationary over the last 2000 years.
Net TOA Imbalance (W/m 2 ) Global Mean Regressed on the AMO Index  2. Net top-of-atmosphere (TOA) flux imbalances in the LMR (from years 0-1850, 1200-1850, and 1850-2000) and the CCSM4 prior, shown for the global mean and regressed on the AMO index.Fluxes out of the Earth system (i.e.upward) are positive.
Surface flux imbalances in the LMR (Table 3) are also qualitatively similar to the TOA flux imbalances, with positive surface flux anomalies characterizing the mean state over all time periods (i.e.energy is leaving the oceans), consistent with millennialscale cooling.When the AMO is in its positive phase, on the other hand, the surface flux anomaly decreases, indicating that energy is being anomalously absorbed by the oceans relative to the mean state.This is consistent with the TOA flux imbalance linked to a positive AMO index, which also decreases relative to the mean.The results from the CCSM4, however, are more energetically consistent in the magnitude of the energy uptake by the climate system when the AMO is in its positive phase (0.19 W/m 2 and 0.18 W/m 2 for the TOA and surface flux imbalances, respectively) compared to the LMR (e.g.0.01 W/m 2 and 0.13 W/m 2 for the TOA and surface flux imbalances, respectively, in the LMR over the years 0 to 1850 CE).-1850, 1200-1850, and 1850-2000) and the CCSM4 prior, shown for the global mean and regressed on the AMO index.Fluxes out of the surface (i.e.upward) are positive.

Time Scales of Variability
We now consider what time scales characterize the AMO index.Tung and Zhou (2013) show a distinct 50-to 80-year peak in the Central England Temperature record (HadCET; see Parker et al., 1992), and suggest that this variability may be related to coherent variability in north Atlantic SSTs.Other studies have synthesized proxy records from the north Atlantic basin over the last 500 years to infer the existence of such a multidecadal oscillation over the last millennium (see, i.e.Delworth and Mann, 2000;Gray et al., 2004).Furthermore, Knudsen et al. (2011) suggest that multidecadal variability, non-stationary and possibly linked to north Atlantic SSTs, is evident in global ice core and sedimentary records over a span of the last 8,000 years.

Motivated by previous research, we now analyze what time scales are evident in the LMR reconstruction of north Atlantic
SSTs.The LMR uses records from proxies over the last two millennia, weighted heterogeneously relative to the prior, and includes the proxy data used by Delworth and Mann (2000), Gray et al. (2004), and others; these records are processed here using the data assimilation algorithm described herein (see Methods, §2; also Hakim et al., 2016).If time scales present in individual proxy record time series are indicative of variability over the large-scale north Atlantic in general, these time scales should also be present in the reconstructed AMO index.
Wavelet analysis of the AMO index is performed on six different LMR reconstructions that use three distinct calibration datasets (MLOST, GIS, and CRU; see Table 1) and two different model priors (CCSM4 and MPI; see Table 1).The reconstructions reveal statistically significant power over different time periods at several different time scales, and the details vary with the reconstruction (Figure 9).For all reconstructions, longer time scale variability (> 40 years) is only evident over the  We have considered the thermodynamics of the AMO over the last 2000 years using the LMR, a paleoclimate reanalysis method that objectively assimilates information from paleoclimate records.Several aspects of our findings regarding the AMO agree with previous model-based and observational studies of north Atlantic SST evolution.We find that a positive AMO index coincides with warmer continents, a two-lobed pattern of warming over the Atlantic, and a warmer Arctic (Kushnir, 1994;Delworth and Mann, 2000;Chylek et al., 2009).In the Arctic, sea ice retreats from the Greenland-Iceland-Nordic seas and thins over much of the Arctic Ocean (Miles et al., 2014).Globally, precipitation increases and the ITCZ strengthens and shifts northward (Knight et al., 2006).
With the LMR, we also find several elements of the AMO that have not been reported in the literature.When the AMO is in its positive phase, southward cross-equatorial energy transport increases, mostly mediated by the atmosphere (the northward shift in the equatorial precipitation); these latter changes are consistent with energy balance requirements as necessitated by stronger warming in the NH than the SH during the positive phase of the AMO, i.e. an effect of the increase in the interhemispheric temperature anomaly (see Friedman et al., 2013).The LMR also shows that when the AMO is in its positive phase, the Earth system loses less (net) energy to space, with much of this excess energy absorbed by the oceans.TOA and surface flux anomalies are consistent with an increase in atmospheric specific humidity and a decrease in low (reflective) clouds.
Our study of the AMO using the LMR demonstrates that climate field reconstruction methods, which assimilate information from the proxy record, can provide valuable information on climate variability in the Earth system.While correlations between north Atlantic SSTs and various thermodynamic fields (temperature, precipitation, and sea ice) are similar in the LMR and the CCSM4 prior, the LMR provides a temporal reconstruction of the AMO index that informs our understanding of north Atlantic SST variability beyond that from the model prior.First, energetic changes associated with the AMO are distinct in the LMR and CCSM4 prior, with the LMR predicting consistent changes in the cross-equatorial energy transport, tropical circulation, and extratropical cloud cover that are not found in the CCSM4 prior.Second, the LMR helps resolve the dominant time scales that characterize the AMO.Since there is little agreement between various GCMs regarding the dominant time scales that characterize the AMO (see Clement et al., 2015), LMR reconstruction of the AMO index is invaluable.Our results suggest that the proxy observations over the last 2000 years, when objectively assimilated, do not exhibit a multidecadal timescale.
The lack of a distinct multidecadal spectral peak in the LMR reconstruction of the AMO is in contrast to other studies that have found such variability in individual observational records (see, e.g.Tung and Zhou, 2013) or limited collections of proxies (see, e.g.Delworth and Mann, 2000;Gray et al., 2004).We point out that such a significant spectral peak in an individual record or collection of records does not necessarily translate into a coherent mode of basin-scale multidecadal variability.While certain records may display oscillations, a basin-scale oscillation requires both spatial coherence and matching time scales in these records over certain regions.The objective assimilation procedure used in the LMR climate field reconstruction utilizes the information provided by the proxy records investigated in these previous studies; the results of this objective assimilation suggest that there is no distinct multidecadal or centennial spectral peak in the AMO index, though there is reddening of the spectra.These results support the null hypothesis presented by Clement et al. (2015), and suggest that there may be little need to consider longer time scale processes when studying the mechanism of AMO over the late Holocene, particularly the preindustrial period.
In our reconstructions, multidecadal spectral power is evident only over the latter portion of the LMR, year 1500 CE forward, and is particularly pronounced following year 1900 CE.While there has been much hypothesized regarding these oscillations in the instrumental record, and whether they are a result of internal climate variability or a result of external forcing (see, e.g.Booth et al., 2012;Zhang et al., 2013;Murphy et al., 2017;Bellucci et al., 2017, and others), we point out two important characteristics of the instrumental period that render it a poor era for studying multidecadal variability in the climate system.
First, the instrumental record is very short; as a result, there is very little that can be said about multidecadal time scale variability over this time period that carries any statistical weight (see, e.g., Wunsch, 1999;Vincze and Janosi, 2011).Second, the instrumental period is one in which the climate system has been strongly anthropogenically forced, suggesting that any variability observed over this period may differ significantly from variability from the preindustrial era when anthropogenic forcing was much smaller.Indeed, Tandon and Kushner (2015) show that in models, the lead-lag relationship between north Atlantic SSTs and the AMOC are very different between the preindustrial and modern periods, suggesting a shift in the mechanism underlying AMV between these two time periods.Both of these limitations suggest that caution must be exercised in extrapolating characteristics of multidecadal variability in the Earth system from the instrumental record.
We also point out that our objective data assimilation approach cannot completely rule out a distinct multidecadal spectral peak in north Atlantic SSTs.While we have shown that objective data assimilation of the proxy record to reconstruct the surface temperature field over the last two millennia does not yield any evidence of multidecadal variability, it is likely that further improvements in the proxy network, the proxy system models, or the data assimilation procedure will improve the reconstruction in such a way as to modify our current conclusions.From this study, however, we can say that the surface temperatures reconstructed over the last 2000 years, which are obtained from assimilating the existing network of proxies using today's state-of-the-art methods, do not provide compelling evidence for multidecadal oscillations over the north Atlantic.
We conclude by pointing out some important limitations of studies of climate variability using the LMR.First, the data assimilation approach itself is Gaussian, which may limit its utility in regimes very different from that of the base climate state.
Second, LMR reconstructions depend on the climate models, proxies, and observation models that are used to create them.These components, in turn, are being continually tuned and refined.Third, it is unknown what effects the limitations in the size and distribution of the proxy network have on the reconstruction made possible by the LMR.In particular, we note that the more limited network from the early portion of the instrumental record (pre-1200 CE) may affect the reconstruction and its variability.Further study will be required to assess sensitivity to the sparsity of the proxy network.Finally, we note that there are other ways to improve paleoclimate reanalyses using proxy records, including the incorporation of linear inverse modeling into the reconstruction methodology to achieve on-line data assimilation (see Perkins and Hakim, 2016).In the future, such refinement of reconstruction methods will help to further study of the Earth system over its long and varied climatic history.

Figure 1 .
Figure 1.(a) Atlantic Multidecadal Oscillation (AMO) index time series and (b) global mean surface temperature (GMST) time series computed from the LMR over the last 2000 yrs, with the yearly time series (light blue, dotted) and 20-year low-pass time series (black; 95% confidence interval shown in orange); and (c) AMO index over the last 100 yrs (bottom) as computed from the LMR (black; 95% confidence interval shown in grey), the 20th Century Reanalysis (20CR, red), the ERA-20C reanalysis (Era-20C, blue), and the NASA GIStemp reanalysis (GIStemp).

Figure 2 .
Figure 2. Low-pass filtered regression of the AMO index on surface temperature (in K) in (a) the CCSM4 prior; (b) the LMR from years 0 to 1850; (c) the LMR from years 1200 to 1850; and (d) the LMR from years 1850 to 2000.Panels (e) through (h) are similar to (a) through (d), but show the unfiltered regression.

5Figure 3 .
Figure 3.As in Figure 2, over the north Atlantic basin.

Figure 4 .
Figure 4.As for Figure 2, panels (a) to (d), but for the low-pass filtered regression of the AMO index on sea ice thickness (colors, in m) and on sea ice concentration (contours, %).

Figure 5 .
Figure 5.As for Figure 2, panels (a) to (d), but for the low-pass filtered regression of the AMO index on precipitation (in mm/day).

Figure 7 .
Figure 7. Low-pass filtered regression of the AMO index on surface energy fluxes (in W/m 2 ) in the LMR (from years 0 to 1850, 1200 to 1850, and 1850 to 2000) and the CCSM4 prior: (a) net (downward) surface flux, (b) net (upward) SW flux, (c) net (upward) LW flux, (d) sensible heat flux, and (e) latent heat flux.

Figure 8 .
Figure 8. Low-pass filtered regression of the AMO index on energy transport (in PW) in the LMR: (a) total energy transport, (b) atmospheric energy transport, and (c) oceanic energy transport for the LMR from years 0 to 1850 (purple lines), LMR from years 1200 to 1850 (blue lines), LMR from years 1850 to 2000 (green lines), and CCSM4 prior (red).

Figure 10 .
Figure 10.Global wavelet spectra of reconstructions shown in Figure 9 (a) for time scales between 0 to 200 years, and (b) zoomed in to show time scales between 2 and 10 years.Time scales that are statistically significant over the red noise background (at p < 0.05) are shown as solid lines, while time scales that are not statistically significant are shown as dotted lines.In units of variance of the AMO index, K 2 .

Table 1 .
Prior and calibration datasets used in the LMR reconstruction.The default datasets are used for all analyses, while additional datasets are only used in the wavelet time scale analyses.

Table 3 .
Net surface flux imbalances in the LMR (from years 0