A 1150-year-long AMV reconstruction suggests early warning for a North Atlantic climate tipping point

The Atlantic Multidecadal Variability (AMV) is a large-scale climate phenomenon with 18 crucial impacts on human societies and ecosystems. Its periodicity and drivers are 19 controversial due to the short temporal extent of instrumental observations and competing impacts of external forcing and internal variability. Here, we use a well-verified set of 21 paleoclimate proxy records and compare four regression methods to perform different 22 reconstructions of the AMV since 850 C.E., built to only reflect North Atlantic internal 23 variability. The best performing reconstruction, when verified both against climate model outputs and independent proxy records is obtained using the non-linear random forest 25 method. It exhibits large multi-decadal variations in the range of 20-90 years, a broader 26 range than the 50-70 years identified in instrumental records. The reconstruction shows 27 that AMV autocorrelation properties have experienced significant changes in the recent 28 decades, suggesting an early warning signal for the proximity of a tipping point in the 29 Atlantic.

have been further supported by a tide gauge-based reconstruction of the ocean circulation 48 intensity in the North Atlantic intergyre region that exhibits strong correlations with the AMV 49 index over the last 60 years 17 . Several long paleoclimate records also indicate the existence 50 of such multidecadal variability on longer time scales 18,19 . Nevertheless, it has also been 51 shown that many features of the observed AMV can be reproduced in slab-ocean models by 52 the integration in the ocean of the stochastic forcing of the North Atlantic Oscillation, thus 53 without any influence of ocean dynamics 20 . It has also been highlighted that periods of cold  The AMV index definition is of paramount importance to tackle correctly the lingering 85 questions concerning its drivers. Since climate is currently in a period with a dominant 86 anthropogenically-induced warming signal, several methods are used to isolate the internal 87 variability inherent to the North Atlantic basin. In this study, we use three yearly AMV indices 88 in which all the externally forced signals (including anthropogenic forcing, solar variations, 89 and volcanic eruptions) in annual NASST have been removed with different techniques, to 90 thus explore the sensitivity to the method. We denote these indices as 8 , 2 and with the target AMV index over the instrumental period. They also do not include any pseudo-120 proxy validation of the reconstruction as a means to verify that the assumed climatic Results: 129 An AMV reconstruction based on optimal regression model approach 130 The availability of large numbers of annually-resolved proxy records sensitive to temperature 131 and precipitation provides a key opportunity to test and develop new objective 132 reconstruction methodologies using advanced techniques. The different reconstructions 133 produced in this study use proxy records significantly correlated at the 95% confidence level 134 with the target AMV indices (cf. Fig. 1). The selection of these proxies is made from a large set 135 of Northern Hemisphere and annually-resolved proxy records that includes the PAGES 2k  numbers of proxy records (Fig. 2b).

185
The nested reconstruction uses a total of 55 Northern Hemisphere proxy records.

