High-latitude precipitation as a driver of multicentennial variability of the AMOC in a climate model of intermediate complexity

Centennial-scale variability of the Atlantic Meridional Overturning Circulation (AMOC) in the absence of external forcing has been identified in several climate models, but proposed mechanisms differ considerably. Therefore, better understanding of processes governing AMOC variability at these timescales is needed. Here, we analyze numerical simulations with PlaSim-LSG, an Earth System Model Intermediate Complexity (EMIC), which exhibits strong multicentennial oscillations of AMOC strength under constant pre-industrial boundary conditions. We identify a novel mechanism in which these oscillations are driven by salinity anomalies from the Arctic Ocean, which can be attributed to changes in high-latitude precipitation. We further corroborate our findings by conducting a set of millennial-length sensitivity experiments, and we interpret the mechanism by formulating a three-box model which qualitatively reproduces regular oscillations of the AMOC. While PlaSim-LSG lacks complexity compared to state-of-the-art models, our results reveal that precipitation minus evaporation (P–E) change in the Arctic is a physically plausible driver of centennial-scale AMOC variability. We discuss how this mechanism might be most relevant in climate states warmer than the present-day, raising questions about the state-dependence of multicentennial AMOC variability.


Introduction
Through its northward transport of heat and salt, the Atlantic Meridional Overturning Circulation (AMOC) plays an important role in governing the climate of the North Atlantic region. Therefore, there has been considerable interest in understanding the variability of the AMOC across timescales, from interannual to multidecadal (Buckley and Marshall 2016) and millennial scales (Lynch-Stieglitz 2017). For example, recently, several studies using comprehensive climate models found unforced millennial-scale AMOC oscillations under glacial boundary conditions (Vettoretti et al. 2022;Klockmann et al. 2020;Kuniyoshi et al. 2022;Romé et al. 2022), which resemble the Dansgaard-Oeschger events found in paleoclimate records.
However, at intermediate, (multi-)centennial timescales, the variability of the AMOC has been studied less extensively, and we will focus on these timescales in the remainder of this study. Although sea surface temperature (SST) proxy records from the North Atlantic region exhibit significant multicentennial variability during the Holocene (Askjaer et al. 2022), the length of the instrumental record, uncertainties and non-stationarity of the AMOC-SST relation (Lozier 2010;Tandon and Kushner 2015;Moffa-Sánchez et al. 2019) and the lack of circulation proxies at sufficient resolution (Lippold et al. 2019) prevent a full characterization of the AMOC at these timescales. Hence, climate models have often been invoked to examine AMOC variability on timescales beyond the instrumental record, which only dates back to 2004 (Cunningham et al. 2007).
In particular, millennial-length integrations with constant pre-industrial forcing allow us to assess internal variability on centennial timescales. Early studies using ocean general circulation models (OGCMs) (Mikolajewicz and Maier-Reimer 1990;Winton and Sarachik 1993) and simplified (1-D and 2-D) ocean models (e.g., Sévellec et al. 2006) suggested the existence of a "loop oscillation" (Winton and Sarachik 1993) in which a salinity anomaly would be advected within the entire Atlantic overturning cell on the characteristic timescales of the thermohaline circulation (that is, multicentennial). While some more complex coupled GCMs support the notion of an oceanic mode of multicentennial variability in the North Atlantic driven by interhemispheric salinity transport (Park and Latif 2008;Delworth and Zeng 2012), others have proposed different mechanisms involving atmospheric or sea ice feedbacks. Vellinga and Wu (2004) proposed that increased precipitation in the Intertropical Convergence Zone (ITCZ) may be a driver of subtropical salinity anomalies impacting AMOC strength on a centennial timescale, a mechanism which Menary et al. (2012) showed to be present in at least two different GCMs. More recently, strong multicentennial AMOC oscillations were discovered in several models of the Coupled Model Intercomparison Project Phase 6 (CMIP6). They were linked to the build-up and release of Arctic Ocean freshwater anomalies moderated by sea ice in both the IPSL-CM6-LR (Jiang et al. 2021) and EC-Earth3 models (Meccia et al. 2022). In contrast, Li and Yang (2022) used a box model to argue that no coupled atmosphere-ocean feedback is required to sustain multicentennial oscillations in CESM1.0 (a CMIP5 model), and that this mode of variability can instead be explained by an eddy-induced ocean mixing feedback.
Because state-of-the-art models are computationally expensive, they generally do not allow for a large number of sufficiently long runs for sensitivity tests of multicentennial AMOC variability. Many proposed mechanisms were thus derived from data analysis of single pre-industrial control simulations, with the exception of Jackson and Vellinga (2013) who also analyzed short (400 years) runs from a perturbed physics ensemble with eight members. To improve our understanding of the mechanisms behind centennialscale AMOC variability, it therefore seems beneficial to trade off reduced model complexity for computational speed using Earth System Models of Intermediate Complexity (EMICs; Claussen et al. 2002), which allow for millenniallength sensitivity runs. While most other EMICs exhibit no or insufficient internal variability (Petoukhov et al. 2005), the group of "simplified comprehensive models" (Claussen et al. 2002) is sufficiently complex to generate internal variability on centennial timescales (e.g., Friedrich et al. 2010;Severijns and Hazeleger 2010).
In this work, we show that one such simplified GCM, PlaSim-LSG, exhibits significant internally driven AMOC oscillations at multicentennial timescales under constant pre-industrial boundary conditions. Using the control simulation and an ensemble of millennial-length sensitivity experiments, we analyze the role of atmosphere-ocean freshwater feedbacks in the high latitudes in driving these AMOC oscillations. Finally, we discuss how our study can serve as a starting point for exploiting the entire model hierarchy for investigating multicentennial AMOC variability.

