Disentangling North Atlantic ocean-atmosphere coupling using circulation analogues

The coupled nature of the ocean-atmosphere system frequently makes understanding the direction of causality difficult in ocean-atmosphere interactions. This study presents a method to decompose turbulent heat fluxes into a component which is directly forced by atmospheric circulation, and a residual which is assumed to be primarily `ocean-forced'. This method is applied to the North Atlantic in a 500-year pre-industrial control run using the Met Office's HadGEM3-GC3.1-MM model. The method identifies residual heat flux modes largely associated with variations in ocean circulation and shows that these force equivalent barotropic circulation anomalies in the atmosphere. The first of these modes is characterised by the ocean warming the atmosphere along the Gulf Stream and North Atlantic Current and the second by a dipole of cooling in the western subtropical North Atantic and warming in the sub-polar North Atlantic. Analysis of atmosphere-only simulations confirms that these heat flux patterns are indeed forcing the atmospheric circulation changes seen in the pre-industrial control run. It is found that the Gulf Stream plays a critical role in the atmospheric circulation response to decadal ocean variability in this model. More generally, the heat flux dynamical decomposition method provides a useful way to establish causality in ocean-atmosphere interactions which could easily be applied to other ocean basins and to either models or reanalysis datasets.


Introduction
Understanding the two-way interaction between extratropical sea surface temperature (SST) variability and atmospheric circulation is a complex problem.It is also an eminently practical one as ocean variability is a major source of prediction skill for forecasts of the atmosphere on subseasonal to decadal time scales (Meehl et al. 2021;Merryfield et al. 2020).
Turbulent heat exchange is the primary process by which the ocean and atmosphere transfer heat between each other (Cayan 1992).On interannual and shorter time scales, the atmosphere largely governs ocean-atmosphere covariability via modulation of air temperature, specific humidity, and nearsurface wind speed, and these, in turn, modify surface turbulent heat fluxes (Q, i.e., latent plus sensible heat fluxes).At decadal time scales and longer, the ocean dominates as it integrates atmospheric variability and responds via changes to ocean circulation and ocean heat transport, hence altering SSTs and Q (Gulev et al. 2013).
The North Atlantic exhibits large decadal SST variability, and there has been considerable debate over the extent to which this is driven by atmospheric circulation (Clement et al. 2015;Cane et al. 2017), internal ocean variability (Zhang et al. 2019;O'Reilly et al. 2016), and external forcing (Booth et al. 2012).Attributing the drivers of low-frequency variability is also complicated by the fact that current climate models likely have an inadequate representation of relevant mechanisms as they show too little decadal variability in both atmospheric circulation (Bracegirdle et al. 2018;Simpson et al. 2018;O'Reilly et al. 2021) and ocean circulation (Kim et al. 2018).On the other hand, some low-resolution models tend to overestimate decadal SST variability in the subpolar North Atlantic (Patrizio et al. 2023).
Relatedly, the atmospheric circulation response to midlatitude SST variability in climate models is likely too weak, leading to overly weak signals in seasonal to decadal predictions (Eade et al. 2014;Scaife and Smith 2018;Smith et al. 2020).The cause of this problem, termed the "signal-to-noise paradox," is not understood, but one leading hypothesis concerns inadequate model representation of ocean-atmosphere coupling (Scaife and Smith 2018).Consequently, exploring the fidelity of models when it comes to ocean-atmosphere coupling may be necessary for making progress on the signalto-noise paradox and improving climate predictions.However, doing so requires evaluation tools which can be applied to both models and observation-based datasets.
To separate oceanic effects on the atmosphere from atmospheric effects on the ocean, studies use a variety of methods.These include 1) low-pass filtering data to isolate time scales at which the ocean dominates, 2) using lagged correlation analysis, and 3) performing atmosphere-only experiments.However, as discussed above, current models are deficient at capturing the response to midlatitude SSTs, which presents challenges with the latter option.Moreover, the shortness of the observed record means that low-pass filtering leaves only a few degrees of freedom and may smooth over the effects of internal ocean variability, which occurs on short time scales.
This study presents a method to separate the ocean-forced component of Q from the atmospheric circulation-forced component.The method has been designed such that it does not require any low-pass filtering and can be applied to both models and observation-based datasets.In this study, the method is only applied to model data as this provides a more controlled setup in which there is a long data record and there is no observational uncertainty associated with variables such as Q.Testing with observation-based datasets, such as reanalysis, is reserved for a future paper.The method involves the use of circulation analogs to identify the component of Q directly associated with atmospheric circulation variability, diagnosing the ocean-forced component as the residual.We apply this method to a 500-yr preindustrial control run.We also test sensitivities of the results by applying it to simulations of the same model with observed external forcings from 1850 to 2014.
The datasets and methods are described in section 2. The method is then applied to a preindustrial control (piControl) simulation in section 3, and the variability of the atmospheric circulation-related and residual Q fields is compared.The leading modes of the decomposition are then examined in section 4, followed by an analysis of any atmospheric circulation response to the Q modes in section 5. Discussion and conclusions are provided in section 6.This is followed by appendixes in which the sensitivity to the choice of parameter values (appendix A) and sensitivity to the presence of variable external forcing and length of the dataset (appendix B) are examined.

a. Data
We analyze simulations made using the Met Office (UKMO) HadGEM3-GC3.1 model (Williams et al. 2018) for which the North Atlantic ocean-atmosphere coupling has been extensively analyzed (e.g., Lai et al. 2022;Khatri et al. 2022).The model consists of coupled ocean, atmosphere, land, and sea ice models.In this study, we utilize the medium (MM) resolution version, which is run with an N216 (grid spacing of approximately 60 km) grid in the atmosphere with a horizontal ocean resolution of 0.258 (ORCA025).The model has 75 vertical levels in the ocean and 85 in the atmosphere.We analyze a 500-yr piControl simulation using HadGEM3-GC3.1-MM.This is a fully coupled simulation which has climatological external forcing.Further details can be found in Menary et al. (2018).The piControl run simulates Atlantic multidecadal variability (AMV) with a 60-80-yr period, consistent with observations, and a slightly weaker Atlantic meridional overturning circulation (AMOC) at 26.58N and at subpolar latitudes with respect to RAPID (Menary et al. 2018) and overturning in the subpolar North Atlantic program (OSNAP) observations (Menary et al. 2020), respectively.
We also briefly analyze an atmosphere-only experiment, known as High Resolution SST-present (highresSST-present), taken from the High Resolution Model Intercomparison Project (HighResMIP) (Haarsma et al. 2016).In this experiment, HadGEM3-GC3.1-MMhas been forced with observed SSTs (taken from HadISST2; Titchner and Rayner 2014) from 1950 to 2014 and with historical greenhouse gas and aerosol forcings.A total of three different ensemble members were run, all with slightly different initial states.Four members of historical coupled simulations from the same model (Andrews et al. 2020) are also analyzed in appendix B. These use historical forcings including greenhouse gas and aerosol variations spanning 1850-2014.

b. Circulation analogs
To attribute Q anomalies to atmospheric or oceanic forcing, we apply a circulation analog method similar to that used by Deser et al. (2016) andO'Reilly et al. (2017).The concept of comparing similar circulation states was first developed in the context of statistical weather prediction by Lorenz (1969) and later van den Dool (1994) and van den Dool et al. (2003).More recently, it has been used to study the degree to which atmospheric circulation trends have played a role in observed temperature trends (Cattiaux et al. 2010;Wallace et al. 2012;Deser et al. 2016).
The circulation analog method attempts to estimate the component of a temporally and spatially varying variable, that is directly and simultaneously associated with changes in atmospheric circulation.In our case, we decompose Q into two components, where Q CIRC is the atmospheric circulation-related component of Q, and Q RESIDUAL is the residual.We interpret Q RESIDUAL to be primarily due to ocean variability and the persistence of SST anomalies forced by atmospheric circulation in previous months.The Q RESIDUAL will also include the effects of small-scale circulation features that do not project clearly on the monthly sea level pressure (SLP) fields.However, one can assume that these are random errors that should cancel out over large enough samples or over long enough periods.Additionally, Q RESIDUAL will vary due to the radiative warming of the atmosphere, for example, through externally forced radiative changes.There is no variable external forcing in the piControl run, but we remove the effects of external forcing when applying the method to historical simulations by regressing out global-mean SST anomalies from the SLP at each grid point.We note that removing the global mean can have the effect of removing internal variability related to the Pacific decadal oscillation (Deser and Phillips 2023); however, this is less of a problem here as our focus is the North Atlantic.The historical runs are discussed further in appendix B. We also perform a linear detrending of piControl anomalies prior to reconstructing SLP in order to remove any drifts in the model.Our method begins by taking the seasonal-mean SLP anomalies for a particular year, say 1901, and calculating the area-weighted Euclidian distance between this and all other years over the North Atlantic region (defined as 208-758N, 908W-08).A subsample of size N s from the N a most similar years is then taken, and a weighted sum of this sample of anomaly maps is used to reconstruct the original (1901) anomaly field.That is, we calculate the values of weights a k , which minimize the following: where SLP ij is the SLP map for year i and grid point j within the North Atlantic region.In this case, i represents the year which we are trying to fit to and k is all other years in the subsample.The term u j is the latitude of grid point j; hence, cos(u j ) is the area-weighting factor.Performing this minimization gives a set of N s weights a k for each of the subsampled years.The term Q CIRC is then calculated by summing the Q anomalies for the same years multiplied by the corresponding weights calculated for the SLP anomalies, i.e., The resampling procedure is repeated N r times to obtain N r reconstructions of SLP and Q CIRC , which are then averaged to find a best estimate of these two quantities.Here, N s , N a, and N r are taken to be 50, 80, and 100, respectively.This entire process is then repeated for all years in the dataset.

c. Linear decomposition of Q anomalies
The term Q is composed of sensible Q S and latent heating Q L terms which can be represented using the bulk formulas Q L 5 rC e LUDH and Q S 5 rC p C H UDT, respectively.Here, r is the air density, U is the near-surface wind speed, C p is the heat capacity of water, L is the latent heat of evaporation, and C e and C H are the transfer coefficients.The DT 5 T s 2 T a is the air-sea temperature difference, and DH 5 H * s 2 H a is the difference between the saturation specific humidity at the sea surface (H * s ) and the near-surface specific humidity (H a ).Here, subscript s represents the sea surface and a is the atmosphere.
A linear decomposition of Q (e.g., Alexander and Scott 1997;Du and Xie 2008;He et al. 2022) yields where overbars represent the time-mean quantities and primes are the anomalies with respect to the time mean.This decomposition assumes that U DH , , UDH andU DT , , UDT , which are both good approximations at the monthly time scale (Alexander and Scott 1997).

d. Indices
The NAO index is calculated as the first empirical orthogonal function (EOF) of December-March (DJFM)-mean SLP over the regions 208-808N and 608W-08, calculated using the Python package "eofs" (Dawson 2016).The AMOC index is defined, following Lai et al. (2022), as the annual-mean Atlantic overturning streamfunction (in depth space) at 458N and 1000-m depth.

e. Significance testing
Many of the time series considered in this study exhibit serial correlation, which reduces the effective number of degrees of freedom.To account for this, following Wilks (2011), we replace the number of time steps N with the effective number of time steps, N Eff 5 N[(1 2 r)/(1 1 r)], when calculating statistical tests.Here, r is the lag-1 autocorrelation.

Evaluating the dynamical decomposition of the
Q method a. Case study year We now apply the circulation analog method described in section 2 to the North Atlantic, over a box bounded by latitudes 208-758N and longitudes 908W-08.An example of the decomposition is shown in Fig. 1 for the winter (DJFM) of model year 1911 in the HadGEM3-GC31-MM piControl simulation.The circulation-related SLP field, marked by a positive NAO-like pattern, is, by construction, almost identical to the full field over the North Atlantic region (Figs.1a,c); however, this is not the case outside of the North Atlantic (not shown).The Q anomalies (defined as the positive upward throughout this study) indicate anomalously high heat loss from the ocean to the atmosphere over a horseshoe-shaped region involving the subpolar North Atlantic and eastern subtropical North Atlantic (Fig. 1b).The dynamical decomposition suggests that a substantial proportion of this is related to the atmospheric circulation, including heat loss over the western subpolar and subtropical North Atlantic (Fig. 1d).The Q RESIDUAL anomalies are of similar magnitude to Q CIRC anomalies and are characterized by ocean heat loss over the eastern subpolar region and heat gain over the western subtropics, with a northward shift of the Gulf Stream (Fig. 1f).
We now examine the contributions of different factors (wind speed, air-sea temperature, and humidity differences) to Q anomalies via linear decomposition (see section 2).Note that the sums of the linearized terms explain the vast majority of the total Q anomalies (Fig. S1f in the online supplemental material).Strong near-surface wind speeds and air-sea temperature and specific humidity differences all contribute to the subpolar heat loss in this year (Figs.2a-c), which is largely driven by atmospheric circulation (Fig. 1e).The Q anomalies related to air-sea temperature differences are largest over the Labrador Sea, where the circulation advects cold air from the North American continent.In contrast, Q anomalies linked to wind speed and latent heating are larger in the eastern subpolar North Atlantic than in the west.The largest contribution to the positive Q anomalies in the Gulf Stream region comes from latent heating, while wind speed variations dominate Q anomalies in the subtropics.Overall, the pattern of Q anomalies related to latent heating is similar in structure to Q RESIDUAL , suggesting that latent heating is the main contributor to Q variability, which is not directly related to atmospheric circulation.

b. Interannual to decadal variability
To test the dynamical decomposition method more systematically, we calculate the correlation, at each grid point, between DJFM-mean SST anomalies and the components of the Q dynamical decomposition.Consider that if a warm SST anomaly is primarily the result of warming by the atmosphere, then the anomalous Q is negative, while a cool SST anomaly will be associated with a positive upward heat flux anomaly.Conversely, if the SST anomaly is warming or cooling the atmosphere, having formed through alterations to ocean circulation or by atmospheric forcing at least a month or two previous, then the sign of the SST anomaly should be the same as that of the anomalous heat flux.That is, a negative As anticipated, the Q CIRC component is negatively correlated with SST across the North Atlantic, suggesting a primarily downward influence (Fig. 3a).In contrast, SSTs are positively correlated with the Q RESIDUAL (Fig. 3b), suggesting a largely upward influence of Q anomalies.For reference, the full Q field shows that over the extratropical North Atlantic, SST variability tends to warm the atmosphere more than the atmosphere warming the SSTs, whereas the influence is generally downward on the edge of the tropical Atlantic (Fig. 3c).The Q variability in the Gulf Stream region stands out as being particularly dominated by the ocean (Fig. 3c).
Further evidence that Q RESIDUAL is forced by the ocean rather than the atmosphere is given by lagged gridpoint correlations between Q RESIDUAL and SSTs in the extended seasons prior to and following DJFM, respectively, August-November (ASON) and April-July (AMJJ), see Fig. S2.In particular, SST anomalies have developed prior to DJFM, suggesting that Q RESIDUAL is driven by persistent ocean variability (Fig. S2a).In contrast, there is no significant correlation between DJFM Q CIRC and ASON SSTs (Fig. S2d).
In addition to the direct Q response to atmospheric circulation variability, ocean surface currents and, therefore, Ekman heat transport will change on submonthly time scales in response to altered surface winds (e.g., Alexander and Scott 2008).This modifies the SSTs and therefore affects Q, particularly in regions of strong SST gradients (Deser et al. 2010).To examine how this process is partitioned in the Q dynamical decomposition, following Deser et al. (2010), we calculate the effective Ekman heat flux.This is given as r 0 c p HU Ek ?=T, where T is the SST, r 0 is the water density, c p is the specific heat capacity of water, and H is the depth of the Ekman layer, taken as 50 m.The Ekman current anomaly, U Ek 5 (U Ek , V Ek ), is given by s ], where t (y) s and t (x) s are the meridional and zonal components of the surface wind stress, respectively, and f is the Coriolis parameter.
Simultaneous gridpoint correlations between the Ekman heat flux and Q CIRC and Q RESIDUAL indicate that the Ekman effect is incorporated in Q CIRC and not in Q RESIDUAL (Figs. 4a,b).We argue that this makes sense on the basis that the anomalous surface Ekman transport is a direct and nearly instantaneous response to atmospheric circulation variability.It is worth noting that this effective Ekman heat flux is a much smaller effect ('20%-30%) compared to total Q variability (Fig. S3).
We now examine spatial variations in the magnitude of Q variability driven by the atmosphere and ocean on different time scales.The majority of the simulated Q variability on interannual to decadal time scales in the subpolar North Atlantic and particularly the Labrador Sea is associated with Q CIRC (Figs. 5a,c,d,f).This is likely due to atmospheric circulation modulating the advection of cold air from the North American continent over the ocean, raising the air-sea temperature contrast (also see Fig. 2b).However, this may not be the only factor determining the large atmospheric influence on Q in this region; for instance, this region is also the location of the maximum mixed layer depth (Buckley et al. 2019).For the eastern subpolar North Atlantic, the decadal Q RESIDUAL variance is larger than Q CIRC (Figs. 5d,e), likely due to the effects of internal ocean variability.In contrast, atmospheric circulation contributes more to Q variance in this region on interannual time scales (Figs.5a-c).The Q RESIDUAL shows a similar magnitude of variability to Q CIRC along the North Atlantic Current (NAC) at interannual time scales, and Q RESIDUAL shows larger variability in the Gulf Stream region on both time scales (Figs.5b,e).However, while the Gulf Stream Q variability is not dominated by simultaneous atmospheric circulation variability, it is possible that atmospheric circulation has a lagged effect on the Gulf Stream through induced changes in ocean circulation (McCarthy et al. 2018).The next section examines the primary modes of Q variability associated with the components of the Q decomposition and relates these to simultaneous and lagged patterns of atmospheric circulation.

Modes of Q variability
To understand the spatial patterns of variability associated with the Q decomposition, we perform an area-weighted EOF analysis separately for Q, Q CIRC , and Q RESIDUAL over the same region used to calculate the decomposition (208-758N, 908W-08).The first EOFs of Q and Q CIRC are both characterized by a tripole pattern (Figs.6a,e), associated with the positive phase of the NAO (Figs. 6c,g).The second EOFs of Q and Q CIRC are again very similar to one another and consist of enhanced Q over the central North Atlantic (Figs. 6b,f), linked to the east Atlantic pattern (Figs.6d,h).The first two EOFs of Q CIRC explain more variance (EOF1: 34.5%, EOF2: 17.5%) than those of Q (EOF1: 25.0%, EOF2: 12.6%), likely because Q also includes variability which is unrelated to atmospheric circulation.
The Q RESIDUAL EOF1 is marked by an anomalous positive Q along the NAC and Gulf Stream with a negative anomaly south of the Gulf Stream (Fig. 6i).The positive Gulf Stream anomaly is located slightly to the north of the negative Q CIRC EOF1 anomaly (Fig. 6e).The second EOF has positive Q anomalies to the southeast of Greenland, with a weaker negative anomaly in the subtropics (Fig. 6j).The structure of Q RESIDUAL EOF2 bears some resemblance to Q CIRC EOF1 except that the subpolar anomaly is located in the eastern rather than in the western subpolar North Atlantic, and the positive component west of North Africa is largely absent (Fig. 6j).The Q RESIDUAL EOFs show no simultaneous correlation with SLP by construction (Figs.6k,l).
The Q RESIDUAL EOF1 is somewhat reminiscent of the "slow" response at 3-4-yr lags to NAO forcing found by Khatri et al. (2022) using the same model but with decadal hindcast data (cf.their Fig. 2).They found that the initial "fast" response to the NAO caused by wind stress and Q anomalies is followed by a slower adjustment to SSTs involving a strengthened overturning circulation.To investigate the relationship of Q RESIDUAL EOF1 with ocean overturning, we regress the Atlantic overturning streamfunction onto the principal component (PC) time series associated with Q RESIDUAL EOF1 in Fig. 7 at various lags.This analysis confirms that Q RESIDUAL EOF1 is associated with a strengthened AMOC (Fig. 7c).The AMOC appears to strengthen about five years prior to the peak of Q RESIDUAL EOF1 (Fig. 7b) and weakens following the Q RESIDUAL EOF1 peak (Fig. 7d).The Q RESIDUAL EOF2 is also linked to AMOC variability, though the peak AMOC strengthening occurs a few years prior to the Q RESIDUAL EOF2 peak (Figs.7e,f).In contrast to Q RESIDUAL , the lagged regression of the Atlantic overturning streamfunction onto the Q CIRC EOF modes indicates that these modes are not preceded by North Atlantic Ocean overturning anomalies (Figs.S4a,b,e,f).Instead at lag 0, a barotropic Ekman response is observed (Figs.S4c,d).This difference in regression patterns between Q RESIDUAL and Q CIRC provides further evidence that Q RESIDUAL is associated with slow ocean variability.
The AMOC variability itself is partly driven by the NAO with slower ocean overturning integrating NAO variations over multiple years (e.g., McCarthy et al. 2015;O'Reilly et al. 2019).This can be seen in Fig. 8a, which shows that the (annual-mean) AMOC is strengthened in the 3 years following a positive winter NAO.The Q RESIDUAL EOF1 is most positively correlated with NAO variability when the NAO leads by 1 year, suggesting that the NAO strengthens the AMOC (Fig. 8a), which drives the transport of warmer water poleward and subsequently alters Q (Fig. 8b).The signature of lower-frequency ocean variability in Q RESIDUAL is also seen in the normalized power spectra of the Q RESIDUAL PCs.The Q RESIDUAL EOF1 shows notable peaks in power for periods of 20 and 40 years, neither of which is seen for Q and Q CIRC (Fig. 8c).The Q RESIDUAL EOF2 also shows a peak at around 40 years (Fig. 8d).These peaks in Q RESIDUAL spectra are significantly different from a red noise model fitted to the data at the 95% level.Note that Q and Q CIRC show considerably more variability than Q RESIDUAL overall, but we examine normalized power spectra to indicate which frequencies relatively dominate each of the EOF PCs.

Circulation response to heat flux anomalies
In this next section, we investigate whether or not there is an atmospheric circulation response to the Q RESIDUAL EOF modes.By construction, the dynamical decomposition removes any internal or remotely forced atmospheric variability as well as any variability which is a response to local ocean variability.Hence, it is not possible to infer using simultaneous regression analysis whether ocean variability drives an atmospheric circulation response.However, if we assume that the atmospheric circulation response to midlatitude ocean variability is relatively weak compared to internal atmospheric variability (e.g., Kushnir et al. 2002), the structure of any true ocean modes of variability will likely be similar to the Q RESIDUAL EOF modes.To assess any impact of the Q RESIDUAL EOF modes on atmospheric circulation, we apply lead-lag regression analysis.
Figure 9 shows the lagged regressions of the Q RESIDUAL EOFs and summarizes the sequence of events preceding and following their peak.Positive Q RESIDUAL EOF1 events are preceded by positive NAO forcing in the previous years, which cools the western subpolar North Atlantic and warms the Gulf Stream region (as indicated by Q CIRC , the unfilled contours in Figs.9k-t).These changes drive a stronger AMOC and warm SSTs along the NAC (Figs. 9a,b).The SST anomalies subsequently propagate toward the eastern subpolar North Atlantic on time scales of 4-6 years (Figs.9b-d).The time scale for the propagation of these anomalies is roughly consistent with observations (Årthun and Eldevik 2016;Årthun et al. 2017) and analysis of an earlier version of HadGEM3 (Menary et al. 2015).The Q RESIDUAL EOF2 is preceded by negative NAO-like anomalies in the 3 years prior to the Q RESIDUAL EOF2 peak.It is plausible that the anomalous heat associated with SST anomalies in the eastern subpolar North Atlantic generated through the atmospheric modulation of Q is stored in the upper ocean as the ocean mixed layer becomes shallower in spring and summer.These heat anomalies then re-emerge with the deepening of the mixed layer in the autumn and winter of subsequent years (e.g., Grist et al. 2019).
In terms of an atmospheric circulation response, following the peak of Q RESIDUAL EOF1, there is a weak Atlantic ridge response (Figs.9m,n).However, p values are greater than the 0.05 threshold for the Student's t test when serial correlation is taken into account.Similarly, there is no statistically significant regression pattern for years following Q RESIDUAL EOF2.This suggests that there is no substantial atmospheric circulation response to these modes in this model.
However, the SST patterns associated with Q RESIDUAL EOFs are considerably weaker in subsequent years than at lag 0 (Figs.9b-d,g-i) and the method removes any lag 0 response by construction.Hence, it is still possible that similar modes could excite a statistically significant atmospheric response.We investigate this further by utilizing a set of high-resSST-present atmosphere-only simulations.These are run using the same model as the piControl simulations but forced with observed SSTs and historical greenhouse gas and aerosol forcings spanning 1950-2014.The SSTs do not react to changes in circulation patterns, and hence, the direction of causality is clear.There are three ensemble members of high-resSST-present, and we take the ensemble average as the SSTs are the same in each case.We project the SST patterns associated with Q RESIDUAL EOF1 and EOF2 onto SSTs in the highresSST-present simulations over a box spanning (208-608N, 608-158W), which is shown in Figs.10a-h).We then regress SSTs, SLP, and zonal wind onto the resulting time series.It should be noted that the piControl and high-resSST-present SST patterns are similar but far from identical (Figs.10a,be,f) as the highresSST-present simulations are forced by observed SSTs, which will have different modes of SST variability and a different mean state.For example, the Q RESIDUAL EOF1 SST projection for highresSST-present is extended further southward than for the piControl (Figs.10a,b).Moreover, the Labrador Sea is substantially more prominent in the SST pattern in highresSST-present for Q RESIDUAL EOF2, while the region of cool subtropical SSTs is shifted westward relative to the piControl (Figs.10e,f).
Like for Q RESIDUAL EOF1 in the piControl, there appears to be a ridge-like response to the Q RESIDUAL EOF1 SST projection onto highresSST-present and a northward shifted jet stream (Figs.10c,d).However, this response is not statistically significant when one accounts for serial correlation.The ridge is also shifted considerably further equatorward (Fig. 10c) compared to the ridge following the peak of Q RESIDUAL EOF1 in the piControl (Figs.9m,n), possibly due to the more southward extension of the Q RESIDUAL EOF1 SST projection onto highresSST-present (Fig. 10b).The Q RESIDUAL EOF2 SST projection appears to show no substantial response (Figs.10d,h) in agreement with Figs.9r-t).Overall, it appears that the Q RESIDUAL modes diagnosed by the dynamical decomposition either do not affect atmospheric circulation or the response is weak.

Conclusions and discussion
This study has presented a method to separate the turbulent heat flux Q into a component directly related to atmospheric circulation Q CIRC and a residual component Q RESIDUAL , which is assumed to be primarily forced by ocean variability.The method uses a circulation analog technique to quantify the circulationrelated component and has been tested using the HadGEM3-GC3.1-MMpreindustrial (piControl) run for the wintertime North Atlantic.The method identifies the Labrador Sea as a region where Q variability is dominated by the atmosphere on interannual and decadal time scales, whereas Q variability in the Gulf Stream region is more ocean-driven on both of these time scales (Figs.5c,f).The Q in other regions, such as the eastern subpolar North Atlantic, shows more of a mix of oceanic and atmospheric forcings, depending on the time scale (Figs.5c,f).
The leading modes of Q RESIDUAL are characterized by a warming of the atmosphere along the Gulf Stream and North Atlantic Current (NAC) for EOF1 and a dipole of Q anomalies with cooling of the atmosphere in the western subtropical North Atlantic and warming in the eastern subpolar region for EOF2.These modes show a substantial low-frequency variability (Figs. 8c,d), and the peak of Q RESIDUAL EOF1 is closely linked with a strengthening of the AMOC (Figs. 7b-d  and 8b).The Q RESIDUAL modes are both driven, at least, in part, by atmospheric circulation anomalies in the previous year.These atmospheric circulation anomalies either alter ocean circulation patterns, as for EOF1, or contribute to SST anomalies in subsequent years via re-emergence of ocean heat anomalies.Lead-lag analysis (Figs. 9m,n,r,s) and examination of atmosphere-only experiments (Figs.10c,g) suggest that any response of atmospheric circulation to the Q RESIDUAL modes is weak.
The dynamical decomposition method requires a few choices of parameters, namely, the number of samples to fit for each year N s , the number of similar years to randomly select from N a , and the number of times that this process is repeated N r .However, the results were not found to be sensitive to these parameters and an analysis of variations in these parameters is given in appendix A. A larger uncertainty concerns the length of the time series used to calculate the decomposition as this affects the number of available analog years, which the algorithm can fit to (Fig. B1).The leading Q RESIDUAL modes are reproduced when the decomposition is applied to shorter simulations, which include variable external forcing, as discussed in appendix B. However, the order of these modes is flipped in terms of the relative variance that these explain and the specific details of the patterns differ, suggesting some sensitivity.
Although beyond the scope of this paper, it could be of interest to modify the method for application to subseasonal time scales.The Q decomposition method, for instance, could be used to examine atmosphere-ocean feedbacks related to the winter NAO (e.g., Kolstad and O'Reilly 2024).It would also be of interest to apply the Q dynamical decomposition method to other regions.For example, atmospheric circulation over the North Pacific, like the North Atlantic, is largely dominated by the internal atmospheric variability and remote forcing.Applying this method would help to identify the key modes of ocean-driven Q variability in this region.However, it may be difficult to interpret results if this method is applied to regions, such as the tropical Pacific, where the ocean strongly forces an atmospheric response.This is because the decomposition may identify the strong atmospheric response as forcing the ocean, rather than the other way round.
It should be noted that the dynamical decomposition method is primarily a diagnostic tool for separating atmospheric-and ocean-driven components of Q.That is to say, the atmosphere responds to the total Q and not Q RESIDUAL .O'Reilly et al. pacemaker experiments.Specifically, restoring tropical Atlantic SSTs toward particular patterns often results in positive Q; though in coupled runs, warm SSTs in this location are usually linked to negative Q.The dynamical decomposition method could be of particular use in comparing Q in free-running models to that in pacemaker experiments and thus establishing whether the SST restoring primarily occurs through atmospheric or oceanic adjustment.
Finally, the low amplitude of decadal variability in North Atlantic atmospheric and oceanic circulation in models compared to observations (e.g., Simpson et al. 2018) suggests model deficiencies in simulating ocean-atmosphere covariability.This dynamical decomposition could be used to diagnose differences between models and observations in the extent to which different regions of the North Atlantic are driven by the ocean and atmosphere.

Sensitivity to Period Length and External Forcing
This study has presented the results of applying the circulation analog method to a piControl run with a long time series and no externally forced variability.It would, however, be desirable to apply this to historical simulations or reanalysis datasets which are shorter and are subject to variations in greenhouse gases and aerosols.This appendix therefore investigates the sensitivity of the method to the length of the period and the presence of variable external forcing.We apply the dynamical decomposition to a fourmember ensemble of runs with observed external forcings spanning 1850-2014 (historical), which have been created using the same model.In each case, the same values of N s , N r , and N a are used as for the original piControl run.For the historical simulations, the influence of external forcing is removed by linearly regressing out the global-mean SST from all variables before performing the dynamical decomposition.To test the sensitivity to the length of the time period, we perform the dynamical decomposition using five 100-yr, nonoverlapping subsets of the piControl run.We then calculate five sets of gridpoint correlations of Q CIRC calculated from each of the separate 100-yr periods with Q CIRC diagnosed from the full 500-yr period.The average of these five sets of correlations is shown in Fig. B1a).The 500-yr period has a larger number of analogs to fit each year to which may affect the magnitude of Q attributed to atmospheric circulation variability.Broadly speaking, Q CIRC calculated from the shorter 100-yr periods is highly positively correlated with Q CIRC from the full 500-yr period Fig. B1a).There is particular agreement in regions which are dominated by atmospheric circulation variability such as the Labrador Sea (Fig. B1a).Performing the same procedure for Q RESIDUAL , the correlations are fairly uniformly around 0.6-0.7 across the North Atlantic, i.e., explaining around 40%-50% of the variance (Fig. B1b).Overall, this suggests that the period length introduces some uncertainty to the decomposition due to the available years to fit analogs to.
To investigate the role of external forcing, we perform the EOF analysis of Q RESIDUAL as before using the historical simulations.Each of these shows similar dominant patterns of Q RESIDUAL variability to the piControl, though with the order of the leading patterns flipped compared to the piControl (Fig. B2).This suggests that variable external forcing and the shorter length of the simulations do not prevent identification of the leading modes.However, it should be noted that the order of those modes is flipped in terms of variance explained and the structure of them differs between ensemble members.Overall, the similarity of these modes suggests that the Q dynamical decomposition could usefully be applied to reanalysis datasets, which have the added complexity of external forcing and are considerably shorter than a typical piControl run.

FIG. 1 .
FIG. 1. Dynamical decomposition of (a)-(c) SLP and (d)-(f) Q anomalies for the winter (DJFM) of the year 1911 in the HadGEM3-GC31-MM piControl run.(a),(d) The full fields, (b),(e) the atmospheric circulation-related components, and (c),(f) the residual components.Note that Q is defined as positive from the ocean to the atmosphere.

FIG. 3 .
FIG. 3. Interannual (DJFM-mean) gridpoint correlations between SST anomalies and (a) Q CIRC , (b) Q RESIDUAL, and (c) Q in the HadGEM3-GC1-MM piControl run.Only correlations with p values less than 0.05 are plotted following the Student's t test with the null hypothesis of zero correlation.

FIG. 5 .
FIG. 5. Variance associated with (a)-(c) interannual and (d)-(f) decadal (DJFM) variability of the components of the Q decomposition.Decadal variance is calculated as the variance of the 10-yr running mean of Q.The variance in (a),(d) Q CIRC , (b),(e) Q RESIDUAL , and (c),(f) the ratio of the variance of Q RESIDUAL to the variance in total Q.Units for (a), (b), (d), and (e) are (W m 22 ) 2 and (c) and (f) are unitless.Regions in which the total Q variance falls below 200 (W m 22 ) 2 in (c) and 20 (W m 22 ) 2 in (f) are shown by gray shading.

FIG. 7 .
FIG. 7. Lagged regressions of the Atlantic, zonal-mean overturning streamfunction (colors) onto Q RESIDUAL PC time series' associated with (a)-(d) EOF1 and (e)-(h) EOF2 as functions of depth and latitude.The mean overturning streamfunction is shown by unfilled contours, with contours drawn every 4 Sv (1 Sv ; 10 6 m 3 s 21 ), beginning at 2 Sv.Hatching indicates the statistical significance of regression coefficients, with p values below 0.05, following the Student's t test.Negative lags indicate the streamfunction state prior to the Q RESIDUAL EOFs' peaks, and positive lags indicate the streamfunction following the Q RESIDUAL EOFs' peaks.

FIG. 8 .
FIG. 8. Lagged correlations for (a) the NAO index with the AMOC and Q RESIDUAL EOF1 PC time series and (b) the AMOC with Q RESIDUAL EOF1 and EOF2 PCs.Filled circles indicate statistically significant correlations with p values below 0.05, following a t test.Normalized power spectra are also shown for the PCs associated with (c) EOF1 and (d) EOF2 of Q (black), Q CIRC (red), and Q RESIDUAL (blue).Dashed curves in (c) and (d) indicate the 95% significance levels for red noise models fitted to the Q PC time series.The colors of these curves are the same as the corresponding power spectrum they are fitted to.
FIG. A1.As in Fig. 6, but for Q RESIDUAL only and for different values of N s and N a .Shown are results for (a),(d) N s 5 50, N a 5 500, (b),(e) N s 5 30, N a 5 80, and (c),(f) N s 5 10, N a 5 80.
FIG. B1.Average of gridpoint correlations between dynamical decompositions calculated for the full 500-yr period and for this same period split into five, nonoverlapping, 100-yr periods.This is shown for (a) Q CIRC and (b) Q RESIDUAL .