186
Their a posteriori weights, given by random forest importance (cf. Methods), and their 187 temporal availability, are presented in Fig. 2c instrumental summer and annual temperatures in this region, although it is observed for 205 some precipitation time series over the historical period ( Fig. 1). Finally, the fact that only five 206 proxy records from Europe are used is mostly due to their relatively reduced presence in the 207 proxy record database we use (<10%). However, we find that one of them has a relatively 208 large weight for the reconstruction (>7%) and covers the entire reconstruction period. It 209 corresponds to a time series of tree rings growth measurements from European Alps 44 , which 210 is strongly correlated with summer temperature over the historical period (r=0.7, p<0.01,

211
Extended Data Table 3). The four other European proxy records are related to either summer 212 or annual temperature and precipitation. The selection of these records for our 213 reconstruction is thus highly consistent with the well-documented fingerprint of the AMV on proxy record used for our reconstruction and their correlations with instrumental data and 216 AMV is given in Extended Data Table 3. Additionally, the same maps of the AMV fingerprints   Table 4). While 308 the reconstructed NASST has a similar response a decade after the eruption than a previous reconstruction also based on NASST 34 , no significant response is found for the AMV 310 reconstruction (Fig. 6). This is partly expected given that we have reconstructed an index of The wavelet analysis in Fig. 7 shows that the AMV reconstruction exhibits important 330 multidecadal variations (Fig. 7), contrary to what has been recently suggested for an                     weakening of the subpolar gyre as trigger of LIA-type episodes. Clim. Dynam. 48(3-4),                   averaged SST anomalies as a proxy for the externally forced signal that is subtracted from 630 NASST, resulting in the index 8 (Fig. 1a). This approach has limitations since it does not 631 account for regional variations in external forcings, such as the distribution of anthropogenic 632 aerosols, which was denser in North America-Europe in the middle of last century and became 633 more prominent in Asia in recent decades 60 . To circumvent this problem, a different index, 634 the (Fig. 1a), is built using climate model historical simulations to isolate the forced The first analyses when preparing this study singled out a proxy record from Asia (named 661 "Asia.MOR1JU" in P2k-ALL), which had abnormally large RF weights (more than 5 times higher 662 than the second) as compared to all the other proxy records. To prevent highly biasing the 663 reconstruction towards this single proxy we have decided to remove it from the database 664 used in this study.
The reconstruction procedure of this study is then based on the resulting database (P2k+). In 666 the next method sections, we describe how different reconstructions are compared and how 667 a final nested reconstruction of the AMV is obtained. These reconstructions will also use 668 correlation tests to select proxies from P2k+ that are significantly correlated with a given AMV 669 index, for a given learning and reconstruction period. This means that only a subset of the 670 most relevant proxy records from the P2k+ database is finally used in each reconstruction.  Computation of a reconstruction and evaluation: 686 We define the reconstruction period as , defined by annual time steps, and the common 687 period of the proxy records and the AMV index as , in this case defined by < annual 688 time steps such that ⊂ . We then define the AMV index as ∈ and the matrix of the available proxy records as ∈ × . We finally denote as ∈ × the sub-matrix of that Where is the length of the testing sample. The weights of the proxy records used for the nested reconstruction are presented in Fig. 2c. 745 Those weights have been calculated using the random forest variable importance 40 . Different importance over the 1001 reconstructions. The same is done in Fig. 3b, but each importance 756 is also averaged over the PPE performed on the 12 members of CESM1-LME used.

758
To reconstruct the NASST, the same procedure as for the AMV is applied to select the optimal 759 reconstruction approach, with the only difference that only one index definition is considered, 760 as no method to remove the forced variability is applied. This means that the optimal model  There is an exception for some ocean records (<10 from Fig 3a)  The "Proxy records database" method section emphasizes the fact that the proxy records do   Early warning signal test: 822 We base our approach on methods for the detection of incoming climate tipping points 50  , under the null hypothesis 856 0 : " = 0"which is tested against the alternative hypothesis 1 : " ≠ 0". The p-857 value (shown in Fig 4b) of the test is deduced from the quantile of this distribution.   Comparison with independent ocean proxy records. a) 37 Ocean 2k proxy records35 signi cantly correlated at least at the 90% con dence with the AMV over the pre-industrial period (i.e., prior to 1870). b) 10-years (black) and 30-years (blue) kernel smooth of the AMV reconstruction.   Comparison of the AMV with NASST and volcanic forcing. a) Final reconstructions of AMV and NASST. b) Superposed epoch analysis47 for responses of the AMV and NASST reconstructions to the ten largest eruptions46 of the last millennium (see Extended Data Table 4). Composite series are performed for 31 years, for which the 11th are the actual years of the eruptions. Each individual response is centered to its values 10 years before the eruption (from N-10 to N-1) before computing the composite time series. 95% con dence levels have been calculated using a Monte-Carlo approach47.

Figure 7
Spectral analysis and early warning signal. a) Discrete wavelet transform of the nested AMV reconstruction from this study. Contours provide a 90% con dence level. The white line and the light white-shaded area below indicate the cone of in uence. The cone of in uence gives the spectrum borders where the edge effect (i.e., the time boundary effect) becomes too important, which cannot be robustly interpreted. b) Early warning signal test50-52 of the nested AMV_F reconstruction (Methods) based on AR1 coe cients, for different window lengths (WL). For each WL, sliding AR1 coe cient are computed and a Kendall τ statistics between time and the sliding AR1 time series are calculated. Signi cance is approximated using Gaussian distributions because of the large length (>50) of the AR1 coe cients (see Methods).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.