AMOC Variability and Watermass Transformations in the AWI Climate Model

Using the depth (z) and density (ϱ) frameworks, we analyze local contributions to AMOC variability in a 900‐year simulation with the AWI climate model. Both frameworks reveal a consistent interdecadal variability; however, the correlation between their maxima deteriorates on year‐to‐year scales. We demonstrate the utility of analyzing the spatial patterns of sinking and diapycnal transformations through depth levels and isopycnals. The success of this analysis relies on the spatial binning of these maps which is especially crucial for the maps of vertical velocities which appear to be too noisy in the main regions of upwelling and downwelling because of stepwise bottom topography. Furthermore, we show that the AMOC responds to fast (annual or faster) fluctuations in atmospheric forcing associated with the NAO. This response is more obvious in the ϱ than in the z framework. In contrast, the link between AMOC and deep water production south of Greenland is found for slower fluctuations and is consistent between the frameworks.

underlie the AMOC variability. Simulations from ocean models driven by prescribed atmospheric forcing are often used to hindcast the AMOC for intercomparing model performance (e.g., Danabasoglu et al., 2016) or to understand the related dynamics such as the role of ocean eddies enhanced by increases in model resolution (see, e.g., Hirschi et al., 2020). However, standalone ocean models miss atmosphere-ocean feedbacks which can modulate modes of climate variability such as the Atlantic multidecadal variability (AMV, e.g., Oelsmann et al., 2020). That is why AMOC variability must be studied also in coupled atmosphere-ocean models, with transient (historical and future scenario) forcing (including solar, greenhouse-gas, and aerosol forcing) or with constant (e.g., preindustrial) forcing. The latter has the advantage that models can be run over long temporal scales, which increases the significance of analyses when it comes to understanding internal variability. Furthermore, analyses of AMOC variability in coupled models can shed light on important interactions in the climate system. Along these lines, Danabasoglu et al. (2012) analyzed the AMOC variability in CCSM4. The authors showed that various AMOC index time series tend to lead sea surface temperature (SST) changes in high latitudes and thus confirmed the driving role of the AMOC. Similar results were obtained by Frankignoul et al. (2013) based on CCSM3. Using 10 climate models, Roberts et al. (2013) identified a link between the trend in the AMOC at 26.5°N and subsurface density anomalies in the Subpolar Gyre while finding a large spread in patterns and magnitudes of SST response. Ortega et al. (2012) analyzed two IPCC scenarios and the last-millennium run in ECHO-G. They found that AMOC variability at different frequencies is controlled by different processes: high frequencies are associated with local changes in the Ekman transport caused by multiple modes of atmospheric variability, whereas low frequencies are linked with deep water production south of Greenland. Xu et al. (2019) analyzed 44 CMIP5 climate simulations and 18 standalone ocean simulations, finding no robust link between the AMOC and the North Atlantic Oscillation (NAO) in CMIP5, in contrast to a stronger link in standalone ocean runs. However, they found a strong link between the deep water production in the western subpolar North Atlantic and the AMOC also in the CMIP5 simulations. Menary et al. (2013) and Menary, Robson, et al. (2020) described the mechanisms of aerosol-forced AMOC variability in the Hadley Center Global Environment Model 2. Using the same model, the recent work by Menary, Jackson, & Lozier et al. (2020) reconciled the modeled and observational OSNAP data and demonstrated that the Labrador Sea may not be the origin of AMOC variability although it correlates with the densities there. Even though tremendous progress in understanding AMOC variability has been made over the past years, model results appear to be dependent on model uncertainties and resolution (see, e.g., Katsman et al., 2018;Menary & Hermanson, 2018;Menary, Jackson, & Lozier, 2020;Reintges et al., 2017;Xu et al., 2019).
There are two caveats in existing studies inquiring into the AMOC variability in the climate system. First, many studies still consider AMOC as a function of depth and latitude (z-AMOC), whereas AMOC as a function of density and latitude (ϱ-AMOC) might be a more appropriate diagnostic (see, e.g., Kwon & Frankignoul, 2014;Johnson et al., 2019;Zhang, 2010). Indeed, ϱ-AMOC is directly connected to surface watermass transformations, which are simply its constituent part (e.g., Walin 1982;Xu et al., 2018). The other part of ϱ-AMOC comes from internal transformations (see, e.g., Xu et al., 2018). The latter are mainly due to horizontal and vertical mixing, and partly due to the nonlinearity of the equation of state. In this respect, ϱ-AMOC explicitly reveals the roles of surface and internal transformations in sustaining AMOC variability, as discussed in Grist et al. (2009). Zhang (2010 and Kwon and Frankignoul (2014) explicitly mention the need in studying ϱ-AMOC for the understanding of the processes in subpolar latitudes. However, AMOC computation in density space is less straightforward and requires either high-frequency model output or online computations, which is rarely done. Hence, it is often compromised by using low-frequency output which varies from several days (see, e.g., Megan, 2018) to one month (see, e.g., Kwon & Frankignoul, 2014). This may bias the diagnostics when mesoscale processes come into play. The second caveat is related to the fact that the AMOC represents the net effect from processes happening at different locations. While numerous studies try to correlate variability in physical processes with that of the AMOC, the direct analysis of sinking through an appropriate depth level, or diapycnal transformation through the isopycnal that corresponds to the center of the main AMOC cell, is much more revealing (e.g., Xu et al., 2018). Such an analysis shows which locations are the main contributors to AMOC variability, and helps to judge the physics of processes that are involved.

10.1029/2021MS002582
3 of 17 The main goal of this paper is to understand local contributions to AMOC variability simulated in a new climate model, AWI-CM, by comparing the z-and density frameworks. Its secondary goal is to draw attention to the utility of local analyses of sinking or water mass transformations that constitute the AMOC when integrated over constant depth or density surfaces respectively. To this end we analyze the output of a preindustrial AWI-CM run over a period of 900 years. The water mass transformations are diagnosed online during the model simulation . We analyze not only the z-AMOC and ϱ-AMOC stream functions and their index time series, but also the spatial patterns of sinking and density transformations through depth levels and isopycnals.
The paper is organized as follows: Section 2 describes the model simulation. Section 3 introduces the density and depth frameworks for AMOC computation. Section 4 addresses the variability of both AMOCs and the associated climate patterns and mechanisms. The last two sections present the discussion and conclusions.

Model Simulation
A preindustrial control simulation (PIC) was generated with the AWI climate model (AWI-CM, Rackow et al., 2016;Sidorenko et al., 2015Sidorenko et al., , 2019 which is built upon the Finite-volumE Sea ice-Ocean Model (FES-OM 2.0; Danilov et al., 2004Danilov et al., , 2015Danilov et al., , 2017Scholz et al., 2019;Timmermann et al., 2009;Wang et al., 2008Wang et al., , 2014 and is available with several atmospheric components in climate setups. Here FESOM was run in combination with OpenIFS (cycle 43R3, Buizza et al., 2017), the atmosphere model developed at the European Center for Medium-Range Weather Forecasts (ECMWF). In terms of global-mean surface heat and freshwater budgets and ocean hydrography, this new climate configuration performs about as well as the configuration described in Sidorenko et al. (2019); details will be provided in a separate paper. AWI-CM was run under preindustrial forcing for 1,000 years and the last 900 years were used for the analysis presented.
FESOM was set up at a resolution which varies from nominal 1° in the interior of the ocean to 1/3° in the equatorial belt and 24 km (1/4° meridionally) north of 50°N. The ocean surface is discretized with about 127,000 grid points, and 46 vertical levels are used. This mesh has been used in the CORE-II model intercomparison project (e.g., Wang et al., 2016aWang et al., , 2016b and Ocean Model Intercomparison Project Phase 2 (OMIP-2, Tsujino et al., 2020). OpenIFS was used in the TCO159 configuration (reduced Gaussian triangular-cubic-octahedral grid truncated at wave number 159 spectral resolution).
The algorithms for the AMOC computation on unstructured meshes are described in  and Sidorenko, Danilov, Fofonova, et al. (2020). We use vertical and diapycnal velocities to compute the stream functions. Vertical velocity is naturally given by the model output. Diapycnal velocity is diagnosed from the divergence of horizontal flow which is conservatively remapped into prescribed density bins. Here the transports in density space were computed during the run time, which overcomes the need of large storage space, which was necessary for almost all previous work that has used this kind of analysis (see e.g., Megann, 2018).
The ϱ bins are chosen according to Megann (2018) (72 levels for a good representation of deep and bottom waters), complemented with additional density levels to include those presented in Xu et al. (2018). Altogether we use 85 density bins spanning the range of 30.0 < ϱ < 37.2 kg m −3 .

AMOC Frameworks
The AMOC derived in the two different frameworks is presented in the upper panel in Figure 1. The traditional computation in depth coordinates, referred to as z-AMOC, is characterized by a middepth cell centered at ∼1,000 m and a bottom cell centered at ∼4,000 m. The maximum of the middepth cell is located at ∼40°N. It is several degrees south of the main regions of deep convection, which is typical for the depth representation in coarse-resolution runs. The region of closed streamlines around the AMOC maximum, including the diagnosed upwelling slightly south of the maximum, is referred to as the recirculation cell. In Figure 1 (upper left panel), it is confined between 50°N and 30°N and is similar to that in the standalone FESOM configuration described in . The isopycnals at these latitudes (near the subpolar/subtropical gyre interface) essentially deviate from the z-level surfaces, which may lead to spurious dissipation in a transport scheme. Hence,  attribute a part of the recirculation on coarse meshes to spurious numerical mixing and demonstrate that it can be significantly reduced if a higher horizontal model resolution is used.
At a certain depth, z-AMOC is a meridional accumulation of zonally integrated vertical velocities. Hence, vertical velocities indicate the locations of positive (downwelling) and negative (upwelling) contributions to z-AMOC. The bottom left panel in Figure 1 depicts the vertical velocity at 1,000 m, where the maximum of z-AMOC is found. The velocities have been conservatively remapped onto a 4° × 4° grid prior to plotting. On the original mesh the vertical velocity exhibits a noisy structure which masks the main signal. This structure is the consequence of stepwise bottom representation and does not disappear with spatial averaging (see patterns in Katsman et al. [2018]). Similar to , the above-mentioned recirculation cell is caused by the downward flux in the Gulf Stream region and in the Eastern North Atlantic, and the upward flux at Cape Hatteras.
The AMOC in density space ( Figure 1, upper right panel), referred further to as ϱ-AMOC, has a middepth cell located at ϱ = 36.7 kg m −3 , associated with the the formation and southward transport of North Atlantic Deep Water and a less expressed bottom cell at ϱ = 37 kg m −3 , associated with the circulation of the Antarctic Bottom Water. The maximum of the middepth cell, as opposed to z-AMOC, is found at ∼55°N which is close to the regions of deep convection. The associated recirculation is confined between 40°N and 65°N, similar to other studies (e.g., Xu et al., 2018). As is the case with z framework, the negative and positive diapycnal velocities here indicate the locations of, respectively, positive and negative contributions to ϱ-AMOC. The diapycnal velocity across ϱ = 36.7 kg m −3 is presented  Figure 1 (bottom right panel). The buoyancy loss (density increases) associated with recirculation of the middepth cell is found in a zonal band between ∼60°N and ∼65°N. The recirculation is closed by buoyancy gain (density decreases) which takes place within the Gulf Stream and along the route of the North Atlantic Current (NAC). As follows from Figure 1, the isopycnal ϱ = 36.7 kg m −3 is deep in the eastern part of the North Atlantic, but becomes rather shallow and outcrops in the western part at the latitudes of the subpolar gyre. In particular, it is shallow in the Labrador Sea, implying that a substantial part of transformations toward larger density in the waters flowing northward in the main cell of ϱ-AMOC should occur along the eastern part of the northward rim of the subpolar gyre or further north, as recently found in Lozier et al. (2019) for the OSNAP section. This clear picture is blurred in the z-AMOC which combines the flow above and below the isopycnal ϱ = 36.7 kg m −3 at latitudes of the subpolar gyre.
Another advantage of using a density framework is that it allows one to distinguish between contributions from (1) surface buoyancy forcing Ψ s , (2) internally caused contributions through mixing and cabbeling Ψ I , and (3) the rate of volume change (or drift) between the density surfaces   V t / . For computational details of these individual contributions we refer to  and Appendix 1 in Sidorenko, Danilov, Fofonova, et al. (2020). It is worth noting that the total AMOC is the sum of all these terms: ϱ-AMOC = Ψ s + Ψ I +   V t / ≈ Ψ s + Ψ I where we neglect (this will be additionally addressed in the discussion part) the contribution from the drift term due to the long-time average used in this paper.
The surface-forced diapycnal water mass transformation (Ψ s ) and the interior transformation (Ψ I ) are shown in Figure 2 (see Sidorenko, Danilov, Fofonova, et al., 2020; for computation details). Ψ s is characterized by three main cells which are all within the upper limb of the AMOC and their difference mostly reflects the fact that they are in different latitudinal circulation regimes. The three cells are centered at ϱ = 30.95, ϱ = 36.52, and ϱ = 36.89 kg m −3 . The maximum of Ψ I (∼15 Sv) is found at 55°N and indicates that the recirculation cell in ϱ-AMOC is largely maintained through internal transformations. It is, however, found at slightly higher densities ϱ = 36.77 kg m −3 (compared to ϱ = 36.7 kg m −3 for the maximum of ϱ-AMOC), which also highlights the role of Ψ s . The surface transformations in the tropics increase buoyancy and are well compensated by transformations toward smaller buoyancy induced by mixing. The net effect is that the water in the upper AMOC which flows across the equator becomes lighter (as described in Xu et al. [2018]). The lower positive cell in the interior stream function connects the two surface-forced cells associated with the deep water production in the Irminger Sea and south of Iceland (ϱ = 36.7 kg m −3 ) and in the Labrador Sea (ϱ = 36.89 kg m −3 ).
One of the prominent differences between the two frameworks is the location of the AMOC maxima which is found at different latitudes. While both frameworks show predominant sinking and buoyancy loss north of 55°N and upwelling and buoyancy gain at the Gulf Stream separation area near Cape Hatteras, there is no coherent signal found at the interface between the subtropical and subpolar gyres. There the diapycnal velocity along the NAC acts toward buoyancy gain and contributes to the decrease of ϱ-AMOC, defining the southern end of the recirculation cell. From inspecting Ψ s and Ψ I in Figure 2, we may conclude that the NAC front is characterized by large internal transformations taking place below ϱ = 36.7 kg m −3 . Contrarily, vertical velocity depicts sinking along the NAC and marks the northern end of the recirculation cell in z-AMOC. The sinking is associated with the inclination of isopycnals at the boundary of the gyres where the along-isopycnal flow allows for non-zero vertical component in z representation. Hence, the position of the NAC matches the latitude where ϱ = 36.7 kg m −3 starts to flatten at a depth of ∼1,000 m toward the south, remaining largely flat south of 30°N (note the depth contours of ϱ = 36.7 kg m 3 in Figure 1, bottom right panel).
In general, we find that the mean AMOC in both frameworks in our coupled model simulation looks qualitatively similar to those in ocean-only simulations (e.g., Xu et al., 2018). The fact that the AMOC maxima are found at different latitudes in the two frameworks points to differences in the underlying mechanisms. One would expect that using ϱ-AMOC is physically more meaningful as it directly accounts for water mass transformations between different density classes (e.g., Johnson et al., 2019;.

Subtropical and Subpolar Maxima
Traditionally, the AMOC variability is quantified with the time series of the stream function maximum at a certain latitude. In some studies, the time series of the subtropical AMOC maximum is used. There the latitude is often chosen to be either 30°N, where the AMOC matches the center of the overturning cell (see e.g., Reintges et al., 2017), or 26.5°N, the location of the observational MOCHA-RAPID array (e.g., Cunningham et al., 2007). Other works address the subpolar AMOC, which is often computed at 40°-45°N, at the location of the z-AMOC maximum. The subtropical and subpolar variabilities differ in amplitude, are lagged with each other by some years and may poorly correlate at year-to-year timescales (e.g., Bingham et al., 2007). In our simulation, the largest variability is found for the subpolar AMOC at the locations of the respective maxima. As we have shown above, the location of subpolar maxima for depth and density frameworks in our experiment are at ∼40°N (∼24 Sv near subpolar/subtropical gyre interface) and ∼55°N (∼26 Sv south of the main convection areas), respectively. The standard deviations of the AMOC maxima reach 1.33 Sv in z-AMOC and 1.68 Sv in ϱ-AMOC. At 30°N, the mean AMOC drops to ∼15 Sv in both frameworks and a smaller variability is detected with a standard deviation of 0.93 Sv. The fact that the subpolar AMOC is stronger than the subtropical one, suggests that there is significant transformation taking place between the two regions. As is the case with variability, only a part of anomalous diapycnal transformations which are largely located in the northern North Atlantic and cause the fluctuations in AMOC is communicated with the subtropical AMOC.
The variability of the subpolar AMOC maxima is shown in Figure 3 (upper panel). From a wavelet analysis (not shown) we found that the ϱ-AMOC has a slightly stronger decadal and multidecadal variability and, counterintuitively, multidecadal peaks in z-AMOC and ϱ-AMOC occur at different times. The correlation between both frameworks at a year-to-year timescale is only ∼0.34 for 0-lag. It is symmetric about 0-lag and increases to 0.47 for ±1-lag and then drops to 0.1 already at ±4-lag. This is an interesting result, illustrating that not only the subpolar maximum and its position but also the variability is affected by the choice of framework. As we mentioned in Section 3, compared to ϱ-AMOC the maxima of z-AMOC are found further south at ∼40°N, near the subpolar/subtropical gyre interface.
The variability of the subtropical AMOC maxima is shown in Figure 3 (bottom panel). The correlation between both frameworks is 0.98. This reflects the fact that the density surface is nearly flat across the basin Subtropical AMOC is computed at 30°N. Subpolar AMOC is computed as AMOC maxima north of 40°N and is located at ∼40°N in z and at ∼55°N in ϱ frameworks. Both maxima have been computed within the entire depth (density) range. For visualization, a 5-year moving average has been applied to the time series prior to plotting (note that the linearly detrended annual time series were used in the analysis elsewhere). south of 30°N (note the depth contours in Figure 1, bottom right panel) and both frameworks coincide. This also agrees with the findings by Zhang (2010) who studied the latitudinal dependence of AMOC in both frameworks.

AMOC Composites
In order to identify the patterns associated with high and low AMOC states we apply a simple composite analysis. Years where the AMOC deviates from the mean by more than one standard deviation have been considered. All fields as well as the time series of the AMOC maxima have been linearly detrended prior to the analysis. Differences between composite patterns of AMOC anomalies computed for the periods of high (AMOC + ) and low (AMOC − ) values of subpolar AMOC maxima are shown in Figure 4 for the two frameworks. The spatial pattern in the z framework takes the form of a recirculation cell of ∼6 Sv centered at ∼40°N. Its vertical extension to the depth of 4,500 m also indicates the engagement of a bottom cell which becomes weaker during AMOC + . The spatial pattern in ϱ framework is also expressed by a positive recirculation cell centered at 55°N within the density range 35.8 kg m −3 < ϱ < 36.97 kg m −3 and has a similar maximum value of ∼6 Sv. Detrending the fields prior to computing the composites appears to be essential for the density framework analysis. Without it the composites will be distorted by the effects which arise from the fluctuations in volume below respective isopycnals during ϱ-AMOC + and ϱ-AMOC − phases. For obtaining annual-mean density transformations this motion of isopycnals shall be subtracted from the total signal. Since the present study was initially meant for AMOC and not transformation analysis, these data have not been stored. Therefore, what we address below as total and internal transformation might be affected by the motion of isopycnals. Yet, the inconsistency was verified not to affect the presented results. We address this issue later in the discussion section. Figure 5 shows the difference between composite patterns of surface buoyancy forced (∆Ψ s ) transformations for different lags (in years) with respect to the subpolar ϱ-AMOC index. For the negative lags (Ψ s leads) we observe positive anomalies which are confined between 36.7 kg m −3 < ϱ < and 36.97 kg m −3 south of 55°N. This indicates that the surface forcing acts to reduce the buoyancy of this density range. At lag 0, when the AMOC reaches its maximum, and thereafter the pattern is reversed, indicating a decrease in buoyancy of water with densities ϱ > 36.7 kg m −3 . Concurrently, the dipole anomaly centered at ∼36.89 kg m −3 (most expressed at lag 0) points to the shift of the respective deep water production cell (seen as a bottom cell in Figure 2, right) toward lighter densities.
Neglecting the model drift and the change of fluid volume under the isopycnals one could write ∆Ψ ϱ = ∆Ψ s + ∆Ψ I . In Figure 6, we show the annual maps for ∆Ψ I at different lags. Note that the sum of ∆Ψ I and ∆Ψ s at lag 0 is equal to the signal in the right panel in Figure 4. Comparing ∆Ψ I and ∆Ψ s reveals that ∆Ψ s alone explains only a small portion of variability in ϱ-AMOC, especially at lag 0. Furthermore, surface transformations affect the variability of ϱ-AMOC in the density classes which are far larger than the density at which the ϱ-AMOC maximum is located. Thus, the interior-mixing-induced transformations are not only   responsible for the existence of the mean recirculation in ϱ-AMOC, forming its subpolar maximum, but also for its variability. This highlights the importance of the interplay between surface buoyancy flux and interior-mixing in the higher density classes. Note that in our analysis we do not separate the effects of the nonlinear equation of state (cabbeling and thermobaricity) when computing Ψ I . Hence, the presented Ψ I combines all diabatic processes which take place.

Composites for Vertical Transport and Diapycnal Transformations
As we have mentioned earlier, AMOC is a meridional accumulation of zonally integrated vertical (z-AMOC) or diapycnal (ϱ-AMOC) velocities. Hence, AMOC variability is directly related to the variability in vertical or diapycnal velocities. A simple composite analysis for the velocities at depths where AMOCs reach their maxima reveals the locations associated with AMOCs fluctuations. Composite patterns of vertical velocity at 1,000 m during subpolar z-AMOC + , z-AMOC − , and the difference between both are presented in Figure 7. Combining these patterns with the z-AMOC variability pattern (left panel of Figure 4) we conclude that subpolar z-AMOC fluctuations are expressed by stronger sinking in the Gulf Stream and NAC regions, and stronger upwelling around Cape Hatteras during z-AMOC + as compared to z-AMOC − . Interestingly, the difference between composites slightly changes sign along the NAC. It indicates that the z-AMOC variability might be partially linked to the shift in the position of the NAC. This in some way matches the findings by Frankignoul et al. (2013) who noticed a southward shift of the Gulf Stream during AMOC intensification.
As we have shown in Section 4.2 the variability of AMOC takes the form of a recirculation cell which in the ϱ framework is largely maintained through internal transformations. We choose two density levels ϱ = 36.7 and ϱ = 36.89 kg m −3 relevant for this variability. In Figure 8, we present the ϱ-AMOC + and ϱ-AMOC − composite maps of diapycnal velocities and their differences across these isopycnals. The variability of the upper cell in ϱ-AMOC (positive pattern in Figure 4, right panel) is associated with the westward shift in the position of buoyancy loss in the Irminger Sea and more buoyancy loss at the southern tip of Greenland during ϱ-AMOC + as compared to ϱ-AMOC − . This cell largely recirculates and is closed by more buoyancy gain along the NAC. The variability of the lower ϱ-AMOC at ϱ = 36.89 kg m −3 is characterized by two anomalies: a positive one (centered at ∼55°N) is associated with the dipole anomaly of reduction of the buoyancy gain at the southern tip of Greenland and the reduction in buoyancy loss in the southern Labrador Sea (LS) and a negative one (continues south of ∼45°N) associated with the anomalous lightward diapycnal velocities along the path of the NAC during ϱ-AMOC + . The reduction in buoyancy loss in the southern Labrador Sea and along the NAC during ϱ-AMOC + is dominant. Hence the lower cell of ϱ-AMOC variability spreads further south. In Section 5, we speculate that this large latitudinal spread can be partially linked to the isopycnal motion which is spuriously imprinted as diapycnal transformation in the analysis.
Interior diapycnal velocities, inducing AMOC cells, redistribute the surface transformations which happen in succession through all density classes (at all levels). Hence, in Figure 9 we present the mean surface transformations at two centers of the ϱ-AMOC pattern. At ϱ = 36.7 kg m −3 , it is expressed by buoyancy gain and loss occurring at different locations. Buoyancy gain is found along the East Greenland Current (EGC) and Labrador Coastal Currents (LCC). It is primarily driven by the freshwater contribution to buoyancy flux (not shown). Concurrently, buoyancy loss occurs east of the buoyancy gain in the Irminger Sea, south  of Iceland, in the Norwegian Sea and to a lesser degree in the Greenland Sea. Heat flux (not shown) is the main contributor to buoyancy change in these regions. At ϱ = 36.89 kg m −3 the surface transformation is mainly expressed by buoyancy loss in the entire LS, in the Norwegian Sea and to a lesser degree in the Greenland Sea. Buoyancy gain is found only in the narrow area of the Denmark Strait. As one would expect, the map of mean surface transformations indicates that the LS is largely responsible for the deep water formation in the North Atlantic.
The differences between composites of surface transformations are presented in Figure 10 for different lags. At ϱ = 36.7 kg m −3 they show only minor anomalies in the LS before lag 0 which indicates that mainly the higher densities are exposed to the surface there. At lags ≥0, at and after ϱ-AMOC reaches its maximum, lower densities are found along the periphery of the LS and the EGC and negative buoyancy anomalies start to appear at these locations. Concurrently, positive buoyancy anomalies are found in the Irminger Current and south of Iceland. Change in the heat flux is the prime driver for these anomalies. Interestingly, a similar pattern is also found in the composites for diapycnal velocities.
A different behavior is seen at ϱ = 36.89 kg m −3 . There, for negative lags we observe a northward shift of the buoyancy loss in the LS which is expressed by large negative buoyancy anomalies in the northern LS and less pronounced positive buoyancy anomalies in the southern LS. At lags 0 and later, the patterns change sign and mainly positive buoyancy anomalies are found along the periphery of the LS. Negative anomalies, however, still exist at lag 0 in the Northwest LS and explain the appearance of the dipole pattern in ∆Ψ s centered around ϱ = 36.89 kg m −3 (see Figure 5 for 0 lag).
Note that the anomaly of buoyancy gain at ϱ = 36.89 kg m −3 and lower densities as imposed by the surface transformations at lags ≥0 may contribute to buoyancy loss at the lighter density classes in the presence of "upward" internal diapycnal transformation. Upward diapycnal velocities at ϱ = 36.89 kg m −3 are persistently present at the southern tip of Greenland (see bottom left and middle panels in Figure 8). The importance of buoyancy loss contribution into the middensities from the lower densities was already addressed in the analysis of the LS mean transformations (e.g., Xu et al., 2018). Here we speculate that this process is also associated with the variability of the AMOC.

AMOC Relation to Sea Level Pressure and Mixed Layer Depth
Variability of the AMOC has been shown to be related to atmospheric forcing and to changes in deep water formation (e.g., Frankignoul et al., 2013;Menary, Jackson, & Lozier, 2020). Here we analyze these links in our simulation. The maps of mean atmospheric Sea Level Pressure (SLP) and Mixed Layer Depth (MLD) maximum for the northern winter (JFM) are shown in Figure 11 and resemble well-known patterns. In SLP the Azores high and Icelandic low-pressure centers are located at ∼30°N and ∼60°N, with an averaged pressure difference of ∼30 hPa between both centers. The MLD pattern points to strong diapycnal mixing broadly along the NAC and Irminger Current regions, and in particular in the LS. Some mixing is also found in the Norwegian and Greenland Seas. In Figure 12, we present the differences between AMOC + and AMOC − composites of SLP at different lags for the two frameworks. Most of the SLP difference patterns resemble the North Atlantic Oscillation (NAO) variability pattern. As one would expect, a positive NAO phase (negative SLP anomalies around Iceland, positive SLP anomalies around the Azores) precedes the appearance of the AMOC maxima. In the z framework the positive NAO pattern is most expressed at lag = −3 but nearly vanishes at lag = −1. In contrast, in ϱ framework the pattern occurs at all negative lags, in particular at lag = −2. Differently lagged SLP anomaly patterns possibly reflect the fact that the AMOC maxima are found at different locations in the two frameworks and that it takes more time for the z-AMOC signal to travel from the area of deep convection to the location of the z-AMOC maximum at ∼40°N. For ϱ-AMOC the behavior is more coherent since the position of the maximum is found at ∼55°N where the convection sites are. Interestingly, the two frameworks show distinct behavior at lag 0, depicting large SLP increase centered at ∼40°N in z framework, similar to a (northward-shifted) positive NAO, and at ∼60°N in ϱ framework, similar to a (northward-shifted) negative NAO.
Differences between composites for MLD are shown in Figure 13. It is worth noting that both frameworks exhibit large similarity which persists for positive and negative lags of several years. This highlights that the MLD variability is associated with a longer than year-to-year timescale. Noticeably, negative lags are associated with a northward shift of the MLD in the LS and positive lags with a southward shift. Hence, extensive deep water production in the northern part of the LS and the reduction of this in the southern LS precedes the appearance of AMOC maxima. This behavior agrees with the differences in composites for the surface buoyancy flux across ϱ = 36.89 kg m −3 , representative for the deep water formation (right panel in Figure 10).
In summary we find that at negative lags, the NAO-like patterns can be directly explained by their influence on buoyancy fluxes, with expected consequences for AMOC in both frameworks. Due to the long memory of the ocean state (we refer to MLD here), several preceding years of positive NAO anomaly continually contribute to the AMOC increase. In contrast, at lag 0 the dominant impact is caused by the immediate response to the anomalous Ekman transport: for z-AMOC the lag-0 SLP pattern is associated with easterly winds (Ekman transport to the north) anywhere south of ∼45°N (where the high-pressure center is located), and for ϱ-AMOC the lag-0 SLP pattern is associated with easterly winds (Ekman transport to the north) anywhere between ∼65°N (where the high-pressure center is located) and ∼40°N (where the low-pressure center is located). Our findings match those by Ortega et al. (2012) who attributed high frequencies to local changes in Ekman transport caused by atmospheric modes of variability, and the low frequencies to the deep water production south of Greenland.

Discussion
The present study aims to illustrate the usability of augmenting the depth and density frameworks with the spatial patterns of sinking and diapycnal transformations for studying the AMOC and its variability. We analyzed these patterns and their variability at a depth level (for z-AMOC) and at an isopycnal surface (for ϱ-AMOC) where AMOC reaches its maximum. The analysis requires horizontal binning, especially for vertical velocities which are too noisy in the main regions of upwelling and downwelling because of stepwise bottom topography. Conservative remapping to 4° × 4° boxes appeared optimal; finer bins still showed a rather patchy structure in regions with steep bathymetry. After the remapping step the main areas of sinking and upwelling, which constitute the AMOC, could be clearly identified. In contrast, the diapycnal velocities are less noisy and could be remapped to the finer meshes, but we still chose 4° × 4° boxes to ensure consistency. Depending on the purpose, remapping could be designed in a more elaborate way, for example by splitting the domain into geographical areas for illustrating their relative roles. The only trivial requirement here is that it should be conservative. It is worth noting that even for the coarse-grained data, the standard deviations in diapycnal and vertical velocities may exceed the local signals associated with the AMOC. This reflects the fact that the exact locations of vertical velocities or where the water masses cross the isopycnals vary considerably from year to year within the relevant region. However, given the large number of years contributing to the composites, the presented patterns are still robust (verification analysis is not presented) and physically plausible.
Since the present study was initially meant for AMOC and not transformation analysis, we did not store the model volume drift under isopycnals (the rate of volume change   V t / above the ϱ surface) as discussed in Section 4.2. Thus, the internal transformations which we present in Figure 6 are biased by these adiabatic fluctuations. For the long-term averages used in this paper, this drift is negligible (less than 1 Sv in amplitude) and therefore does not change the mean transformation fields. However, it may bias patterns of their variability. Although the analyzed patterns as well as the AMOC time series have been detrended prior to the analysis we cannot exclude that a part of discussed variability can be caused by adiabatic fluctuations. It is also worth noting that internal transformations redistribute the surface buoyancy fluxes which have been aggregated over some period of time and the partition of ∆Ψ ϱ onto ∆Ψ s and ∆Ψ I becomes (due to possible lags) less apparent if made instantaneously or on a year-to-year basis. Although this partitioning works for mean fields, some deeper look into consistent ways for partitioning the anomalies is required.
It is known that more insight can be gained from the AMOC analysis in density framework as it allows one to distinguish between surface buoyancy forced and internal transformations (e.g., Xu et al., 2018;Zhang, 2010). Our deeper analysis in the density framework reveals that looking at spatial patterns of transformations only at the locations where AMOC reaches its maximum may be sufficient to learn about the mean AMOC, but is insufficient to address its variability. Indeed, the places where transformations through selected isopycnals are large do not imply that the ϱ-AMOC is being modified right there. Surface transformations happen in succession through all density classes (at all levels) and are further redistributed by interior diapycnal transformations. Hence, for addressing the AMOC variability we augmented the composite maps of ϱ-AMOC and its constituents at different isopycnal levels. Only the aggregated analysis allows one to associate deep water production with AMOC changes.
We note that subpolar AMOC maxima in both frameworks are located at different latitudes and hence not only differ in mechanisms of underlying sinking and diapycnal transformations but also disagree on temporal delays. Although the latitudinal propagation of AMOC perturbation signals might be of high interest, it is known to be prone to errors due to model uncertainties and resolution (e.g., Frankignoul et al., 2013;Katsman et al., 2018;Menary & Hermanson, 2018;Menary, Robson, et al., 2020). We plan to analyze latitudinal propagation of AMOC-related anomalies in future work by employing an eddy-resolving setup of the North Atlantic.

Conclusions
We analyzed the AMOC variability in a 900 years climate run by using density and z frameworks. The AMOC variability is nearly identical in both frameworks south of ∼30°N where the isopycnals are nearly flat. In the northern North Atlantic, where isopycnals are sloping, the two frameworks show substantial differences. First, the recirculation cell and hence the maximum of the subpolar AMOC are located at different latitudes for the two frameworks. Second, the variability of the subpolar AMOC maxima is correlated with a coefficient of only ∼0.4 on the annual timescale (0.64 after applying a 5-year moving average), implying that the associated patterns of variability (difference between composites) are different. The ϱ-AMOC emphasizes the role of internal transformations.
The individual analysis of internal and surface-forced constituents of the AMOC reveals that variability is largely driven by internal transformations. These in turn are triggered by the anomalous surface buoyancy flux, that is, the internal transformations modify the water masses which are ventilated through the surface buoyancy flux. Negative heat flux anomalies in the northern LS, reduced freshwater export from the Nordic Seas, and an associated northward shift of MLD there precedes AMOC maxima. Interestingly, AMOC maxima are followed by reversed MLD anomalies. This suggests that, after the gradual build-up of an anomalously thick mixed layer over several years due to positive NOA anomalies, the finally Ekman-transport driven maximum overturning acts to degenerate the mixed layer, leaving behind negative MLD anomalies that persist over several subsequent years.
Another interesting detail is that the anomalous upward diapycnal velocities along the path of the NAC, corresponding to the southern end of the middepth recirculation cell in ϱ-AMOC and related to its variability, occur in the absence of surface transformations in this area. We suggest that it is partially caused by the dissipative truncation errors in horizontal advection which lead to diapycnal mixing in places where isopycnals are sloping. Hence, as in , we attribute this behavior to spurious numerical mixing due to sloping isopycnals that essentially deviate from level surfaces.
We found that the surface buoyancy transformations which precede AMOC maxima are associated with NAO-like SLP anomaly patterns. Furthermore, a more coherent response to atmospheric change is found in the density framework rather than in z. Yet, at the year with high AMOC the dominant impact in both frameworks is caused by the instantaneous Ekman transport. Given that no association between SLP and AMOC is found at positive lags (AMOC leads), we conclude that the atmosphere is the main driver for