Model experiments
Coupled atmosphere-ocean simulations are performed using PlaSim-LSG, an EMIC which couples the Planet Simulator (PlaSim; Fraedrich et al. 2005), an atmospheric GCM with a spectral dynamical core and simplified physics (Lunkeit et al. 2011), to the Large-Scale Geostrophic Ocean model (LSG;Maier-Reimer et al. 1993). The use of a geostrophic model, in which the nonlinear terms of the Navier-Stokes equation are neglected, is motivated by the scale analysis of Hasselmann (1982) for climate variability in a coarseresolution ocean model. Since AMOC variability on timescales of decades and longer is largely attributed to changes in geostrophic circulation (Buckley and Marshall 2016), we believe that LSG is well-suited for studying AMOC variability on centennial timescales.
PlaSim-LSG has fully interactive components for the atmosphere, ocean, sea ice and the hydrological cycle, while ice sheets and vegetation are prescribed. The threedimensional atmosphere and ocean components are coupled through the surface fluxes of momentum, heat and freshwater (Lorenz 2006) without flux corrections. Coupling of heat fluxes between PlaSim and LSG is performed within a shared slab ocean with a constant depth of 50 m, which also acts as the uppermost layer of LSG. This slab ocean serves as an intermediary between PlaSim and LSG, filtering high-frequency noise and damping short-lived perturbations, which-together with time-implicit integration-allows for an ocean and coupling time step of 10 days compared to 45 min for the atmosphere. Sea ice is formed and melted thermodynamically in the slab ocean based on the zero-layer model of Semtner (1976). No transport of sea ice is considered. To stabilize the large-scale ocean circulation and to isolate the effect of P − E , runoff into the oceans is redistributed globally in the coupling step. The local surface freshwater flux Φ surf into the ocean is therefore given by where ⟨R⟩ is the global sum of runoff into the oceans divided by the total ocean surface, and ḋ sice is the time derivative of sea ice liquid water equivalent. (1) PlaSim is configured at T21 resolution (corresponding to about 5.6 • × 5.6 • on a latitude-longitude grid) with 10 vertical levels in the atmosphere. LSG uses an "E"-type semi-staggered grid (Arakawa and Lamb 1977) with an effective horizontal resolution of 3.5 • × 3.5 • and 22 vertical layers on z-coordinates with a spacing between 50 m in the upper ocean and 1000 m in the deep ocean. The main difference with respect to the original LSG version (Maier-Reimer et al. 1993) is the introduction of the Farrow and Stevens (1995) predictor-corrector scheme for advection (cf. Prange et al. 2003), which is less diffusive than the original upstream scheme. As a consequence, ocean vertical diffusivity A v can be controlled explicitly via the parametrization of Bryan and Lewis (1979): The large-scale characteristics of the AMOC in PlaSim-LSG strongly depend on the chosen values for this parametrization, especially on the vertical diffusivity in the upper ocean layers. In a preliminary study (Angeloni et al. 2020), we kept the bottom diffusivity A v (6000 m) as well as = 4.5 × 10 −3 m −1 and z * = 2500 m fixed, and identified different AMOC regimes for different values of A v (0) in PlaSim-LSG ( Fig. 1): For low upper ocean diffusivities (A v (0) ≲ 0.2 cm 2 s −1 ), the AMOC collapses. For intermediate values, the model exhibits a relatively constant AMOC strength of about 17-19 Sv (1 Sv = 10 6 m 3 s −1 ), before a state with multicentennial oscillations emerges for A v (0) = 0.8 cm 2 s −1 , corresponding to a * = 1.0479 cm 2 s −1 and a range = 0.1673 cm 2 s −1 . Finally, for even higher values of the upper ocean diffusivity, these oscillations disappear again and the Atlantic is in a state of very strong overturning (about 30 Sv).
In this study, we analyze a 3000-year simulation (after 1000-year spinup) of PlaSim-LSG which exhibits multicentennial AMOC oscillations (A v (0) = 0.8 cm 2 s −1 ) in more detail. We use constant pre-industrial boundary conditions (pCO 2 = 285 ppm) and will refer to this simulation as the control simulation. In addition, we perform a set of sensitivity experiments to isolate the influence of Arctic freshwater anomalies onto AMOC oscillations, in which the amplitude of Φ atm = P − E + ⟨R⟩ anomalies over the Arctic Ocean is scaled. To this end, we diagnose the monthly climatology Φ atm,i from the control simulation, compute anomalies Φ � atm,i with respect to this climatology during each coupling step, and multiply the anomalies by a factor c, such that the scaled freshwater flux is given by Φ c atm,i =Φ atm,i + c Φ � atm,i . Here, the index i ∈ {1 … 12} denotes the month to emphasize the difference from annual mean climatologies and anomalies defined later in this section. Because the scaling is only applied to the Arctic Ocean, its effect on the global freshwater balance is negligible and no correction is applied.
This is in line with Rahmstorf et al. (2005), who argued that compensating for freshwater forcing in a different ocean basin had negligible effects even when the added freshwater flux was more than an order of magnitude stronger than in our study. Each sensitivity experiment is initialized from year 2000 of the control simulation and is integrated for 2000 years.
While PlaSim-LSG reproduces large-scale patterns of the hydrological cycle reasonably well even at an atmospheric resolution of T21 (not shown), there are several shortcomings in its climatology, as can be expected from a simplified GCM with very coarse resolution. During the spinup, there is a strong drift in the Southern Ocean similarly to earlier coupled versions of LSG (e.g., von Storch et al. 1997). Following this drift, Southern Ocean surface temperatures show a strong warm bias of about 10 K and virtually all Antarctic sea ice disappears. Nevertheless, observed zonal mean temperatures are reproduced well in the northern hemisphere, which is the focus of this study. Mean Arctic sea ice concentrations in the PlaSim-LSG control simulation are lower than in observations of the late 19th century (Fig. S1a, b; Rayner et al. 2003) and in the piControl simulations of most CMIP6 models (Fig. S2a), especially in boreal summer. The simulated mean Arctic Ocean salinity of 33.8 psu (2). The diffusivity at the bottom of the ocean A v (6000 m) is kept at 1.3 cm 2 s −1 throughout the ensemble, while the parameters a * and a range are adjusted to match the given diffusivity at the surface A v (0) . In all simulations, = 4.5 × 10 −3 m −1 and z * = 2500 m is also significantly lower than in observations (Fig. S1c, d;Zweng et al. 2018). However, the AMOC mean of 19.5 Sv (Sect. 3.1) agrees with observations (Cunningham et al. 2007) and lies well within the CMIP6 range (Bellomo et al. 2021).

