Investigating the impact of CO2 on low-frequency variability of the AMOC in HadCM3

This study investigates the impact of CO 2 on the amplitude, frequency, and mechanisms of Atlantic me- ridional overturning circulation (AMOC) variability in millennial simulations of the HadCM3 coupled climate model. Multichannel singular spectrum analysis (MSSA) and empirical orthogonal functions (EOFs) are applied to the AMOC at four quasi-equilibrium CO 2 forcings. The amount of variance explained by the ﬁrst and second eigenmodesappearsto be small(i.e., 11.19%); however, the results indicate that both AMOC strength and variability weaken at higher CO 2 concentrations. This accompanies an apparent shift from a predominant100–125-yr cycleat 350ppm to 160yr at 1400ppm. Changes in amplitude areshown to feed back onto the atmosphere. Variability may be linked to salinity-driven density changes in the Greenland–Iceland– Norwegian Seas, fueled by advection of anomalies predominantly from the Arctic and Caribbean regions. A positive density anomaly accompanies a decrease in stratiﬁcation and an increase in convection and Ekman pumping, generatinga strong phaseof the AMOC (andvice versa). Arctic anomalies maybe generated via an internal ocean mode that may be key in driving variability and are shown to weaken at higher CO 2 , possibly driving the overall reduction in amplitude. Tropical anomalies may play a secondary role in modulating variabilityand arethoughtto be moreinﬂuentialat higherCO 2 , possiblydue to an increasedresidencetimein the subtropical gyre and/or increased surface runoff driven by simulated dieback of the Amazon rain forest. These results indicate that CO 2 may not only weaken AMOC strength but also alter the mechanisms that drive variability, both of which have implications for climate change on multicentury time scales.


Introduction
The Atlantic meridional overturning circulation (AMOC) is an important system of ocean currents that transport significant quantities of heat across the North Atlantic. It has important implications for Northern Hemisphere climate, with the potential to influence surface air temperatures (SATs), sea surface temperatures (SSTs), and precipitation (e.g., Vellinga and Wu 2004;Knight 2005;Frankcombe and Dijkstra 2010;Delworth and Mann 2000).
Both observational (Clark et al. 2002) and modeling studies (Swingedouw et al. 2014) have shown the potential for the AMOC to fluctuate on strength on a range of time scales. In the absence of long observational datasets (Wunsch and Heimbach 2006), climate models are used in order to understand longer time scales of variability (e.g., von Storch et al. 2000). Simulations forced with current and preindustrial concentrations have yielded a range of spectral time scales from decadal to centennial depending on the model used (see Swingedouw et al. 2014). A number of mechanisms have been attributed to driving this variability. Models such as IPSL (Msadek and Frankignoul 2009), CCSM (Danabasoglu 2008;Danabasoglu et al. 2012), and ECHAM5 (Timmermann et al. 1998;Zhu and Jungclaus 2008) indicate a two-way coupling between the atmosphere and ocean. Studies using the GFDL model (Delworth et al. 1993;Delworth and Greatbatch 2000;Delworth et al. 1997) and ECHAM5-MPI-OM (Jungclaus et al. 2005) concluded that fluctuations are oceanic in origin but forced by the atmosphere. External forcing events such as volcanic eruptions (Otterå et al. 2010;Zanchettin et al. 2012;Iwi et al. 2012;Stenchikov et al. 2009) and salinity anomalies induced by El Niño (Mignot and Frankignoul 2005;Schmittner et al. 2000) have also been suggested.
The potential for the AMOC to influence regional and global climate has prompted a range of modeling studies investigating how it may respond to anthropogenic climate change. These have commonly focused on the impact of CO 2 on overall AMOC strength, with a number of studies showing that higher concentrations will weaken the AMOC, albeit to varying degrees ). This has primarily been linked to a reduction in surface density in sinking regions, via surface heat loss and/or a change in the freshwater flux (e.g., Gregory et al. 2005;Swingedouw et al. 2007; Thorpe et al. 2001;Dixon et al. 1999;Mikolajewicz and Voss 2000;Bakker et al. 2016b). In contrast, there has been less focus on the possible effects of CO 2 on the amplitude and frequency of AMOC variability.
This study will investigate how CO 2 concentration impacts low-frequency variability of the AMOC in the HadCM3 coupled climate model. Although there are clear limitations to using one model, it does allow thorough investigation of the potential mechanisms that would not be possible with a multimodel study. Despite a number of drawbacks (see section 2), HadCM3 has been shown to simulate the AMOC, ocean heat transport and aspects of the freshwater oceanatmosphere cycle accurately compared to observations (Gordon et al. 2000;Pardaens et al. 2003). Furthermore the resolution of the model permits millennial-scale simulations in order to study low-frequency variability.
Previous studies with HadCM3 have highlighted the relationship between density anomalies in the Greenland-Iceland-Norwegian (GIN) Seas and AMOC variability. Variations in density have commonly been linked to salinity anomalies, which are advected into the region primarily from the Arctic and/or Caribbean regions. Vellinga and Wu (2004, hereafter VW04) identified a coupled ocean-atmosphere mechanism in the tropics, with periods of strong AMOC driving an increased equatorial SST gradient, a northward shift in the position of the ITCZ, and an increase in precipitation north of the equator. This generates a fresh anomaly that is advected north into the region of downwelling initiating a reversal in AMOC strength and vice versa. This observation follows earlier work by Manabe and Stouffer (1997), who showed via a hosing experiment that the AMOC was sensitive to subtropical freshening. Menary et al. (2012) showed that this mechanism was also present in the Kiel Climate model and highlighted proxy records that give a comparable periodicity of between 50 and 150 yr as evidence for the mechanism. Hawkins and Sutton (2007, hereafter HS07) used 3D empirical orthogonal functions (EOFs) to highlight the importance of the Arctic, which releases salinity anomalies into the GIN Seas, increasing convection and possibly driving an increase in AMOC strength. The formation of these anomalies was hypothesized to be a result of an internal ocean mode controlled by salinity gradients between the Arctic and the GIN Seas, with southward advection in Arctic surface waters compensated by deep northward advection of an opposing anomaly from the GIN Seas.
Jackson and Vellinga (2013, hereafter JV13) concluded that both the mechanisms of VW04 and HS07 might be present, with the Arctic playing a more predominant role and the Caribbean region modulating variability when it is present. In contrast to HS07, they deduced that Arctic anomaly formation is a coupled atmosphere-ocean feedback, responding to stochastic sea level pressure (SLP) variations. This ageostrophically alters the geostrophic balance in the Beaufort Gyre, driving upwelling and advection of salinity into the GIN Seas. A perturbed physics ensemble that altered surface fluxes and winds showed that the background climate might also be important in determining both the frequency and amplitude of variability.
Although the amplitude, frequency, mechanisms, and impacts of AMOC variability have been well researched, there has been less focus on how these may be altered at higher CO 2 concentrations. A number of previous studies that have investigated time-dependent dynamic systems such as the AMOC have focused on the framework introduced in the Stommel box model (Stommel 1961), which indicates that a decline in the strength of the system, such as that induced by higher CO 2 , would increase instability, pushing the system closer to a point of collapse (i.e., a bifurcation point) and so increase variability (Tziperman 1997;Tziperman et al. 1994;Wiesenfeld and Mcnamara 1986;Ditlevsen and Johnsen 2010;Rahmstorf 1995;Monahan 2002). However, results from MacMartin et al. (2016) were not consistent with this hypothesis. They investigated the impact of CO 2 on the AMOC at a range of latitudes in two simulations of GFDL-ESM2M. They showed that at 4 3 CO 2 a decrease in AMOC strength was accompanied by a decrease in midlatitude variability, with the region of greatest variability shifted to a limited area at higher latitudes. This was linked to a change in ocean dynamics, with an increase in stratification possibly reflecting a weaker and shallower mean overturning at low latitudes.
In this study we aim to answer two key questions: 1) Does CO 2 impact the amplitude/frequency of AMOC variability? 2) What mechanisms are responsible for these changes? Section 2 outlines the model and methods and section 3 gives an overview of the climatology and impact of CO 2 on AMOC strength. Section 4 investigates spatial and time series variability, and section 5 discusses the climatic implications of variability. Section 6 examines the possible mechanisms before a discussion and summary is presented in section 7.

