Is natural variability really natural? The case of Atlantic Multidecadal Oscillation investigated by a neural network model

Is Atlantic Multidecadal Oscillation a genuine representation of natural variability in the climate system? Or perhaps is it strongly forced by external drivers? In this paper, a data-driven attribution investigation has been performed for the Atlantic Multidecadal Oscillation (AMO) behaviour in the past via a machine learning technique, NN modelling. We clearly see a forced nature of AMO in the last 150 years, with a strong contribution of the forcing coming from anthropogenic sulphates, which induces its typical oscillating behaviour. The following original application of our model to future predictions of the AMO behaviour shows that it shall probably lose its oscillating characteristic features. The only way to recover them is to consider an unrealistic increase in anthropogenic sulphates in the future under a strong mitigation scenario, and possibly a low-power solar regime. Due to the established influence of AMO on climate and meteorological phenomena in several regions of the world, our results can be important to better understand the past and envisage several future scenarios.


Introduction
As well known, the climate system is endowed with a complex dynamics, characterised by both internal natural variability and interactions with external (natural and anthropogenic) forcings. Attribution studies of the recent global warming show that the global temperature (T) has been mainly driven by external anthropogenic forcings, even if internal variability can have a role in determining temperature values at short time scales (interannual to multidecadal ones).
If El Niño Southern Oscillation (ENSO) clearly contributes to interannual variability, not only in global temperatures (Timmermann et al. 2018), decadal 'modulations' of the global T curve have not been strictly attributed to atmospheric or oceanic (possibly coupled) patterns. Nevertheless, an unforced internal component that varies on multidecadal time scales has been identified in the temperature record (DelSole et al. 2011) and it has been noted that the warming and cooling of this component matches that of the Atlantic Multidecadal Oscillation (AMO).
It has been shown (see, for instance, Knight et al. 2006) that AMO influences climate and meteorological phenomena in several regions of the world, e.g. Sahelian rainfall, Atlantic hurricanes and European temperatures in summer. But, could AMO-as an internal natural oscillation in atmosphere-ocean-land-cryosphere interactions-have a role in modulating global T? Or, very differently, are there some external forcings which drive the AMO behaviour? This dilemma can be hopefully solved by attribution studies of the AMO behaviour itself.
As a matter of fact, many studies faced this problem (Knight et al., 2005;Jungclaus et al. 2005;Otterå et al. 2010;Booth et al. 2012;Zhang et al. 2013;Knudsen et al. 2014;Clement et al. 2015;Bellucci et al. 2017;Cane et al. 2017;Murphy et al. 2017;Bellomo et al. 2018;Kim et al. 2018;O'Reilly et al. 2019;Athanasiadis et al. 2020). Even if several papers emphasised the role that external forcings (sulphate aerosols in particular) could have in driving the AMO curve, at least in the last century, others claimed to a role of ocean dynamics and internal factors of the climate system.
The first explanations for the AMO behaviour have focused on the role of naturally occurring changes in ocean circulation, primarily the Atlantic Meridional Overturning Circulation (AMOC), basing on climate model simulations 1 3 with constant external forcing that exhibit multidecadal climate variability with pattern and amplitude that resemble the observed AMO: see, for instance, Knight et al. (2005) and Jungclaus et al. (2005).
Then, Otterå et al. (2010) found that the phasing of the multidecadal fluctuations in the North Atlantic during the past 600 years is, to a large degree, governed by changes in the external solar and volcanic forcings: volcanoes play a particularly important role in their results.
In the study by Booth et al. (2012), anthropogenic aerosol emissions and periods of volcanic activity explained 76% of the simulated multidecadal variance in detrended 1860-2005 North Atlantic sea surface temperatures. Aerosols were recognized as a prime driver of twentieth-century North Atlantic climate variability. Zhang et al. (2013) found discrepancies which cast doubts on the claim by Booth et al. (2012), but no alternative explanation for the AMO behaviour has been supplied. Knudsen et al. (2014) found that external forcing played a dominant role in pacing the AMO after termination of the Little Ice Age. Furthermore, this paper (based on statistical analyses of high-resolution proxy data) suggests that the AMOC is important for linking external forcing with North Atlantic sea surface temperatures, a conjecture that reconciles the previous opposing theories concerning the origin of the AMO.
The results by Clement et al. (2015) show that AMO is not driven by ocean circulation and AMOC, but is forced by mid-latitude atmospheric circulation.
Basing on a multimodel analysis, the study by Bellucci et al. (2017) suggests that anthropogenic aerosols and greenhouse gases (GHGs) might have played a key role in the 1940-1975 North Atlantic cooling. Cane et al. (2017) performed a modelling study which concludes that the ocean is influent on AMO, but not that it is important to the simulation of the climate variability as represented by the surface temperature. Murphy et al. (2017) analysed the AMO in the preindustrial and historical simulations from the Coupled Model Intercomparison Project Phase 5 (CMIP5) to assess the drivers of the observed AMO from 1865 to 2005. They concluded that there is an essential role for external forcing in driving the observed AMO and that it is very unlikely that internal variability is sufficient to drive multidecadal changes resembling the AMO observed over the last century. Bellomo et al. (2018) found that phase changes of the AMO over the years 1854-2005 can be explained only in the presence of varying historical forcing. Furthermore, they found that internal variability is large in North Atlantic SST at timescales shorter than 10-25 years, but at longer timescales the forced response dominates. Moreover, single forcing experiments show that GHGs and tropospheric aerosols are the main drivers of the AMO in the latter part of the twentieth century. Kim et al. (2018), from observations and model simulations, support a key role for ocean dynamics, rather than external forcings, in Atlantic multidecadal variability during the last half century. O'Reilly et al. (2019) used proxy data over a period of more than 300 years to investigate the internally generated AMO and the externally forced AMO components. In addition, they analysed a long ensemble simulation over the same period to estimate the external forcing of the AMO in the proxy datasets. The proxy AMO is found to closely follow the accumulated forcing of the North Atlantic oscillation (NAO) over almost the entire analysis period, referred to as an 'internal' source of AMO.
Finally, in a study about predictability, Athanasiadis et al. (2020) provided evidence that ocean dynamics may drive part of the observed decadal atmospheric variability. In their vision, the NAO partly drives ocean circulation anomalies and AMO.
Obviously, one can argue that the signal of AMO has been estimated for centuries and millennia well before the industrial era: see, for instance, Knudsen et al. (2011). Thus, the studies cited above seem quite limited for establishing the unforced or forced nature of AMO. Recently, however, Mann et al. (2021) addressed this problem and showed this oscillation can be due to pulses in volcanic activity in the last millennium. The great majority of these studies has been performed by runs of Global Climate Models (GCMs).
In this framework, following a first attempt briefly described in Pasini et al. (2017), here we address this problem via a neural network (NN) model specifically developed for modelling relationships amongst variables in small datasets (Pasini 2015). This nonlinear data-driven model allows us to perform an attribution activity by NNs in which we consider data of external forcings (solar radiation, volcanic activity, GHGs and black carbon (BC) radiative forcings (RFs), sulphates forcing) as inputs/predictors and the AMO index as target/predictand. Once the forced nature of AMO in the past is discovered, for the first time (at our knowledge) this model allows us to predict its future behaviour in specific scenarios, too.
Looking at the title of this paper, in our approach we would understand if AMO is driven by external forcings or is a manifestation of internal variability (possibly driven by other natural patterns and cycles): in this latter case, there is no direct effect of external forcings, even if this could be mediated by other variability patterns influenced by them. Thus, the possible influence of NAO on AMO is not explicitly considered, but only evaluated on the residuals of our NN reconstruction by external forcings: see Section 4. In particular, this choice allows us to perform also future projections of AMO, basing on several scenarios of forcings (obviously, nobody can envisage future scenarios for NAO).
Incidentally, we must stress that in the AMO index adopted here (Enfield et al. 2001) (see also Section 2) the signal of Atlantic SSTs has been linearly detrended in order to 'clean' it from the strong influence of the GHG forcing, but, first of all, this forcing is not really linear, and secondly there are other forcings (both anthropogenic and natural) which can influence this variability pattern.
In general, attribution studies are extensively performed by ensemble runs of GCMs. However, it has been recently shown (Pasini and Mazzocchi 2015;Pasini et al. 2017;Mazzocchi and Pasini 2017) that in a complex system such as the climate an alternative and complementary strategy can be adopted in order to test the reliability of attribution results and possibly increase their robustness. In particular, NNs and other data-driven models (as completely independent models) can be applied to attribution studies. In the case of attribution of the recent global warming, their results corroborate those obtained by GCMs, and permit to explore new aspects of the problem, too. In the case of attribution of AMO behaviour, these new models could clarify its forced or non-forced nature.

Data
In this investigation, we base on annual data related to the AMO index and to the following radiative forcings: the warming part of the anthropogenic forcing RFWARM (GHGs + BC RFs), its cooling part given by sulphates (RFSOX), RF from solar activity (RFSOLAR) and RF from volcanoes (RFVOL).
In particular, we consider the AMO index calculated with the method of Enfield et al. (2001) and download this record from www. esrl. noaa. gov/ psd/ data/ times eries/ AMO. Data about anthropogenic radiative forcings are downloaded from the freely available dataset collected at http:// www. stern davidi. com/ datas ite. html, referred to a paper by Stern and Kaufmann (2014). Here the data on GHG concentrations are taken from the NASA/GISS website and the related RF calculations are performed through standard formulas (Ramaswamy et al. 2001;Kattenberg et al. 1996). Data about sulphates come from the global estimates of their emissions in the past (Smith et al. 2011;Klimont et al. 2013) and the calculation of direct and indirect RFs as in Stern and Kaufmann (2014), which is based on slight modifications of previous studies (Wigley and Raper 1992;Boucher and Pham 2002). These last data are available until 2007: in order to consider a prolonged time series but not to do a too long extrapolation, we continue this series of RFSOX with constant data until 2011. Past data about RF of black carbon come from the RCP8.5 scenario (Meinshausen et al. 2011).
As far as the data of natural radiative forcings are concerned, solar irradiance is approximated by an index previously assembled (Lean 2000) and available at https:// data. giss. nasa. gov/ model force/ solar. irrad iance/. The conversion from solar irradiance to RFSOLAR is obtained in a standard way (Kattenberg et al. 1996). Volcanic activity of dust emission is considered through optical thickness data (Sato et al. 1993) and downloaded from https:// data. giss. nasa. gov/ model force/ strat aer/; RFVOL is − 27 times the optical thickness (Stern and Kaufmann 2014).
In this paper, we consider also two indices of natural variability: The Southern Oscillation Index (SOI), related to ENSO (Ropelewski and Jones 1987;Allan et al. 1991;Können et al. 1998), which is a well-known factor influencing the interannual variability of the global T series. The data are available at www. cru. uea. ac. uk/ cru/ data/ soi/ soi. dat; The NAO index (Jones et al. 1997), downloaded from https:// cruda ta. uea. ac. uk/ cru/ data/ nao/ Interannual variability at the surface could be affected by changes in ocean heat content, too. Thus, we use also data of the quantity of heat stored in the ocean until the depth of 700 m (OHC700), from http:// www. nodc. noaa. gov/ OC5/ 3M_ HEAT_ CONTE NT/ basin_ data. html: see Levitus et al. (2012) for technical details.
Finally, Coupled Model Intercomparison Project Phase 6 (CMIP6) scenarios are considered here as future estimates of the anthropogenic forcings (Eyring et al. 2016;O'Neill et al. 2016;Gidden et al. 2019;Smith et al. 2020). In particular, we based our analysis on recent data from the model REMIND 2.1 of the Potsdam Institute for Climate impact research (Leimbach et al. 2010), downloaded from https:// data. ene. iiasa. ac. at/ ar6/#/ docs

Methods
Feedforward NNs (Hertz et al. 1991;Bishop 1995) are models which permit to perform multiple nonlinear regressions. In our case, we can assess the possible influences of external forcings-considered inputs of the networks-on the AMO behaviour (the target to be 'approached' by the networks' output). But in a framework of short records available such as our case, one has to avoid overfitting problems and handle the sources of variability (e.g. the random choices for initial weights and validation sets). Our tool (Pasini 2015), adopted here, properly maximises the information for training, minimises overfitting problems and 'averages away' the NN model variability by multiple ensemble runs. This tool has been recently applied to several climate-driven problems: see, for instance, Pasini and Modugno (2013), Pasini et al. (2017), Pasini and Amendola (2019), Pasini et al. (2020). Here, we briefly sketch the main characteristic features of the tool, leaving further details to its presentation paper (Pasini 2015).
The NNs considered in this tool are feedforward with one hidden layer, endowed with hyperbolic-tangent transfer functions at the hidden level and a linear function at the output neuron. It is worthwhile to stress that in the present study the training technique of this tool is modified by considering a quasi-Newton backpropagation method-the so-called Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (LeCun et al. 1998)-which is more suitable for handling small datasets with respect to the standard backpropagation adopted previously. This leads to differences in AMO attribution results when compared with those of our previous attempt .
In order to obtain a nonlinear relationship between inputs and targets (to be approximated by outputs) which can be generally valid, a training-validation-test procedure must be applied. Usually, the free parameters of the NNs (i.e. the connection weights) are fixed on the training set by stopping the training phase when the error on the validation set begins to increase; then, the generalisation performance is measured on a third set, unknown to the NNs, the so-called test set. In our tool, developed for analysing small datasets, a generalised leave-one-out procedure is adopted for training, validation and test, which works as follows.
In a dataset formed by annual inputs-target patterns, starting from the first year we extract a single inputs-target pattern from the total available dataset and consider it as test set. Then, a validation set is randomly chosen and the remaining patterns constitute the training set. At the end of a NN run, the connection weights are fixed and we obtain a transfer function from inputs to output. But this result can be influenced by the specific random choice of the initial weights and of the members of the validation set. Thus, multiple runs are performed in an ensemble approach by choosing different random values for weights and members of the validation set. This gives the chance to calculate the ensemble mean of the outputs and 'average away' the intrinsic variability of NN results.
After this reconstruction of the first (annual) value of the target, the procedure can be followed also for the other years; each of their patterns becomes, sequentially, the test set. In this way, one can achieve the estimation of all output values and ensemble means at the end of this generalised leaveone-out procedure.
In this paper, we use also standard multilinear regressions in order to understand the specific value added by the use of NNs in our analysis. In doing so, the same approach to training is adopted: for each single pattern, the coefficients of the linear regression are fixed on the other data (the union of the training and validation sets used in the NN method). Obviously, an advantage is given to the linear model, so that, as we will show in the following section, the better performance of the NN model is even more notable. Of course, in the linear model no ensemble strategy is required.

Analysis and attribution results
In the application of the tool to our problem, the validation sets are composed by the 10% of the total available patterns and 20 multiple runs are performed in our ensemble approach.
Furthermore, in this paper on the attribution of AMO behaviour we follow the rationale of dynamical-modelling attribution via GCMs. This way of acting is well represented in the case of attribution for global T, where one chooses an ensemble of validated GCMs and simulates the behaviour of global T by supposing that certain forcings remain fixed at their preindustrial levels. Analogously, here we take an ensemble of NN models that are able to well reconstruct the time series of the AMO index once all observed values of forcings as inputs are considered, then apply their transfer functions (the validated NN models) to new inputs, in which we mimic the fact that some of their values show no trend since the preindustrial period, finally obtaining new outputs in terms of a 'simulated' AMO index. In doing so, one can appreciate the roles of the real changes in the different forcings on the AMO behaviour.
As cited above, each model ensemble is composed by many NNs. Here we adopt 20 NNs which supply reconstructions of the AMO index values for each year considered test set in an iterated generalised leave-one-out training-validation-test procedure. The networks are endowed with 4 inputs (representing the anthropogenic and natural forcings for each year: RFWARM, RFSOX, RFSOLAR, RFVOL), 4 hidden neurons in a single hidden layer and one output (to be compared with the annual value of the AMO index, i.e. the target). The small architecture of the networks and the specific training method adopted here allow us to avoid overfitting.
The first result (illustrated in Fig. 1) shows that the ensemble mean (blue line) of NN outputs is able to reconstruct the 'shape' of the AMO curve when all the inputs are considered with their real observed values, with a high variance explained (R 2 = 0.591, RMSE = 0.116). At the same time, a multilinear regression model has a quite poor performance (R 2 = 0.265, RMSE = 0.156) and does not show the AMO cycle almost at all: see Fig. 2. In general, even if the ensemble mean tends to average away nonlinear influences, the NN approach clearly shows that nonlinearities between forcings and AMO are crucial to reconstruct its behaviour, which appears to be strongly driven by these external forcings.
The spread of single NN runs (red lines) in Fig. 1 generally spans a range of less than ± 0.1 °C, even if a high variability is visible in 1884, a year following a big volcanic eruption. This is understandable because of the 'spot-like' time series of RFVOL, so that the reconstruction of this year requires a big extrapolation from the values of RFVOL presented to the NNs in the training set.
It is worthwhile to note that the ensemble mean of NN outputs well follows decadal and interdecadal 'modulations' of the AMO curve, but interannual variability is not very well caught. This high-frequency variability could be due to other influences not included in our forcings, or to natural climate variability which manifests itself at higher frequencies. Here we consider the time series of SOI and NAO indices and try to see if they could explain a part of the residuals, i.e. the time series of observed AMO minus reconstructed AMO (ensemble mean) since 1866, initial date of SOI estimations. Then, since 1955, initial date of OHC700 estimations, we consider also the insertion of a linearly detrended OHC700 index as predictor in our NNs. The results (not shown) are quite poor: 20 NNs endowed with SOI and NAO as inputs and other 20 NNs with SOI, NAO and OHC700 as inputs show ensemble mean performance in terms of R 2 of less than 0.1. We hope to do a more accurate investigation on interannual variability of AMO in a future work.
Once this good reconstruction of the AMO oscillations in the recent past is obtained, a proper attribution activity can start to understand which forcings are the most influential in driving its behaviour. We mimic the situations in which each forcing in turn does not exhibit any trend from the beginning of our time series and consider these new values as inputs to the NNs endowed by the weights fixed through the previous reconstruction stage. In doing so, only one forcing in turn is kept fixed or stationary; the other ones are considered with their real values.
If RFVOL is fixed to the low initial value of its record, the final reconstruction result is very similar to the full reconstruction runs: see Fig. 3 (and also  values of two indices of reconstruction performance and an index of statistical significance are reported). The 'shape' of AMO record is completely preserved. This is not surprising due to the 'spot-like' nature of this forcing, which influences just few years of our time series.
As far as the values of RFSOLAR are concerned, they show the peculiar 11-year cycle but also a clear increase in the period 1910-1950, a transition from a low-power regime to a high-power one. Thus, in our attribution runs we cannot consider a constant value for RFSOLAR, but create a synthetic stationary series (RFSOLSTAT) which shows the 11-year cycle but no more this transition. This is done by a first-order Fourier series built on the first decades of data. See Fig. 4 for the plot of these two series.
When RFSOLSTAT is considered and transferred through the structure of the networks (with weights fixed by the previous reconstruction step) towards the output, the new reconstruction substantially shows a distorted oscillation, with progressively underestimated values of the AMO index starting from the first decades of the last century. This preserves the oscillating characteristic period of the curve and affects only the single values of the last 100 years (see Fig. 5). The performance indices in Table 1 show however a quite consistent decrease.
Thus, the features of the AMO time series are not substantially affected by changes of the natural forcings in the   But, what about the influence of anthropogenic forcings? Through the same method described above, we can evaluate the driving role of anthropogenic forcings on AMO behaviour. If we consider RFWARM fixed at its initial value and constant throughout all the period, a new attribution experiment can be performed. The results (presented in Fig. 6 and Table 1) show that the indices of performance decrease more strongly than in the case of RFSOLSTAT, the oscillation between 1920 and 1970 shows a lower maximum and the final increase in AMO index after 1970 disappears.
Finally, we consider RFSOX constant at its value of the initial year of our time series. The results of this attribution experiment are shown in Fig. 7 and Table 1. The indices of performance show a strong decrease, but the most impressive result comes from Fig. 7: the oscillating behaviour of AMO is now completely destroyed.
Then, if sulphates are considered constant and fixed at the initial value of their record, the oscillating signal of AMO completely disappears, so claiming to a big role of these anthropogenic aerosols for driving AMO behaviour. At the same time, however, when GHGs are considered constant at their preindustrial values only a small oscillation is visible in the simulated AMO signal.
We can conclude this attribution investigation with some clear results: as a matter of fact, the AMO signal is modulated especially by the sulphate forcing, but a component of GHG forcing is needed for determining well its observed form. No substantial role has been detected for external natural forcings as far as the oscillating characteristic features of the AMO time series are concerned.
In short, reminding the question asked in the title of this paper and following our investigation, for the case of AMO natural variability seems not to be really natural.

Future scenarios and projections
In our basic NN runs with real values for external forcings, we obtained 'transfer laws' which strictly link these external forcings to AMO behaviour in the past (via the weights of the NNs fixed at the end of our generalised training-validation-test procedure). Now it is interesting to investigate what happens when these relationships are used for predicting the future. We can do so in a very simple manner with a forward step, that is by propagating the values of external forcings coming from future scenarios towards the output of the networks, which now simulates the annual value of AMO index in the future.
For anthropogenic forcings, we use their future values contained in the CMIP6 data and the Shared Socioeconomic Pathways (SSP) scenarios (Eyring et al. 2016;O'Neill et al. 2016;Gidden et al. 2019;Smith et al. 2020). For natural forcings, we consider a future with a constant low value of RFVOL (it is difficult to hypothesise frequency and strength of future eruptions) and two different scenarios for the solar activity, as shown in Fig. 8. In the first scenario, we consider that until the end of this century RFSOLAR shall remain in a high-power regime, and in the second one, we hypothesise its return to the values of the low-power regime which characterised the period before 1910: this is done by a first-order Fourier series with a superimposed decreasing trend up to 2050 and a return to values comparable at the low-power regime of the first 60 years  of the observed dataset for the last 50 years (2050-2100).
We consider the structure of the networks fixed in the last step of our procedure, when the value of the AMO index for 2011 was reconstructed. Then, we apply these transfer functions to the new inputs coming from the cited scenarios. The results are as follows.
Our NN ensemble mean projections under four SSP scenarios are presented in Fig. 9 for the high-power solar  Fig. 10 for the low-power one. In the first case, after a slow initial increase, a quasi-constant behaviour of the AMO index is visible for SSP2-4.5 scenario and a slight decreasing tendency for SSP1-2.6, while SSP3-7.0 and especially SSP5-8.5 show clear increasing tendencies (see Fig. 9). In any case, the oscillating AMO behaviour seems to be destroyed in the future.
If the low-power solar scenario is considered (Fig. 10), SSP2-4.5 shows still a quasi-constant behaviour, while SSP3-7.0 and SSP5-8.5 are still increasing, even if more slowly than in the previous case. Only SSP1-2.6 appears more sensitive to this change in the solar forcing and, starting from the 2030 decade, its projection shows a decrease which, however, does not achieve the minimum values of the AMO time series. Furthermore, the decreasing tendency is quasi-linear and does not show any periodic cycle.
Thus, the future behaviour of AMO shows no periodic oscillations in our NN projections and seems to be only partially affected by changes in the solar forcing. This last consideration is consistent with the results obtained in our reconstructions for the past, which show the little impact of changes in natural forcings on AMO behaviour and lead to the result that its major driver was instead the radiative forcing of sulphates. In the SSP scenarios, RFSOX is assumed to decrease in absolute value from the present radiative forcing of about − 1.0 to − 0.5 W/m 2 in 2100, almost linearly (the negative sign of RFSOX influence indicates its cooling effect). This is due to the realistic hypothesis that the environmental legislation (which regulates sulphates emissions, dangerous for health) will be applied more strictly in the future.
We have no doubt that this hypothesis has a high probability to be correct, but, just to show a counter-example, one can hypothesise a different tendency for sulphates in the future and look at what can happen to the AMO index. Suppose that RFSOX shall pass linearly from the present value of − 1.0 to − 1.4 W/m 2 in 2100 and put these data in input to our NNs; the projection results are shown in Figs. 11 and 12 for the hypothesised future cases of a high-power solar regime and a low-power one, respectively.
In Fig. 11, the behaviours of SSP5-8.5 and SSP3-7.0 show increasing values of the AMO index, but with Fig. 9 NN ensemble mean projections of the future behaviour of the AMO index under four SSP scenarios (SSP1-2.6, blue line; SSP2-4.5, yellow line; SSP3-7.0, brown line; SSP5-8.5, red line) in the case of a high-power solar radiation hypothesised for the future Fig. 10 As in Fig. 9, but in the case of a low-power solar radiation hypothesised for the future reduced increasing rates in the last half of the century (and an almost null rate in the last decades); SSP2-4.5 projections show small increases followed by small decreases. Finally, SSP1-2.6 shows an initial linear tendency to decrease after 2030, with a plateau in the last decades which achieves the minimum values of the historical AMO, mimicking the reappearance of a cycle, even if with a period longer than that experienced in the past.
In Fig. 12, while in the two highest SSPs an initial tendency to a slow increase and a final asymptotic behaviour are visible, SSP2-4.5 shows a decrease after 2050; the projections of SSP1-2.6 show a more rapid decrease and achieve the minimum values of the historical AMO, also showing the completion of a cycle, with a period which resembles that of the historical AMO.
Thus, in our NN model the only way to recover an AMO cycle in the future is to hypothesise a big increase of sulphate emissions in the atmosphere and the correspondent strong negative radiative forcing under a strong mitigation scenario. In the meantime, a low-power regime of the Sun must be hypothesised, too, if we want to recover a cycle with a period resembling that of the historical AMO.

Conclusions
In this paper, a data-driven attribution investigation has been performed for the AMO behaviour in the past via a machine learning technique, NN modelling. The results show that a NN model with the values of anthropogenic and natural forcings as inputs permits to reconstruct the oscillating characteristic features of AMO. The following NN attribution study-which follows the way of acting of more classical attribution ones via GCM ensembles-shows that this oscillating behaviour strongly depends on anthropogenic sulphates.
Once the forced nature of AMO in the last 150 years is seen, the original application of our model to future predictions of the AMO behaviour shows that it shall probably lose its oscillating characteristic features. The only way to Fig. 11 NN ensemble mean projections of the future behaviour of the AMO index under four SSP scenarios, when RFSOX passes linearly from the present value of − 1.0 to − 1.4 W/m. 2 in 2100. In this case, we consider a high-power solar radiation scenario for the future (SSP1-2.6, blue line; SSP2-4.5, yellow line; SSP3-7.0, brown line; SSP5-8.5, red line) Fig. 11, but when we consider a low-power solar radiation scenario for the future recover them is to consider an unrealistic increase in anthropogenic sulphates in the future under a strong mitigation scenario, and possibly a low-power solar regime.

Fig. 12 As in
In the framework of the recent literature on the attribution of AMO behaviour (widely described in Section 1), the present paper explores this topic by a data-driven method, which is complementary to more classical dynamical studies by GCMs and permits to investigate the unforced or forced nature of AMO in an independent way. This can be important in itself, as the behaviour of complex systems often benefits of analyses from different 'angles', especially for comparing the new results to previous ones achieved with different methods. In particular, our study corroborates the previous dynamical studies which evidenced the importance of anthropogenic forcings (sulphates, in particular) for the reconstruction of the past behaviour of AMO. Here, the oscillating behaviour is completely lost if we do not consider the real changes in sulphate RF.
Furthermore, our method allows us immediately to predict the future behaviour of AMO in different SSP scenarios, achieving the original result that its oscillating behaviour shall be lost in each scenario, if we exclude an artificial scenario (built ad hoc by us) in which anthropogenic sulphates shall show a big increase in concentration and the solar RF possibly shall present a low-power regime.
Due to the established influence of AMO on climate and meteorological phenomena in several regions of the world (see, for instance, Knight et al. 2006), the results presented in this paper can be important to better understand the past and envisage several future scenarios.