Diagnostics
We define AMOC strength as the maximum of the meridional overturning streamfunction in the Atlantic between 40 • N and 60 • N, as the maximum of the overturning cell is located in this range (Fig. S3). Composites for the strong and weak AMOC phases are obtained by averaging over 21-year intervals around the maxima and minima of the AMOC time series. For the increasing and decreasing phases, composites are centered around the midpoints between these minima and maxima. To define the extrema, a 100-year running mean (black line in Fig. 2) is applied to the AMOC time series solely for the purpose of peak detection. All subsequent diagnostics, including composites and lagged regressions, are computed on unfiltered annual mean time series.
To decompose the competing effects of salinity S and potential temperature on density anomalies, we perform a Taylor expansion of the equation of state (cf. Vellinga and Wu 2004). is expanded at the local climatological mean (s(x, y, z),̄(x, y, z)) , from which salinity and potential temperature deviate by a small s ′ and ′ , respectively: Here and in all of the following analysis, all results refer to annual mean quantities (x) if not stated otherwise. The term "climatological mean" x refers to the mean of x over the entire simulation, while "anomalies" x ′ are the deviation of annual means from this climatology. Following (3), we can decompose density anomalies (to a very good approximation) into a salinity contribution ′ s and a temperature contribution ′ .
When examining the salinity budget of the Arctic Ocean, it is useful to express transport through the liquid freshwater flux (e.g., Lique et al. 2009) where u is the velocity across a section of area dA , whose normal vector is defined to point into the Arctic Ocean. The integral is taken over one horizontal and the vertical dimension, either over the full ocean depth or over the upper 300 m. Here, we choose s 0 to match the simulated Arctic Ocean mean salinity of 33.8 psu. While some authors (e.g., Schauer and Losch 2019) have criticized the use of ocean freshwater fluxes because of their nonlinear dependence on the "nonunique" s 0 , our choice of s 0 is physically motivated and makes Φ liq readily interpretable: In the climatological mean, the net export of liquid freshwater from the Arctic Ocean through its gateways is approximately balanced by the positive freshwater flux at the surface. Note that in the LSG setup used here ( Similarly to the equation of state, u and s in the integrand of (4) can be expanded into a mean and an anomaly term: obtained from the median-smoothed spectrum with a smoothing window Δf smooth = 0.05 year −1 following Appendix A2 of Mann and Lees (1996) Integrating (5), its second term is interpreted as the contribution to the freshwater flux due to advection of salinity anomalies by the mean current and its third term as the contribution to the freshwater flux due to transport of mean salinity by current anomalies The residual term is not small everywhere, but it has a weak dependence on AMOC strength ( Fig. S4), such that it is neglected in the analysis below. When computing (lagged) regression coefficients, we test their significance using the "random phasing" method of Ebisuzaki (1997) to take into account the strong autocorrelation of many quantities. To this end, we construct 1000 surrogate time series of the regressor, which have identical Fourier spectra but differ in their randomly chosen phases for each frequency. After repeating the regression for each of these surrogate time series, we consider regression coefficients significant at the (twotailed) 95% confidence level if they are larger than the 97.5th or smaller than the 2.5th percentile of the resulting distribution.

Life cycle of salinity and circulation anomalies
The 3000-year time series of AMOC strength of the control simulation is shown in Fig. 2a. AMOC strength has a mean of 19.5 Sv and varies on multicentennial timescales with a peak-to-peak amplitude of about 3-4 Sv. AMOC strength at the maximum of the overturning cell is in phase with the AMOC at 26.5 • N, although oscillations are weaker there. Multicentennial variability of the AMOC is characterized by regular, sinusoidal oscillations with similar amplitude throughout the control simulation. Their mean period is about 270 years as determined from the first maximum of the autocorrelation function. In the power spectrum (Fig. 2b), the oscillations are represented by a remarkably high peak in the range between 200 and 400 years, which exceeds the 99% significance threshold by nearly two orders of magnitude. It is the only peak in the spectrum to exceed the threshold, indicating the absence of spectrally consistent unforced variability on shorter timescales.

Salinity-driven density anomalies
The AMOC variability in our control simulation is accompanied by strong near-surface density changes in the midand high latitudes. Figure 3 shows composites of density anomalies for the four phases of an AMOC oscillation, each obtained from a 21-year interval per oscillatory cycle. The strongest positive density anomalies occur in the Labrador Sea at the AMOC maximum, as well as in the Norwegian Sea and the Barents/Kara (BK) sea during the increasing AMOC phase. The two former regions are close to the main areas of convective activity in the North Atlantic (Fig. 4d). Compared to these regions, in the South Atlantic and in the Southern Ocean typical density (Fig. S5) and salinity anomalies (Fig. S6) are about one order of magnitude smaller. Therefore, it appears unlikely that these regions in the southern hemisphere play a significant role in driving variability in the North Atlantic, and we focus on processes in the northern hemisphere in the following.
We decompose density anomalies into a salinity contribution ′ s and a temperature contribution ′ using (3); then, we average them regionally and over the upper 300 m. Lagregression analysis (Fig. 4a-c) shows that the salinity contribution is almost in phase with the total density anomaly in the regions with the strongest density changes. The temperature contribution is weaker and has the opposite phase. Hence, in all three regions analyzed in Fig. 4, near-surface density changes are driven by salinity rather than temperature changes. In the Arctic Ocean, density changes are completely governed by salinity, while the temperature contribution can be neglected, as can be expected at sea surface temperatures near freezing (Aagaard and Carmack 1989). Therefore, we focus on salinity changes as a driver of nearsurface density in the following.
The maximum of ′ s leads the AMOC by about a quarter of a period (73 years) in the BK sea, while it occurs slightly after the AMOC maximum (lag + 2 years) in the Labrador Sea. In the Norwegian Sea, the maximum occurs at lag −82 years. While we expect salinity in the Labrador Sea to be roughly in phase with the AMOC because an enhanced AMOC means that more salty water is transported here from the lower latitudes, the phase lags in the Nordic Seas and the Barents Sea are less straightforward to interpret. In particular, salt transport is determined not only by local salinity but also by changes in circulation, which are significant during PlaSim-LSG oscillations. This interaction will be investigated in the following section.

Decomposition of freshwater transport
To disentangle the effects of salinity and circulation changes on freshwater export from the Arctic Ocean, we decompose anomalies of the liquid freshwater flux at the Fram strait and the Barents section into advection of salinity anomalies by the mean current Φ � s (6) and transport of mean salinity by current anomalies Φ � u (7). We integrate only over the top 300 m because of the shallow depth of the Barents sea and because circulation changes in the Fram strait are mostly barotropic (Fig. S7).
Lag regression of Φ � s and Φ � u onto AMOC strength is shown in Fig. 5. Freshwater transport change in the Fram strait is almost completely determined by Φ � u , while Φ � s and Φ � u have a similar amplitude in the Barents section. Φ � s in the Barents section is in phase with the salinity anomaly in the entire BK sea (cf. Figs. 4a and 5b). For both sections combined, Φ � s leads the AMOC by 62 years and Φ � u lags the AMOC by 23 years, while the total freshwater flux is in phase with the AMOC. Since the mean freshwater flux is negative ( −30 mSv for the top 300 m of Fram strait and Barents section combined), this means that freshwater export from the Arctic is at its minimum during an AMOC maximum. The phase of Φ � u can be explained to a first order by geostrophic flow, with the gradient of sea surface height (SSH; Fig. S8), which is driven by freshwater anomalies in the same way as near-surface salinity, determining circulation anomalies in the upper ocean. The SSH and salinity gradients are especially large in the eastern Fram strait, where circulation anomalies indeed appear to run parallel to density anomaly isolines (Fig. 3).
From Fig. 5c, it may seem that the contribution of salinity and current anomalies cancel out at about lag −70 years. However, this is only true locally. South of the Fram strait and Barents section, salinity anomalies advected by the mean current and current anomalies transporting mean freshwater are exported to different locations. Hence, the two  (5) to each gridpoint. Figure 6 visualizes salinity anomalies overlaid by the mean current, while Fig. 7 shows mean freshwater S 0 − S overlaid by current anomalies. Approximate local values for Φ � s and Φ � u can be obtained by multiplying the two fields in each plot. Lag −70 years is within the increasing AMOC phase. Here, current anomalies drive an enhanced transport of freshwater southward to the Norwegian Sea (Fig. 7b), while the positive salinity anomalies in Fig. 6b are advected to the Denmark strait along the East Greenland Current, causing an export of salt from the high latitudes (i.e., an input of freshwater into the Arctic region). The different export locations of Φ � s and Φ � u are crucial for the salinity cycle described in Sect. 3.1.1, since no salinity reinforcement to the AMOC could be provided in the Labrador Sea if Φ � u and Φ � s were balanced out everywhere.
This pathway along the mean East Greenland Current is supported by the presence of a pronounced salinity-induced density anomaly in the Denmark strait at lag −30 years (Fig. S9). This salinity maximum is located below the surface layers because the mean current has a negative vertical component of the order of 10 −6 m s −1 in the northern Greenland Sea, where convection is practically absent. This means that it is four orders of magnitude smaller than the horizontal current, equivalent to a downwelling of about 200 m along a distance of 2000 km between the Barents section and the Denmark strait. Finally, Figs. 6b and 7b show that the salinity anomaly off the coast of Norway does not significantly contribute to the Arctic-North Atlantic salinity cycle, since

P-E changes as driver for PlaSim-LSG oscillations
We demonstrated that near-surface salinity anomalies in the BK amplify AMOC oscillations, which would otherwise taper off. We now investigate the role of net surface freshwater flux in driving these salinity anomalies.

Role of Arctic P-E in the control simulation
First, we examine the contributions of P − E + ⟨R⟩ and sea ice thickness changes to salinity in the BK sea in the control simulation. To this end, we diagnose salinity tendencies ̇s related to P − E + ⟨R⟩ and to changes in sea ice volume at each timestep online within LSG. In the model, the surface freshwater flux only affects salinity in the uppermost layer, but subsequently interacts with deeper layers through advection, diffusion or convection. Figure 8 shows integrated annual mean anomalies of the diagnosed salinity tendencies (total tendencies, those related to P − E + ⟨R⟩ , and those related to sea ice) over an interval Δt = 70 years . Integration is performed to isolate the contribution to centennial-scale variability, which is clearly visible in the integrated time series. Hence, Δs can be interpreted as "low-frequency salinity changes".
In the BK sea, Δs P−E+⟨R⟩ is very closely related to Δs total (r = 0.83, p = 0.001), while sea ice-induced low-frequency salinity changes are weaker and tend to oppose the total salinity change (r = −0.42, p = 0.08). The residual between the total and the P − E + ⟨R⟩-related salinity changes is not significantly correlated with Δs total (r = 0.19, p = 0.34), making P − E + ⟨R⟩ a plausible driver for salinity changes in the BK sea.
Assessing P and E separately, the magnitude of precipitation changes is larger than that of evaporation changes over the entire Arctic Ocean, but especially in the BK sea (Fig. 9a, b). This is in contrast to the convective regions of the North Atlantic, where P − E changes are evaporationdriven (Fig. 9c, d). Because P − E is determined by the convergence of the moisture flux Q (P − E = −∇ ⋅ Q) on annual and longer timescales (Peixoto and Oort 1992), the freshwater anomaly over the Arctic Ocean can be directly related to an anomaly in moisture transport towards the Arctic. During an AMOC maximum, this moisture transport strengthens  especially over the Irminger and Greenland Seas, fueled by positive evaporation anomalies over the central North Atlantic (Fig. S10). These evaporation anomalies are in turn associated with North Atlantic SST anomalies of up to 2 K at the AMOC maximum. In addition, we computed a simple moisture budget over the region north of 75 • N following the method of Schär et al. (1999) to evaluate the importance of evaporation changes within the Arctic on modulating Arctic precipitation. The budget reveals that, in the annual mean, 16% of Arctic precipitation is sourced from evaporation over the same region, and this value does not differ significantly between different phases of the AMOC (not shown). Hence, a significant influence of sea ice-driven evaporation anomalies on more localized precipitation anomalies (such as in the BK sea) appears to be unlikely.

Sensitivity to Arctic P-E changes
While the data analysis suggests that Arctic P − E changes are the key atmospheric feedback in driving multicentennial AMOC variability in PlaSim-LSG, we seek for a more rigorous way to test this hypothesis by performing a set of sensitivity experiments. To this end, we varied the scaling factor c for monthly P − E anomalies in the Arctic Ocean as described in Sect. 2.1 between 0 and 4, with c = 1 corresponding to the original simulation analyzed above.
The 2000-year AMOC time series of these sensitivity experiments are shown in Fig. 10. For small scaling factors, i.e., approaching a relaxation towards the P − E climatology, the AMOC slowly transitions to a very strong (around 30 Sv) state without multicentennial oscillations. For large scaling factors (here: c = 4 ), AMOC strength rapidly decreases and the overturning cell collapses north of about 45 • N within 150 years. The stability of these two regimes can be assessed by resetting c to 1 after the system has approached its new equilibrium following the initial perturbation. While the strong AMOC state appears to be unstable and the AMOC returns to the attractor of the control simulation, the AMOC does not recover from its collapsed state within 2000 years (Fig. S11). This behavior is a strong indicator of bistability.
Multicentennial oscillations occur for a wide range of intermediate scaling factors, between 0.5 and 2 in the set of experiments performed here. Within this range, both the period and the amplitude decrease significantly with an increasing amplitude of Arctic freshwater forcing (Fig. S12). While the change in amplitude is probably related to slow feedback processes, longer periods for smaller c are consistent with the mechanism of P − E-driven salinity anomalies: Fig. 9 Lag regression of precipitation, evaporation and P − E + ⟨R⟩ onto AMOC strength for different ocean regions. Here, we use the sign convention that downward fluxes are positive (precipitation is positive, evaporation is negative). Thick lines indicate significant regression coefficients at the 95% confidence level. Negative (positive) lag means that the freshwater flux time series leads (lags) the AMOC. Note that the y-axes have different scales it takes longer for Arctic salinity anomalies to build up under weaker low-frequency freshwater forcing. The second timescale at play is illustrated by the lag correlation between salinity and AMOC strength for different scaling factors (Fig. 11). Near-surface salinity in the BK sea and sub-surface salinity in the Denmark strait robustly lead the AMOC by 65-85 years and 10-50 years, respectively, with no apparent dependence on c. This supports the assumption that the reinforcement for the AMOC is provided by salt advection from the Arctic Ocean via the Nordic seas and that the corresponding timescale is determined by the mean circulation in this region. The maximum of BK sea near-surface salinity regression onto AMOC strength is also remarkably similar across different ensemble members for 0.5 ≤ c ≤ 1.5 , indicating good predictive power of BK sea salinity for AMOC strength.

Interpretation using a three-box model
According to our analysis of the control simulation and to the sensitivity experiments, the central element of the mechanism for AMOC oscillations in PlaSim-LSG is the feedback between high-latitude P − E , Arctic Ocean salinity and AMOC strength. To test whether this mechanism is sufficient to explain the existence of oscillations, we illustrate the main processes at play with a simple three-box model based on the canonical two-box model for the thermohaline circulation by Stommel (1961). In the Stommel model, the two boxes symbolize well-mixed subtropical and North Atlantic ocean basins, with a temperature difference T = T s − T NA and a salinity difference S = S s − S NA ( T, S > 0 , i.e., the subtropical Atlantic is warmer and more saline than the North Atlantic). The density-driven flow Ψ = T − S between the two boxes is identified as AMOC strength, which is positive for a thermally driven AMOC and negative for a salinity-driven AMOC (see e.g. Chapter 3 of Dijkstra (2005) for a thorough description of the original Stommel model and its solutions).
For a minimal extension of Stommel's model to represent the mechanism described here, we introduce a third box representing the Arctic Ocean, which has a salinity anomaly S a . The two drivers of this salinity anomaly are Arctic P − E anomalies and the transport of Arctic salinity anomalies into the North Atlantic; the latter corresponding to the Φ � s term in the PlaSim-LSG analysis. Following our reasoning from above, P − E anomalies are in turn driven by the integral of moisture flux anomalies at the boundary of the Arctic region. For simplicity, we assume that these moisture flux anomalies, and therefore Arctic P − E , follow a linearized thermodynamic scaling (in the sense of Held and Soden 2006) proportional to North Atlantic temperature anomalies Here, T NA denotes the time mean of T NA , such that T s −T NA can be expressed in terms of a single constant T c . To capture the difference in timescales between the Arctic and the North Atlantic, we formulate the three-box model as a fast-slow system: where is the atmospheric coupling strength of the Arctic hydrological cycle to the North Atlantic and < 1 is the timescale of salinity changes in the Arctic, i.e., the AMOC is the fast and the Arctic Ocean the slow component. We also carry over the parameters i of the original Stommel (1961) model, where 1 can be interpreted as the strength of the temperature forcing at the surface, 2 as the strength of the salinity forcing at the surface, and 3 as the ratio between the timescales for salinity and temperature restoration (Dijkstra 2005). Assuming T s and therefore T c to be constant is in good agreement with PlaSim-LSG, where subtropical ocean temperatures only change by some tenths of a degree during one oscillatory cycle. Several sets of parameters can be found for which oscillations of S a and Ψ occur in the three-box model. To be consistent with the present-day AMOC, which is thermally driven (Rahmstorf 2002), we impose that Ψ must be positive in the original Stommel model, i.e., for → 0 . In addition, is chosen such that the Arctic Ocean timescale is larger, but still on the same order of magnitude as the AMOC timescale. For smaller , we would obtain a relaxation oscillator more reminiscent of millennial variability in a glacial climate (Crucifix 2012), while the oscillations in PlaSim-LSG are practically symmetric with respect to their ascending and descending phases. In Fig. 12a, we show trajectories for different values of using one set of parameters that fulfills these criteria ( 1 = 2, 2 = 0.6, 3 = 0.3, = 0.25, T c = 1.8). In this configuration, the box model exhibits three AMOC regimes (Fig. 12b): for = 0 , we recover the original Stommel model with a non-oscillating, thermally driven AMOC. As increases, the globally attractive fixed point gradually shifts towards smaller values of Ψ before the system undergoes a Hopf bifurcation and oscillations occur. For sufficiently high values of , another stable fixed point appears for which Ψ < 0 , corresponding to a state without deep-water formation in the North Atlantic. The limit cycle and this stable fixed point on the lower branch coexist for ≳ 1.8 , but the basin of attraction of the limit cycle becomes smaller with increasing , leading to different attractors for = 2 and = 2.5 in Fig. 12, which are both in the bistable regime (Fig. 13a). For = 2 , the fixed point and the limit cycle, which are reached from different initial conditions, are depicted in Fig. 13b.
Our simple three-box model demonstrates that the interplay of AMOC strength, Arctic Ocean salinity and highlatitude precipitation can in theory be sufficient to explain oscillatory behavior of the AMOC. In addition, it captures two main features of PlaSim-LSG oscillations. First, salinity in the Arctic Ocean leads T − S (Fig. 12c). Second, for a given initial condition, oscillations can only be maintained if the amplitude of Arctic P − E changes is not too small or too large. The collapsed state and the oscillating state coexist in the box model for a wide range of parameters, analogously to the bistable behavior of the AMOC in PlaSim-LSG. To explain other features like the change in amplitude and periodicity in the sensitivity experiments, more complex models are likely needed.

Discussion and conclusions
EMICs are an attractive tool for studying centennial-scale AMOC variability, as trading off model complexity for computational cost allows to probe physical mechanisms thoroughly. In this study, we have shown that regular multicentennial AMOC oscillations occur in one such EMIC, PlaSim-LSG (Fig. 2). Combining analysis of the control simulation and sensitivity experiments, we identified lowfrequency variations in high-latitude P − E as the main atmospheric feedback driving these oscillations. P − E variations over the Arctic Ocean can be linked to changes in moisture transport to the Arctic from lower latitudes (Fig. S10). All of high-latitude precipitation, moisture transport to the Arctic and evaporation in the North Atlantic are lowest during a weak AMOC phase (Sect. 3.2.1), when a cold anomaly in near-surface temperatures persists across the northern mid-and high latitudes. The ensuing negative Arctic P − E anomaly leads to the build-up of a positive salinity anomaly in the Arctic Ocean, particularly in the BK sea (Fig. 8). This salinity anomaly is transported by the mean current to the Greenland Sea and reaches the North Atlantic within approximately 70 years (Fig. 6). Here, the salinity anomaly provides the reinforcement that strengthens the AMOC. When the AMOC transitions to a strong phase, atmospheric temperatures rise and the opposite phase of the cycle starts. We proposed a simple Stommel-type three-box model to demonstrate that this mechanism may explain regular oscillations of the AMOC in a physically plausible way. The results of our sensitivity experiments with PlaSim-LSG underlined the robustness of the mechanism at play and unveiled additional features like the characteristic timescales and bistability.
To our knowledge, this is the first study in which highlatitude precipitation is identified as a potential atmospheric feedback for multicentennial AMOC oscillations. Nevertheless, it shares many elements with previously proposed mechanisms. In particular, freshwater anomalies in the Arctic Ocean have also been identified as the central driver of multicentennial AMOC oscillations in IPSL-CM6A-LR (Jiang et al. 2021) and EC-Earth3 (Meccia et al. 2022), although freshwater anomalies there are driven by changes in sea ice in contrast to P − E in PlaSim-LSG. It is likely that these different drivers can be attributed to the different background climate state in PlaSim-LSG compared to these state-of-the-art models. In particular, sea ice concentration in the PlaSim-LSG control simulation is much lower than in CMIP6 piControl simulations and more akin to that of the last interglacial ( Fig. S2b; Otto-Bliesner et al. 2021), when global mean temperatures were about 2 • C higher than during the preindustrial period (Turney and Jones 2010). This makes it unlikely that the mechanism proposed here plays a significant role in the preindustrial or present-day climate. However, it highlights a possible mechanism for maintaining multicentennial AMOC variability in warmer climate states with a lower mean and variability of sea ice, in which the relative importance of moisture transport variations could become more important.
A precipitation-salinity-AMOC feedback similar to the one proposed here has previously been discussed by Vellinga and Wu (2004), but in the subtropical Atlantic. In their model, the subtropical precipitation anomaly signal due to an ITCZ shift dominated over a precipitation signal in the Nordic Seas (their Fig. 6c). In contrast, the absence of clear salinity anomalies linked to the ITCZ in PlaSim-LSG might be explained by differences in resolution or convective parametrization. Overall, a mechanism involving advection of salinity anomalies from the subtropical Atlantic, the South Fig. 13 Bistability in the three-box model: a sign of the real part of the eigenvalues of the Jacobian at the fixed points of (9)-(11) as a function of and T c . The bistable range is shown in dark blue, where one stable fixed point ( − − − ) and one unstable fixed point with complex conjugate eigenvalues ( − + + ) coexist. b Ensemble of 1000 trajectories starting from different initial conditions for T c = 1.8 and = 2 . Red dots mark the state of each solution at t = 400 . All other parameters are as in Fig. 12a Atlantic or even the Southern Ocean seems very unlikely in PlaSim-LSG due to the comparatively small amplitude of salinity and temperature changes in these regions (Fig. S6). While the absence of Antarctic sea ice, which is crucial for sustaining low-frequency oscillations in the Southern Ocean (Park and Latif 2008), in our simulation may over-emphasize the role of the northern hemisphere, our results add to a growing body of literature (e.g., Vellinga and Wu 2004;Jiang et al. 2021;Waldman et al. 2021;Li and Yang 2022;Meccia et al. 2022) which demonstrates how centennialscale AMOC variability can be driven without a significant contribution from the southern hemisphere, even in models with a realistic Antarctic sea ice climatology. This is consistent with the recent results of Askjaer et al. (2022), who showed that multicentennial surface temperature variability in both proxy records and transient climate model simulations of the Holocene is most pronounced in the northern hemisphere high latitudes.
Aside from the sea ice climatology, the main limitations of PlaSim-LSG are its low resolution which affects the North Atlantic storm track (Dong and Valdes 2000), and more importantly that local runoff into the oceans and sea ice dynamics are not (adequately) represented. In our simulation, runoff into the Arctic Ocean tends to be in phase with high-latitude precipitation such that we would expect it to amplify the P − E anomalies, strengthening the atmosphere-ocean feedback outlined above. While the potential effect of sea ice dynamics is harder to gauge, we would not expect it to alter the mechanism proposed here significantly, since it can only affect salinity anomalies in the Arctic Ocean indirectly by controlling the availability of sea ice (Meccia et al. 2022). However, we showed that the sea ice contribution to these salinity anomalies in PlaSim-LSG is less important than the P − E contribution (Fig. 8). Finally, we note that-similar to CMIP6 models but in contrast to some EMICs-PlaSim-LSG does not consider any coupling to ice sheets, whose freshwater discharge has been suggested to amplify multicentennial climate variability (Bakker et al. 2017).
Since centennial-scale AMOC oscillations have previously been reported in various versions and setups of LSG (Mikolajewicz and Maier-Reimer 1990;Pierce et al. 1995;Timmermann et al. 1998;Hertwig et al. 2015) involving different mechanisms, we strongly suspect that LSG has features which favor such oscillations. In particular, LSG is known to be highly diffusive (Maier-Reimer et al. 1993), even though the original upstream advection scheme has been replaced in the current version, and our control simulation used an even higher value for upper-ocean vertical diffusivity than the original Bryan and Lewis (1979) scheme. While other studies have pointed out the importance of the oceanic mixing parametrization for unforced AMOC oscillations (Peltier and Vettoretti 2014) and for AMOC hysteresis (Prange et al. 2003), we demonstrated how strongly the upper ocean vertical diffusivity can control not only the mean state, but also low-frequency variability of the AMOC (Fig. 1). This highlights the need to investigate the role of (vertical) mixing on multicentennial AMOC variability further in more complex models. In addition, this parametrization dependence of AMOC variability serves as a reminder to be cautious when relating model-specific phenomena to the real world.
In conclusion, our study has at least two implications for the study of multicentennial AMOC variability in state-ofthe-art models. First, the parallels between the mechanisms proposed by Jiang et al. (2021) and Meccia et al. (2022) and the one described here make us confident that PlaSim-LSG can serve as a testbed for advancing the understanding of multicentennial AMOC variability in CMIP6 models. For example, PlaSim-LSG could be used to design targeted sensitivity experiments for computationally more expensive models, while its apparent bistability could serve as a starting point to explore the interplay between AMOC stability (Weijer et al. 2019) and the existence of an oscillating AMOC state. Second, we provided evidence that highlatitude P − E can be a plausible driver of multicentennial AMOC oscillations. It appears that a warm background climate state in which the Arctic has significantly less sea ice than in the preindustrial climate would be required for such a P − E-driven mechanism to maintain AMOC variability, since the sea ice contribution to low-frequency freshwater flux variations dominates over the P − E contribution in models with a more realistic preindustrial sea ice climatology (e.g., in Jiang et al. 2021). For example, the mechanism proposed here could have acted during the last interglacial, but it is also a candidate to maintain low-frequency AMOC oscillations once sea ice variability decreases under global warming. Investigating this state-dependence of multicentennial AMOC variability, which has largely been unexplored so far, is an intriguing avenue for future work with models of different complexity.