a. Model
The Hadley Centre Climate Model, version 3 (HadCM3), is a coupled Earth system model consisting of 3D dynamical atmosphere and ocean components (Valdes et al. 2017). The atmospheric component has a horizontal resolution of 3.758 3 2.58 with 19 vertical levels. The ocean has a horizontal resolution of 1.258 3 1.258 with 20 vertical levels that have a finer resolution toward the surface. The ocean uses a ''rigid lid'' approach so volume remains constant, with runoff from the land converted to salinity flux where river outflow points reach the ocean. The mixed layer model in the ocean is that of Kraus and Turner (1967) with vertical mixing of tracers using a simplified version of the Large et al. (1994) scheme, which is discussed in detail in Gordon et al. (2000). Horizontal mixing of tracers (e.g., salinity) via eddies uses the Gent and McWilliams (1990) isopycnal parameterization with isopycnal mixing using the Redi (1982) scheme implemented by the method of Griffies et al. (1998). HadCM3 does not use flux adjustment, which can impact deep-water formation and influence the AMOC (Marotzke and Stone 1995).
The drawbacks of the model include issues related to ocean resolution that affects some significant channels including the North Atlantic overflows (Roberts et al. 1996). The Canadian Archipelago, including the Hudson and Davis Straits, is closed and there is zero barotropic flow through the Bering Strait and so no net volume transport. This is in disagreement with observations that indicate that freshwater import into the Arctic via the Bering Strait is significant (Cattle and Cresswell 2000). To compensate for this, all freshwater export from the Arctic is via the Atlantic sector into the Nordic seas. This may have implications for the AMOC as the Bering Strait region has been shown to moderate AMOC strength (Hu and Meehl 2005). The Denmark Strait and Greenland-Scotland ridge have been deepened in order to produce a mean outflow that matches observations (Gordon et al. 2000). Other problems include a weakness in the wind stress in the North Atlantic storm track, which potentially impacts gyre strength (Gordon et al. 2000). Furthermore, peak flow of the AMOC in HadCM3 is at approximately 800 m (see Fig. 4), shallower than the 1000 m observed by RAPID-MOCHA (Smeed et al. 2015). The cause of this is likely to be surplus surface salinity in the North Atlantic, driven by an excess in net evaporation in the topics and subtropics and to a lesser extent insufficient subtropical runoff (Pardaens et al. 2003). Consequently the Atlantic is too stratified and stable, reducing the depth of maximum overturning. This has also been shown to be the case in other models, including GFDL CM2.1, NCAR CCSM3, and the MPI models (Roberts et al. 2013;Msadek et al. 2013).
The results presented here are from 2000-yr quasiequilibrium simulations at four CO 2 concentrations; 350, 700, 1050, and 1400 ppm. We will refer to these experiments as 1x, 2x, 3x, and 4x, respectively. All analysis is conducted on the final 1000 yr of each simulation, a time period that is suitable for examining low-frequency variability. It also permits a 1000-yr spinup in order for the ocean to reach a relative state of equilibrium as indicated by the volume averaged upper ocean temperatures (not shown). The full time series of the AMOC index (i.e., mean AMOC strength between 408 and 508N at 800 m; see section 4b) for each experiment are shown in Fig. 1, showing that the AMOC has reached a relative state of equilibrium in each simulation after 1000 yr. To remove a small climate drift the AMOC streamfunction and climate variables have been detrended by subtracting a liner least squares fit prior to analysis.

1) SINGULAR SPECTRUM ANALYSIS
To analyze variability of the AMOC we have employed singular spectrum analysis (SSA) and its multivariate counterpart (MSSA). This is a nonparametric technique used to decompose a time series into a range of periodicities and noise in order to extract a signal of interest and more clearly understand the dynamical nature of the system being studied. The MSSA technique can be used to analyze gridded datasets. The methodologies and applications of SSA and MSSA are outlined in Vautard et al. (1992), Elsner and Tsonis (1996), and, more recently, Golyandina and Zhigljavsky (2013), with a thorough mathematical overview for the application to climatic time series in Ghil et al. (2002). Numerous studies have applied SSA and MSSA to climatological data (e.g., Ghil et al. 2002;Plaut et al. 1995;Moron et al. 1998), including the AMOC (Alvarez- Garcia et al. 2008;Huang et al. 2012). The SSA and MSSA in this paper have utilized the R package ''Rssa,'' which is extensively outlined in Golyandina and Korobeynikov (2014) and Golyandina et al. (2015).
The SSA algorithm consists of two basic stages, decomposition and reconstruction. First a covariance matrix is generated from the time series with a specified window length representing the number of lags M. Consequently the covariance matrix is of dimensions M 3 M. MSSA differs from SSA in the way this trajectory matrix is constructed (see Golyandina et al. 2015). The matrix undergoes singular value decomposition (SVD) to produce corresponding eigenvectors and eigenvalues. The projection of the original time series onto the EOFs yields the principal components (PCs), with the first PCs (PC1s) representing the largest amount of variance. These PCs can be ''reconstructed'' and put back on the scale of the original time series by projecting back onto the corresponding eigenvectors. Combination of all the individual reconstructed PCs would generate the original time series; however, partial reconstruction permits the identification of potential oscillatory components and the removal of noise.
An important consideration in SSA and MSSA is the window length M. Ghil et al. (2002) state that the size of M is a compromise: a large M allows more information to be extracted from the original time series; however, it will permit fewer repetitions, which consequently reduces statistical confidence. Other studies state that M should be approximately divisible by the time scale of a potential oscillation in the data (Golyandina 2010;Golyandina et al. 2001). In this study, we tested AMOC variability for a range of window lengths ranging from 100 to 300 yr. The results were similar for each value of M indicating stability in our results; the values used are outlined in sections 5a and 5b. All climate variables used in this study have been decomposed using the MSSA technique with a window length that corresponds to the appropriate CO 2 concentration.

2) SIGNIFICANCE
To calculate the statistical significance of the anomalies and correlations we have used the moving block bootstrap technique (Wilks 1997) with 95% confidence limits calculated using the bootstrap percentile method (see Hall 1988). If the bootstrap confidence interval passes through zero the anomaly/correlation is also deemed insignificant. Power spectra have been generated using discrete Fourier transform to reveal periodicities and their relative strengths. We have tested the significance of these peaks against the 95% confidence level determined by a red noise spectrum of a first-order autoregressive (AR1) process.

Impact of CO 2 a. Climatology
The impact of CO 2 on a range of climate variables is shown in Fig. 2. SATs and SSTs are amplified across the region, the former predominantly at high latitudes and over land and the latter to the north of Scandinavia and east of the United States. There is a complex pattern of precipitation change, with a northward shift and/or weakening of the ITCZ driving a decrease across much of the tropics and an increase at high latitudes. This pattern of temperature and precipitation changes results in an increase in salinity across the tropics and a decrease at the poles. There is a dramatic decrease in sea ice, with almost complete disappearance at 4 3 CO 2 across much of the Arctic. Many of the climate impacts simulated by HadCM3 are consistent with those from other climate models, as shown in the CMIP5 experiments (e.g., Knutti and Sedlacek 2013) and the IPCC Fifth Assessment Report (AR5; Flato et al. 2013).
The mixed layer climatology gives a guide as to the key location of deep-water formation. The GIN Seas region shows the greatest depth of overturning extending to mean depths of 247 m off the coast of Scandinavia. In contrast to other models (such as CCSM), HadCM3 exhibits shallow overturning depths in the Labrador Sea. There is a reduction in the depth of the mixed layer in the GIN Seas of up to 87 m at 4 3 CO 2 . Mean depth profiles for salinity, temperature, and density for the GIN Seas (see Table 1; Fig. 11a) are shown in Fig. 3. For reference these are shown with the NODC Levitus observational dataset spanning 1900-92. Within the model, cooler, fresher, and less dense surface waters overlie warmer, saltier, and denser subsurface waters. Ocean velocity maps (not shown) show that this stratification is likely to be a result of a warmer northeastward flow of water via the North Atlantic Current (NAC), which sits below colder fresher water transported southward by the East Greenland Current (EGC). CO 2 alters these profiles by increasing temperature and density throughout the water column. Salinity decreases at the surface and increases in the subsurface, the former possibly reflecting the modeled increase in precipitation in the high northern latitudes and the latter reflecting saltier Atlantic inflow at higher CO 2 . Furthermore the water column is more stratified, indicated by greater contrast between surface and subsurface temperature and salinity.

b. AMOC strength
The mean Atlantic AMOC streamfunction for 1x and the anomalies due to increasing CO 2 are shown in Fig. 4. At 1x, the upper North Atlantic Deep Water (NADW) cell reaches a peak of 16.6 Sv (1 Sv [ 10 6 m 3 s 21 ) at 358N at a depth of 800 m. The mean strength compares well to the 16.9 Sv observed by the RAPID-MOCHA (between April 2004 and October 2015 at 26.58N) but is too shallow compared to observations (Smeed et al. 2015). There is a reduction in AMOC strength across this cell of up to 5.0 Sv (30.2%) at 2x, 4.95 Sv (29.8%) at 3x, and 5.6 Sv (34%) at 4x, in addition to a potential shallowing in the NADW cell. There is also a decrease in the strength of the northward moving Antarctic Bottom Water (AABW) cell and surface waters (08-308S and 308N) of up to 3.8 Sv.
The decline in AMOC strength has been commonly linked to changes in surface heat loss and/or the freshwater flux that reduces density in downwelling regions. (e.g., Gregory et al. 2005;Swingedouw et al. 2007; Thorpe et al. 2001;Dixon et al. 1999;Mikolajewicz and Voss 2000). Figure 3 shows the reduction in density in the GIN Seas at higher CO 2 . Thorpe et al. (2001) and Thorpe (2005) have previously investigated CO 2 -induced weakening of the AMOC in HadCM3. They concluded that both CO 2 -induced warming and a decrease in salinity at high latitudes drive this change, with the former playing a more prominent role accounting for approximately 60% of weakening.
It is also important to highlight the increase in stratification in the GIN Seas at higher CO 2 (Fig. 3), which would consequently decrease convection and may weaken overturning. This is highlighted by the Brunt-Väisälä (or buoyancy) frequency for the top FIG. 3. Mean depth profiles for potential temperature (solid; 8C), salinity (dashed), and density (kg m 22 ; dotted) in the GIN Seas (see Table 1) for the different simulations. The NODC Levitus observational dataset is included for reference, showing mean values spanning 1900-92 (data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, from their website at http://www.esrl.noaa.gov/ psd/). The GIN Seas region in the Levitus data is defined as 60.6258-79.3758N, 13.758W-11.258E. Large values indicate a greater density gradient, more stable conditions, and a reduced rate of convection. At 1x the lowest values and so a region of weaker stability is found around the GIN Seas region with values increasing at higher CO 2 specifically in the center of the region, indicating greater stability and thus reduced convection. This may act to reduce crossisopycnal flow and drive the apparent reduction in AMOC strength.

a. Spatial variability
To understand the spatiotemporal patterns of AMOC variance, we have applied MSSA (see section 2b) to the AMOC streamfunction for the region 308S-908N. This procedure was originally outlined in Plaut and Vautard (1994) and subsequently used in Huang et al. (2012). The AMOC streamfunction was first decomposed using MSSA with a range of M values from 100 to 300 yr (see section 2), with similar results produced. At 1x, 2x, and 3x values, the first and second eigenmodes bear the resemblance of a sine and cosine pair with a similar period that is out of phase (i.e., a phase quadrature relationship) and represent an oscillation on the order of approximately 110-125 yr. We chose an M value of 115 yr for 1x and 125 yr for 2x and 3x. At 4x the oscillation for the first and second is on the order of 150-160 yr and we chose an M value of 160 yr. These values for M were also chosen based on the initial spectral analysis of the AMOC index (as shown in Fig. 7 and discussed in section 4b).
The leading two eigenmodes are used to reconstruct the AMOC streamfunction. Together they account for 11.19% (6.86% and 6.69%), 7.69% (5.36% and 4.99%), 5.51% (3.84% and 3.8%), and 3.88% (2.8% and 2.3%) of the total variance for 1x, 2x, 3x, and 4x respectively. Although these values appear to be small, at 1x they are similar to those presented in previous studies that have derived oscillating MSSA modes for the AMOC (Huang et al. 2012) and for SST variability of the world oceans (Moron et al. 1998). The similarity of the first and second modes and 2D graphs of the eigenvectors (not shown) indicate a degenerate pair at all concentrations albeit more weakly at higher CO 2 concentrations. There is also a relatively coherent degenerate pair oscillation for a combination of eigenmodes between 5 and 9 depending on the concentration of CO 2 . This is at a higher frequency and may indicate a weaker oscillation of between 15 and 20 yr, possibly the frequency identified in the study of Dong and Sutton (2005). These eigenmodes will be investigated in a future study. We have applied this MSSA method and these values of M to all climate variables used in this study.
To identify the key areas of variance, EOF analysis is applied to the MSSA AMOC streamfunction. The values of the EOFs and corresponding PCs are arbitrary, so we have regressed the MSSA time series onto the normalized PCs in order to put them on the same scale and allow comparison. The first and second EOFs and their corresponding scaled PCs are shown in 1x, 2x, and 3x, equivalent to 111-yr cycles, and possibly five oscillations at 4x, equivalent to 160 yr.
At 1x, the first EOF explains 87.2% of the variance of the MSSA data, which decreases to 30.8% at 4x. The first EOF maps for all CO 2 concentrations are characterized by a singular deep overturning cell southward of 608N that corresponds to basinwide variability of the AMOC. This peaks between 408 and 508N at a depth of 800-1000 m for all CO 2 concentrations. The northward boundary of the cell is limited by the location of the sill on the southern border of the GIN Seas. North of this there appears to be two cells of different sign at 1x that may correspond to a north-south shift in variability within the GIN Seas and Arctic. The second EOF explains 3.2% of the variance of the AMOC at 1x increasing to 19.5% at 4x. This EOF is more complicated, showing positive and negative anomalies that may reflect variability in both the NADW and AABW cells.
At higher concentrations of CO 2 the depth, southward extent, and strength of variance in the main cell decreases, accompanied by a reduction north of 608N. The region of greatest variance remains at approximately 408-508N and 800-m depth for all simulations, in contrast to MacMartin et al. (2016), who identified a northward shift in variance at 4 3 CO 2 . The PC1 and PC2 time series indicate a decline in the clarity of oscillations and a possible shift toward lower frequency. At 1x and 2x AMOC amplitude and frequency is relatively coherent; however, at 3x there is a decline in amplitude between 375 and 800 yr and at 4x both frequency and amplitude appear relatively intermittent. This would indicate that higher concentrations of CO 2 influence not only AMOC strength and depth of overturning, but also the amplitude, frequency, and coherence of variability.

b. Time series variability
To measure temporal variability, it is common to define an AMOC index [meridional overturning index (MOI)] at a latitude and depth where variability is greatest. Using the first EOFs as a guide, MOIs are calculated as the mean strength of the AMOC between 408 and 508N at a depth of 800 m for all CO 2 concentrations. The MOI was decomposed using SSA (see section 2), using the same values for M (section 4a). As before, the procedure was tested with a range of window lengths and the results remained stable. Initial decomposition shows that the first eigentriple accounts for greater than 99% of the variance and corresponds to the overall trend of the time series. This trend does not contain oscillatory components and represents variability on longer time scales than those of interest in this study. To isolate higher-frequency variability, the overall trend is reconstructed and extracted and the resultant time series decomposed for a second time. The EOFs, eigenvalues, and 2D graphs of the eigenvectors (not shown) can be used to identify eigenvectors that may constitute an approximate sine wave with similar frequencies and a phase shift (see Golyandina and Korobeynikov 2014). The first and second eigenvectors are extracted and reconstructed for all CO 2 concentrations and are shown with the original MOIs and the first and second PCs in Fig. 6. Together they account for 39.53% (21.18% and 18.35%), 31.02% (16.81% and 14.21%), 26.09% (14.33% and 11.76%), and 17.56% (9.15% and 8.41%) of the total variance for 1x, 2x, 3x, and 4x respectively. The correlation values between the MOI and PC1 are strong for each simulation, indicating that the MOI successfully represents variability within the overall AMOC streamfunction.
Power spectra for the original and the decomposed MOI time series are shown in Fig. 7. Spectral peaks are consistent for both MOIs, although variability at higher frequencies has been removed for the SSA MOI spectrum. The peaks are significant at 1x, 2x, and 3x, but at 4x the peak is only just discernible against the 95% confidence interval of red noise. At 1x, the key time scale of variability is on the order of 100-125 yr, in agreement with the studies of VW04, JV13, and Menary et al. (2012). With increasing CO 2 there is a decrease in the power of this variability and an apparent shift to lower frequency at 4x. Spectral peaks are on the order of 112, 125, and 165 yr for 2x, 3x, and 4x respectively. This reiterates what is shown by the EOF analysis; higher concentrations of CO 2 act to weaken and possibly decrease the frequency of AMOC variability.
These results show the potential for CO 2 to influence the variability of the AMOC, a similar conclusion to that found by MacMartin et al. (2016) with the GFDL model. There are, however, disparities in the results between the two models, namely that MacMartin et al. (2016) focused on a decadal 15-20-yr oscillation and showed that the amplitude of variability remained consistent yet shifted northward from 258-508 to 608N at 4x. We do not see this northward shift and see a decline in the amplitude of variability and a possible change in frequency. However, both studies challenge the theory of the Stommel (1961) box model, which stipulates that a weaker and thus less stable AMOC should exhibit greater variability.

Feedback of variability on climate
To understand how the AMOC feeds back on climate, Fig. 8 shows composite anomaly maps for a number of variables during periods of maximum AMOC (AMax) compared to their mean state. Periods of AMax are defined as the years where PC1 (Fig. 6) for each CO 2 concentration is more than one standard deviation above the mean. A 10-yr running mean has been applied to the climate variables followed by decomposition and reconstruction via MSSA. Composite maps showing AMOC minimum (AMin) are not shown but broadly show the opposite pattern.
AMax coincides with a broad increase in SSTs north of the equator compared to the mean, peaking at approximately 0.88C at 1x and 0.38C at 4x. Anomalously warm temperatures are focused within and to the north of the GIN Seas, ranging from 0.68 to 0.88C at 1x. In the Arctic and to the south of the equator, AMax coincides with a decrease in SSTs on the order of 0.18-0.48C, Periods of AMax are associated with a significant northward shift in the position of the ITCZ over the Atlantic at 1x, increasing (decreasing) precipitation north (south) of the equator. This shift was first observed in HadCM3 by VW04, who linked it to the generation of salinity anomalies in the tropics (see section 6b). There is also a weak trend for an increase in precipitation across much of the North Atlantic extending into the GIN Seas. The northward shift in the position of the ITCZ is significantly weakened at higher CO 2 concentrations, specifically at 2x. The anomalously warm SSTs and SATs at AMax are likely to drive the dipole pattern observed in SLP anomalies that display a strong area of low pressure north of the equator and high pressure to the south. The strongest low pressure anomaly is on the order of 42 Pa located within and to the south of the GIN Seas. This low pressure anomaly is also present at 2x but it diminishes in strength and expands into the high northern latitudes at 3x and 4x.

a. Convection in the GIN Seas
The degree of stratification and convection in the GIN Seas region is likely to play a key role in driving AMOC variability. As shown in Fig. 5, there is a reduction in the Brunt-Väisälä frequency in the GIN Seas with increasing CO 2 indicating weaker convection and so a reduced overall AMOC. We can also investigate the degree of convection between periods of AMax and AMin. Figure 9 shows the average depth profiles in the GIN Seas for temperature, salinity, and density anomalies during AMax and AMin, and the GIN Seas in Fig. 3, the pattern of anomalies during AMax is broadly opposite to the mean water profiles, indicating a less stratified water column at AMax. The opposite effect is seen at AMin. This reduction in stratification acts to enhance convection at AMax as shown by an increase in wind stress curl in the GIN Seas (Fig. 10), particularly at 1x. The increase in temperature and salinity and consequent negative SLP anomaly (Fig. 8) acts to drive anticyclonic winds, which, in addition to the Coriolis force, drives a divergence in the Ekman layer and upward Ekman pumping, reducing stratification and subsequently increasing convection. Rather than the initial driving force, this process may be a positive feedback cycle that responds to an initial increase in the AMOC strength and acts to further enhance convection and strengthen the AMOC.
At higher concentrations of CO 2 there is smaller contrast in temperature, density, and salinity between AMax and AMin, particularly in the top 1000 m. As a result the water column is more stratified during AMax, which, as shown in Fig. 5, weakens convection at higher CO 2 . There is also a possible shift in convection to a region just north of Iceland at 2x and 3x, an impact also shown in Fig. 5. This lessened convection may reflect a weaker overall AMOC at higher CO 2 , which in turn reduces the impact on SSTs and SATs at AMax and consequently the increase in Ekman pumping, which is required to drive enhanced variability.
b. Impact of regional salinity anomalies Although Ekman pumping may enhance convection in the GIN Seas that subsequently drives a period of AMax, temperature-and salinity-driven density changes are likely to initiate the shift between AMax and AMin. As highlighted by Thorpe et al. (2001), the mean strength of the AMOC is primarily influenced by temperature-driven density changes. However, previous studies have commonly focused on salinity anomalies as the key driver of variability (Jungclaus et al. 2005;Delworth et al. 1997) including within HadCM3 (Hawkins and Sutton 2008;JV13;VW04;Pardaens et al. 2003;Menary et al. 2012). Figure 11 shows the no-lag correlation of MSSA salinity averaged over the top 666 m with the MOI. At 1x there is strong correlation in the North Atlantic and the GIN Seas and a strong negative correlation in the Arctic and the southwestern North Atlantic basin. At increasing CO 2 , this general pattern remains, with strong negative correlation in the tropics and Arctic, yet there is a weakening in the positive correlation in and around FIG. 11. Correlation (no lag) of MSSA decomposed salinity anomalies with the MOI, for (a) 1x, (b) 2x, (c) 3x, and (d) 4x. In (a) the boxes outline the GIN Seas, Caribbean, EGC, and Arctic regions that are used in the lagged correlation analysis. Only anomalies that are considered 95% confident are shown (see section 2). the GIN Seas. This is also highlighted in the lower panels in Fig. 10, which shows a positive relationship between Ekman pumping and salinity in the GIN Seas at 1x; this weakens at higher CO 2 with an apparent southward shift to a region just north of Iceland. This section will investigate the changing roles of salinity anomalies at higher CO 2 .
As discussed, the studies of VW04, HS07, and JV13 concluded that salinity anomalies are created and advected primarily from the Arctic and Caribbean regions via the EGC and NAC respectively and are the primary method by which salinity anomalies are created in the GIN Seas. JV13 showed via a salinity budget analysis that regional processes such as sea ice fluctuations and surface processes play only a small role. A number of previous studies have highlighted the potential for Southern Ocean processes to impact AMOC variability. Bakker et al. (2016a) showed that fluctuations in the discharge from the Antarctic ice sheets in response to subsurface ocean temperatures influenced centennial-scale climate variability. Studies have also linked AMOC variability to circumpolar wind stress (Toggweiler and Samuels 1993) and to Weddell Sea processes driven by sea ice fluctuations and deep ocean convection Park and Latif 2008;Latif et al. 2013). We performed a lagged correlation analysis of Southern Ocean salinity with the MOI (not shown) but did not find any significant correlations; therefore we will focus predominantly on the Caribbean, Arctic, and GIN Seas region in this study. Figure 12 shows the lagged correlations of salinity with the MOI for the GIN Seas, EGC, Caribbean, and Arctic (Table 1 and Fig. 11). Positive anomalies in the GIN Seas and EGC occur at approximately zero lag (i.e., at AMax), and salinity in the Caribbean region peaks at approximately 30-53 yr prior to AMax and in the Arctic between 37 and 50 yr (Table 2). Regional correlation shows consistent oscillatory patterns, with time periods between peaks of 110 and 130 yr, a comparable time scale of variability to the MOI (Fig. 7). The strength of correlations weakens at higher CO 2 specifically for the Arctic; however, the Caribbean remains relatively consistent throughout.
To determine the time scale of salinity fluctuations, an EOF analysis was applied to the MSSA decomposed salinity in the different regions and the leading PCs extracted and regressed onto salinity. Figure 13 shows the MOI time series and the first PCs for the Arctic and Caribbean regions. Here the PCs have been lagged with the year where salinity in the region reaches a peak prior to AMax as shown in Fig. 12 and Table 2. The corresponding spectral analysis for each region is shown in Fig. 14. Salinity fluctuations within the Arctic have a greater magnitude than in the Caribbean and show strong correlation with the MOI at 1x and 2x that declines at higher CO 2 particularly at 4x. At 3x, the decline in the extent of variability beyond 400 yr is accompanied by a reduction in variance within the Arctic. The power spectra for Arctic salinity show strong similarities with both the GIN Seas and the overall MOI spectra (Fig. 7), indicating that this region has a strong relationship with the overall time scale of AMOC variability. In contrast, the strength of variability within the Caribbean region remains relatively consistent and shows a small increase at 4x. The power spectrum shows a significant peak at FIG. 12. Lagged correlations of MSSA decomposed salinity anomalies with the MOI for regions identified in Table 1 and Fig. 11a, for (top)-(bottom) 1x, 2x, 3x, and 4x. Thick dashed lines highlight 95% confidence for all time series (see section 2). Salinity anomalies lead the MOI for negative lags. TABLE 2. The lag values (yr) for the generation of positive Arctic and Caribbean salinity anomalies prior to AMax. They have been identified from the lead-lag analysis shown in Fig. 12. They have been applied to the PCs for Arctic and Caribbean anomalies in Fig. 13.

Simulation
Arctic lag (yr) Caribbean lag (yr)   1x  37  43  2x  45  30  3x  48  38  4x  50  53 approximately 160 yr at 4x, consistent with the MOI spectra. Variability at lower CO 2 concentrations shows peaks that are compatible with the MOI but are small in comparison to the Arctic. The PC and spectral patterns may indicate that the Arctic is key in setting the amplitude and time scale of variability while the Caribbean plays a secondary role, an idea similarly put forward by JV13. At higher concentrations of CO 2 , weaker variability in Arctic salinity may drive a reduction in the overall amplitude of AMOC fluctuations. In contrast, variability of salinity in the Caribbean remains consistent and may increase, indicating that this region is comparatively more important for driving variability at higher CO 2 .

1) ARCTIC ANOMALIES
To understand more clearly the development of salinity anomalies in the Arctic, Fig. 15 shows the lagged correlation of Arctic salinity with the MOI prior to and following AMax. At 260 yr during periods of AMin, an isolated positive anomaly develops in the Beaufort Gyre in the center of the basin. This anomaly strengthens to a peak at lag 240 yr (as shown in Fig. 12) before being advected out of the gyre as shown in lag 220 yr. At this time, a negative anomaly develops north of Siberia that intensifies and spreads across the Beaufort Gyre at lag zero and 120 yr during periods of AMax. This negative anomaly is advected out of the basin at lag 140 yr and may contribute in pushing the AMOC into a negative phase. At this time, a positive anomaly develops off the Siberian coast that intensifies at lag 160 yr and restarts the cycle. JV13 showed a similar pattern of Arctic salinity anomalies across the Beaufort Gyre in HadCM3. This cyclical generation of anomalies is present but weakens at higher concentrations of CO 2 , concurrent with that shown by the EOF, spectral, and lead-lag analysis.
There remains uncertainty as to how salinity anomalies are generated in the Arctic and consequently how this mechanism may weaken at higher CO 2 . A salinity budget analysis may help to identify the role of different inputs and how they are influenced at higher CO 2 ; however, we do not have the correct diagnostics to carry out this type of analysis. JV13 performed such an analysis on a preindustrial simulation of HadCM3, concluding that both sea ice processes and surface flow play only a minor role in the generation of Arctic anomalies. We cannot state that this is the case at higher concentrations of CO 2 ; however, we do see only trivial fluctuations in sea ice cover between AMax and AMin and almost zero ice cover at 4 3 CO 2 , which may indicate that sea ice plays only a small role.
JV13 linked the formation of these Arctic anomalies to stochastic variations in SLP that alter the geostrophic balance in the Beaufort gyre. A period of high (low) sea level pressure drives an increase (decrease) in anticyclonic wind stress, strengthening (weakening) the gyre and resulting in downwelling (upwelling) and freshening (salinification) in the center of the basin and upwelling (downwelling) and salinification (freshening) at the coasts. We correlated MSSA SLP with the MOI (not shown); however, correlations were weak with no consistent pattern in the results to agree with this hypothesis, specifically at 2x and 3x where the pattern of salinity anomalies remain strong.
The study of HS07 stipulated that fluctuations in the AMOC might act as a lagged positive feedback on generating salinity anomalies in the Arctic. They highlighted an opposing pattern of salinity anomalies seen at depth in the Arctic basin that may be generated by deep advection of an opposing water mass from the GIN Seas. This migrates to shallower depth following advection of the surface anomaly. A lead-lag correlation analysis of the MOI with Arctic salinity at a range of depths is shown in Fig. 16. It indicates a relatively clear out-ofphase relationship between anomalies at depth relative to those at the surface. At 1 3 CO 2 , a positive (negative) anomaly at 5 m (300 m) follow approximately 40 yr (30 yr) later from those at depth. This oscillatory nature is coherent at 1x, 2x, and 3x but weakens at 4x. This weakening is consistent with the idea that this is a positive feedback and a direct consequence of weaker overall variability; that is, weaker variability is due to weaker oscillations in GIN Seas salinity, which in turn FIG. 13. The MOI (black) and PC1 for salinity anomalies from the Arctic (blue) and Caribbean (red), for (top)-(bottom) 1x, 2x, 3x, and 4x. The PCs have been lagged with values identified in the lead-lag analysis in Fig. 12 and shown in Table 2. drives weaker northward advection of anomalies into the Arctic and thus inhibits the formation of surface anomalies that are then advected south and required to maintain variability. It remains uncertain, however, what the overall driver is that weakens this mechanism. It may be a response to an overall reduction in strength of the AMOC; however, further experiments are required in order to test this hypothesis.

2) CARIBBEAN ANOMALIES
The study of VW04 previously linked the generation of negative (positive) anomalies in the Caribbean to the northward (southward) shift of the ITCZ and an increase (decrease) in precipitation north of the equator during AMax (AMin). The resulting anomaly is advected north and pushes the AMOC into an opposing phase. This northward shift in the ITCZ is apparent in the precipitation anomalies in Fig. 8, but it is less pronounced at higher concentrations of CO 2 . Despite this, evidence in Figs. 11-14 indicates that variability in Caribbean salinity remains consistent and may even increase in amplitude at 4x.
A possible driver of this may be linked to the residence time of salinity anomalies in the Caribbean region. VW04 concluded from a tracer experiment that salinity anomalies advected from the Caribbean had a long residence time in the subtropical gyre, with subduction to subsurface waters prior to northward advection into the GIN Seas. Higher concentrations of CO 2 may influence this process in two ways; first the weaker AMOC may slow the rate of northward advection. Second, increased salinity in the subtropical gyre at higher CO 2 is accompanied by an increase in stratification, which may reduce the rate of vertical mixing, a similar impact to that observed in the perturbed physics ensemble of JV13. Consequently the residence time in the subtropical gyre is increased and, because of relatively high evaporation, the salinity anomalies are intensified. This may not only produce the comparatively large anomalies at higher CO 2 , but also the increased advection time may result in the apparent shift toward lower frequencies (Fig. 14).
In addition to this, Caribbean salinity anomalies may also be influenced by changes in surface runoff rate at higher concentrations of CO 2 . Figure 17 shows the surface runoff anomaly during periods of AMax for the different simulations. There is a greater impact at higher concentrations of CO 2 , where there is a positive anomaly toward the north of the continent and a negative anomaly in the central regions. This general pattern at 3x and 4x reflects the change in precipitation rates over FIG. 14. Power spectra of the PC1 of the MSSA decomposed salinity anomalies averaged over the top 666 m, in different ocean regions. These regions are identified in Table 1 and Fig. 11a. Note the change in scale for the Arctic. The dashed lines represent the 95% confidence interval where all spectra are considered significant. this region at AMax as shown in Fig. 8. However, at these higher CO 2 concentrations there is an increase in the proportion of this precipitation that is converted to runoff simulated by the model (Fig. 17b). This may reflect shifting vegetation dynamics across the Amazon region (Fig. 17a), with a change from predominantly broadleaf and C 4 grasses at 1 3 CO 2 to bare soil at higher CO 2 concentrations. This dieback of the Amazon rain forest in HadCM3 is a feedback that has been investigated in a number of previous studies (Betts et al. 2004;Cox et al. 2000;Boulton et al. 2013) and reflects a ''tipping point'' at which drying out and increasing temperatures overcome CO 2 fertilization and inhibit vegetation growth. Removal of vegetation alters the partitioning of water, reducing evapotranspiration and increasing infiltration driving a positive surface runoff anomaly. The river routing scheme (Gordon et al. 2000) channels this surface runoff instantaneously into the ocean via river outflow points.
The key river catchment in this region is the Amazon River, with an outflow point located in northeastern Brazil at approximately the region of strongest negative correlation in Fig. 11. The increase in runoff rate at this point may help to drive the fluctuations in Caribbean salinity that influence AMOC variability at higher CO 2 , despite the reduced impact on the ITCZ. This potential mechanism indicates that in addition to tropical atmosphere-ocean feedbacks, AMOC variability may also be influenced by terrestrial-ocean feedbacks brought about by atmospheric changes (i.e., precipitation). This mechanism is strengthened at higher CO 2 due to Amazon dieback, which may act to counter the weaker role of the Arctic in driving AMOC variability. The potential for Amazon runoff to impact salinity and consequently the AMOC was also identified by Latif et al. (2000), albeit caused by changes in precipitation rather than vegetation. We are aware that Amazon dieback is an artifact of the model that reflects regional trends in climate change, specifically precipitation, that vary between GCMs (Schaller et al. 2011). HadCM3 is therefore an extreme case, with other GCMs more consistently simulating a drier seasonal environment that remains forested (Malhi et al. 2009). However, the potential feedbacks of the land surface may have implications for predicting variability in climate models under increasing CO 2 , a factor that has been largely overlooked. Furthermore, it highlights a teleconnection between relatively localized terrestrial land surface cover and ocean variability that subsequently impacts climate across the North Atlantic region.

Summary and discussion
This study has used the HadCM3 coupled climate model run for millennial time scales to investigate the strength and variability of the AMOC at four quasi-equilibrium CO 2 concentrations; 350, 700, 1050, and 1400 ppm (experiments 1x, 2x, 3x, and 4x, respectively). The AMOC streamfunction and an assigned AMOC index (MOI) have been decomposed and reconstructed with multivariate singular spectrum analysis (MSSA) and analyzed with EOFs in order to remove noise and isolate variability. The key findings are the following: d AMOC strength is shown to decrease at higher concentrations of CO 2 with a reduction of 30.2% at 2x, 29.8% at 3x, and 34% at 4x. This is likely to be driven by a reduction in density in the Greenland-Iceland-Norwegian (GIN) Seas, an increase in stratification and so reduced convection that inhibits overturning.
d Only a small proportion of variability is explained by the first and second eigenmodes of the MSSAanalyzed AMOC streamfunctions (i.e., 11.19% at 1x). This is increased for the SSA of the MOI (i.e., 39.53% at 1x).
d Analysis indicates that low-frequency variability of the AMOC is also weakened at higher CO 2 . Significant frequencies are on the order of 100-125 yr for 1x, 2x, and 3x, which increases to 160 yr at 4x, with a consecutive decrease in power with increasing CO 2 .
d The decline in the extent of AMOC strength and variability has consequent impacts on climate during periods of maximum and minimum AMOC.
d AMOC variability is likely to be driven primarily by salinity-driven density changes in the GIN Seas, with anomalies advected into the region predominantly from the Arctic and Caribbean regions.
d A positive salinity anomaly acts to increase density, which reduces stratification and consequently increases convection; this then acts to increase downwelling and consequently AMOC strength, increasing SSTs and SATs and decreasing SLP, which enhances Ekman pumping; and this in turn acts as a positive feedback, further enhancing convection and AMOC strength.
d In HadCM3, salinity anomalies advected from the Arctic may be key in driving variability, with the Caribbean playing a secondary role.
d Higher concentrations of CO 2 are associated with a weakening in the amplitude and coherence of Arctic salinity, which in turn may be responsible for reducing the amplitude of AMOC variability. In contrast, Caribbean fluctuations remain consistent and may increase, indicating that this region may play a more predominant role at higher CO 2 .
d Arctic salinity anomalies may be generated by deep northward transport from the GIN Seas as indicated by a contrasting pattern of salinity anomalies seen at depth in the Beaufort Gyre. This seesaw, depthvarying motion of anomalies may act as a positive feedback to reduce variability at higher CO 2 (i.e., decreased variability results in weaker northward advection, weaker anomaly generation, and consequently weaker variability). This may be a response to the overall reduced strength of the AMOC; however, further study incorporating a salinity budget and/or tracer analysis is required to test this hypothesis.
d Tropical salinity anomalies may be formed by a shift in the position of the ITCZ at 1x as identified in the study of VW04. However, at higher CO 2 an increase in the residence time and/or change in surface runoff may play an important role in maintaining the amplitude of anomalies. The former may be a response to a weaker overall AMOC and an increase in stratification, which increases the residence time of anomalies in the evaporative subtropical gyre. The latter may be linked to a CO 2 -induced vegetation dieback in the Amazon FIG. 16. Lagged correlations of the MOI with MSSA decomposed Arctic salinity anomalies at a range of depths for (top)-(bottom) 1x, 2x, 3x, and 4x. Thick dashed lines highlight 95% confidence for all time series (see section 2). Salinity anomalies lead the MOI for negative lags region, which increases surface runoff as a proportion of precipitation. These drivers may contribute in maintaining the generation of tropical salinity anomalies that helps drive AMOC variability despite weaker Arctic fluctuations.
The modeled reduction in amplitude of variability that accompanies a weaker overall AMOC that we see in this study is in contradiction to theoretical analysis of timedependent systems, which predict an increase in variability as a system weakens and becomes more unstable (e.g., Tziperman 1997). This was first highlighted in the study of MacMartin et al. (2016) using the GFDL model and we show that this is also the case with HadCM3.
The results here highlight that AMOC variability is driven by a combination of mechanisms in HadCM3: a possible internal ocean mode in the Arctic region and a coupled atmosphere-ocean mode in the Caribbean region. Higher concentrations of CO 2 act to alter these mechanisms, weakening both the Arctic and Caribbean modes and possibly driving a shift to an atmosphere-terrestrial-ocean feedback in the tropics.
This is a single-model study and some of the mechanisms presented here and how they may be influenced by FIG. 17. Land surface cover and runoff in the Amazon region for (top)-(bottom) 1x, 2x, 3x, and 4x. (a) Land surface cover, (b) the proportion (%) of runoff relative to precipitation for each grid box (%), and (c) composite runoff anomaly maps (mm day 21 ) calculated during periods of AMax averaged over periods when the PC1 of the AMOC streamfunction is one standard deviation above the mean. Only anomalies that are considered 95% confident are shown (see section 2). CO 2 are likely to be highly dependent on the model used. Inaccuracies in the simulated overflow channels and issues pertaining to resolution may have implications for the AMOC. These subsequently lead to the modeled AMOC being too shallow compared to observations, which may have implications for the time scale and the mechanisms that drive variability. Despite this, it is important and interesting to understand how variability in models may respond to climate forcing. Applying the methods outlined in this paper would elucidate how variability in other models responds to climate change. Understanding model internal variability is crucial in order to differentiate impacts that are thought to be anthropogenically forced (i.e., a CO 2 -induced weakening) relative to those that are a ''natural'' response of the modeled system. This may help us improve our ability to more accurately predict the AMOC and thus model decadal-centennial climate.