Low-frequency variability enhancement of the midlatitude climate in an eddy-resolving coupled ocean–atmosphere model. Part I: anatomy

This study investigated the coupling of the wind-driven ocean gyres with the atmospheric westerly jet using an idealised, eddy-resolving, coupled model. An empirical orthogonal function analysis of the low-pass filtered data showed that the ocean gyre variability is dominated by meridional shifts of the western boundary current extension (WBCE) and changes in the strength of the subtropical inertial recirculation zone. On the other hand, the atmospheric potential vorticity (PV) variability is dominated by the growth of standing Rossby wave patterns, while its pressure variability is dominated by a zonally-asymmetric meridional shift of the atmospheric jet. Damping sea surface temperature (SST) variability in the atmosphere was shown to weaken its PV variability and reduce the zonal asymmetry of the jet-shift mode. Singular value decompositions revealed a positive feedback between meridional shifts of the WBCE and the growth of standing Rossby wave disturbances in the atmospheric jet. The atmosphere’s response is controlled by shifts in the meridional eddy heat flux over the SST front which triggers the growth of baroclinic instabilities. This instability growth eventually leads to a large-scale, barotropic pressure response over the eastern ocean basin, or an aforementioned meridional shift of the atmospheric jet. Reduction in the atmospheric resolution inhibits the ability of atmospheric eddies to resolve length scales associated with meridional shifts of the SST front and WBCE. The lack of resolution consequently weakens the influence of ocean gyre variability on the atmospheric jet and reduces the strength of the positive feedback.


Motivation
Western boundary currents (WBCs), such as the Kuroshio in the North Pacific and Gulf Stream in the North Atlantic, are fast, narrow currents that reside on the western boundary of every major ocean basin. These currents, and their associated WBC extension (WBCE) and inertial recirculation zones, show vigorous variability ranging from interannual to decadal time-scales. In addition, The convergence of warm and cold waters by the WBCs create sharp SST gradients, or a front, leading to persistent wintertime surface heat fluxes into the overlying atmosphere. Recent research suggests that resolving the so-called 'atmospheric mesoscale', on the order of 10-100 km, over these regions may lead to improvements in modelling of the midlatitude climate variability and strengthen SST-induced forcing of the upper-troposphere (Czaja et al. 2019). This claim is backed up by current observational estimates which show that the influence of WBC variability on the troposphere is considerably stronger than earlier estimates (Frankignoul et al. 2011;Taguchi et al. 2012;Révelard et al. 2018).
Although coupled general circulation models (GCMs) and observations are vital for the progressions in our understanding of the climate, idealised, quasigeostrophic, coupled models provide a unique, and often under-utilised, perspective into the midlatitude climate variability, owing to their relatively cheap cost to run and transparency under analysis. However, studies using these class of models are yet to be run under resolutions that adequately resolve the mesoscales in both the ocean and atmosphere. Furthermore, advances in 1 3 modelling of the quasigeostrophic wind-driven circulation in the past decade or so now allow for much more realistic, turbulent circulations to be achieved (e.g., Deremble et al. 2011;Berloff 2015). This warrants a significant re-investigation into the midlatitude climate variability using such models.
We will begin with a brief background on the subject, followed by a description of the models and methods used in Sect. 2. Section 3 will show the reference solutions to the benchmark modelled climate and, Sect. 4 will decompose its variability. The lagged covariability and ocean-atmosphere interactions are investigated in Sect. 5, with a final discussion of results given in Sect. 6.

Background
In the past, attributing drivers of the midlatitude coupled climate variability has proven to be complex in nature, with much of the early literature in disagreement about the role played by the midlatitude ocean circulation, even those using the same or similar models (e.g., see Barnett 1994, 1996;Schneider et al. 2002, for the classic example). It is argued that most of the midlatitude oceans act as a whitenoise integrator of high-frequency atmospheric forcing through a first-order Markov process (Hasselmann 1976;Frankignoul and Hasselmann 1977). Low-frequency sea surface temperature (SST) anomalies damp surface heat fluxes, which consequently increases the surface temperature variability of the ocean and atmosphere (Battisti et al. 1995;Barsugli and Battisti 1998;Saravanan and McWilliams 1997;Kushnir et al. 2002). For example, the basin-scale tripole SST anomalies observed in the North Atlantic and Pacific are often partly attributed to this process (Mantua et al. 1997;Okumura et al. 2001;Alexander 2010;Schneider and Fan 2012). Indeed, earlier review studies of the midlatitude coupled climate variability did not find evidence for a strong response of the extratropical atmosphere to large-scale SST anomalies (Kushnir et al. 2002;Kwon et al. 2010). In these studies, it was proposed that SST variability is excited by basin-scale modes of atmospheric variability which then weakly feedback on the atmosphere.
In recent times it has become more clear that the whitenoise integrator hypothesis is insufficient in explaining climate variability over regions of the ocean such as the WBCs and WBCEs (Smirnov et al. 2014). For starters, these regions of the ocean show persistent wintertime SST anomalies which increase the surface heat flux variability (Tanimoto 2003;Taguchi et al. 2012). Furthermore, recent reanalysis data shows the influence of this SST variability on the atmosphere is now significantly stronger than previous estimates (Frankignoul et al. 2011;Taguchi et al. 2012;Révelard et al. 2018).
The variability of WBCs and WBCEs have been suggested to be controlled by: (1) large-scale atmospheric forcings that trigger decadal responses in the ocean gyres, either through advective effects (Dewar 2003) or westwardpropagating baroclinic Rossby waves (Taguchi et al. 2007;Sasaki et al. 2013; 2) intrinsic interannual variability due to chaotic, turbulent behaviour of WBCEs (Nonaka et al. 2016;Bishop et al. 2017;Small et al. 2019;Laurindo et al. 2019). This time-scale dependence of WBC variability has been shown by Nonaka et al. (2020). These changes in strength and position of WBCEs then impact local coupling with the overlying storm track and atmospheric jet through shifts in the low-level baroclinicity (Nakamura et al. 2004;Brayshaw et al. 2008;Sampe et al. 2010;O'Reilly and Czaja 2015). The dynamics of these regions of the ocean are strongly controlled by mesoscale ocean dynamics, which must be dynamically-resolved by the model in order to fully diagnose their effects.
Furthermore, there is now research to suggest that models which resolve the so-called 'atmospheric mesoscale', on the order of 10-100 km, may also be required, rather than just resolving the Rossby deformation scales, which are on the order of 1000 km (Ma et al. 2016(Ma et al. , 2017Siqueira et al. 2021). It is postulated that as the atmospheric mesoscales are better resolved, oceanic forcing of the upper-troposphere will be enhanced (Czaja et al. 2019). For example. Famooss Paolini et al. (2022) showed that increases in atmospheric resolution leads to a westerly jet response to SST anomalies that is more similar to reanalysis data. This improved response of the downstream westerly jet was to move in the same direction as the SST front displacement. Such a response is opposite to what is expected from the linear theory (Hoskins and Karoly 1981), implying the anomalous diabatic heating associated with the SST front shift is mainly balanced by transient eddy forcing (Peng et al. 1997;Peng and Whitaker 1999). We note here that despite recent improvements in atmospheric GCMs, eddy-mean flow feedbacks still show complicated, nonlinear responses to SST anomalies (Seo et al. 2017) and increasing atmospheric resolution does not necessarily lead to circulations that are closer to observations (Czaja et al. 2019).
Within the framework of eddy-resolving, idealised ocean models driven by fixed winds, the importance of mesoscale eddies in ocean gyre variability is well-known through detailed analyses of the underlying potential vorticity (PV) dynamics (e.g., Berloff and McWilliams 1999a;Hogg et al. 2005;Berloff et al. 2007b;Deremble et al. 2011;Berloff 2015;Shevchenko et al. 2016;Kurashina et al. 2021). Furthermore, the importance of mesoscale ocean eddies in the midlatitude climate variability has been shown through studies using idealised coupled models. For example, Dewar (2001) showed that mesoscale eddy variability feeds back strongly on ocean Ekman pumping, peaking at decadal time-scales, and thus controls the SST response and atmospheric variability. Time-scale control of atmospheric variability has also been shown to be important by Hogg et al. (2006) and Berloff et al. (2007a), with preferred decadal time-scales selected by nonlinear adjustment of the ocean gyres (Dewar 2003). This is able to set the consequent timescale of decadal coupled oscillations through transitions between pairs of preferred atmospheric states (Kravtsov et al. 2007). For example, Kravtsov et al. (2007) found that as the westerly atmospheric jet switches from a poleward to an equatorward state, nonlinear adjustments of the WBCE to anomalous forcing lead to a persistent SST anomaly forming. This SST anomaly anchors the atmosphere to remain in its equatorward latitude state through strong surface heat fluxes until the anomaly dissipates. After the SST anomaly dissipates, the atmosphere returns to its more preferred poleward state, forcing the WBCE back poleward. Suppressing mesoscale ocean eddies had a weakening effect on the coupled variability, further supporting the view that mesoscale ocean eddies are vital in not only in controlling ocean variability but also the overlying atmosphere. Furthermore, Kravtsov et al. (2006) and Farneti (2007) showed a westward-propagating, coupled, phase-locked Rossby wave which acted on an interannual time-scale associated with the basin crossing time. This coupled interaction was shown to not require the atmosphere to transition between two preferred states (Kravtsov et al. 2006) and was akin to the positive feedback mechanism proposed by Goodman and Marshall (1999). However, an important point to note is that none of these studies have modelled climates which adequately resolve both the ocean and atmosphere mesoscales. Furthermore, none have investigated the influence of atmospheric resolution on the modelled climate, nor have they detected modifications in atmospheric modes of variability by the ocean circulation. This paper will re-establish the anatomy of the low-frequency climate variability in a quasigeostrophic coupled model run under much higher resolutions than those in previous studies (Dewar 2001;Hogg et al. 2006;Kravtsov et al. 2006Kravtsov et al. , 2007Berloff et al. 2007a). Our main aim in this study is to then see whether any novel coupled behaviours are observed and attempt to explain the mechanisms involved.

The Quasi-Geostrophic Coupled Model (Q-GCM)
The model we will be using for this study is version 1.5.0 of the Quasi-Geostrophic Coupled Model (Q-GCM) which was first described and implemented in Hogg et al. (2003). It is designed to mimic aspects of the midlatitude climate system and consists of a wind-driven, double-gyre box ocean coupled to a periodic, channel atmosphere through mixed layers that allow for transfers of heat and momentum. The model is similar to the one used in Kravtsov and Robertson (2002) and Kravtsov et al. (2007) but also includes a dynamically-active atmospheric mixed layer.

Notation
Before we move into describing the model, we will introduce the notation used in this paper which is identical to the one used in Hogg et al. (2003). Left-superscript, given by 'o' or 'a', indicates ocean or atmosphere variables and rightsubscript, given by 'm' or integer i, indicates the mixed layer or i-th isopycnal layer, respectively. Layers are counted away from the ocean-atmosphere interface.

Model Geometry
The box ocean model consists of 3 isopycnal layers with depths o H 1 = 350 m, o H 2 = 750 m and o H 3 = 2900 m . The benchmark solution has a nominal ocean resolution of 5 km for a basin with lateral dimensions 4800 km × 4800 km . This gives a discretised grid of size 961 × 961 points representing cell vertices in the longitudinal and latitudinal directions, respectively. Within the upper-isopycnal ocean, there is a mixed layer of fixed depth o H m = 100 m . The channel atmosphere also consists of 3 isopycnal layers with depths a H 1 = 2000 m, a H 2 = 3000 m and a H 3 = 4000 m . The benchmark solution has a nominal atmosphere resolution of 20 km for a channel with lateral dimensions 30 720 km × 7680 km . This gives a discretised grid of size 1537 × 385 points representing cell vertices in the longitudinal and latitudinal directions, respectively. Within the lower-isopycnal atmosphere, there is also another mixed layer of variable depth with an unperturbed thickness of a H m = 1000 m.
Parameters in Tables 1 and 2, such as eddy viscosity and boundary condition coefficients, have been chosen such that a realistic and turbulent circulation is achieved. For example, it is well known that different choices of these parameters have large and profound effects on the ocean circulation (e.g., Haidvogel et al. 1992;Berloff and McWilliams 1999a, b;Hogg et al. 2005;Deremble et al. 2011), so care has been taken to pick the correct parameter choices.

The quasigeostrophic (QG) layers
The QG equations for the ocean and atmosphere may be written in terms of PV anomalies and dynamic pressure anomalies. For example in the ocean, they are denoted in vector form as o q = ( o q 1 , o q 2 , o q 3 ) T and o p = ( o p 1 , o p 2 , o p 3 ) T , respectively. Ocean PV anomalies are defined as where f 0 = 9.37 × 10 −5 s −1 is the Coriolis parameter and A is a matrix which defines interactions between layers through dynamic pressure anomalies or buoyancies with g ′ i as the reduced gravities. In the ocean, g � 1 = 0.0222 ms −2 , g � 2 = 0.0169 ms −2 which correspond to Rossby deformation radii of 40.0 km and 20.6 km , respectively. The definition of PV anomaly in the atmosphere is the same as the ocean apart from the addition of land orography which is given in Appendix 1.
We use (x, y) as our zonal and meridional coordinates where x is always measured increasing eastwards from the western boundary of the model domain; y is measured increasing polewards from the central latitude y 0 .
Dynamic pressures in the i-th isopycnal ocean may be converted to the transport streamfunction through the relation Zonal and meridional velocities are computed through taking derivatives of o p i Atmospheric velocities may be computed in an identical manner.
We may now write the quasigeostrophic equations for the ocean as where o A 2 = 50 m 2 s −1 and o A 4 = 2.0 × 10 9 m 4 s −1 are the Laplacian and biharmonic diffusion coefficients, respectively. The Jacobian operator is given by J(f, g) = f x g y − f y g x and the horizontal Laplacian and biharmonic operators are given by ∇ 4 H and ∇ 6 H , respectively. The forcing in the QG layers is defined by the 3 × 4 matrix B and the 4 × 1 entrainment vector o e . The QG equations for the atmosphere are similar except for the aforementioned forcing terms which we will now describe.

Entrainment vectors
In the ocean, the matrix B is given by In the atmosphere, we simply have a B = − o B . The entrainment vectors o e and a e for the ocean and atmosphere, respectively, are given by  These forcing vectors correspond to vertical velocities into each QG layer. Within the upper-isopycnal ocean layer, we combine the adiabatic and diabatic forcing as a single term to simplify our subsequent analysis. We call this term the upper-ocean total entrainment 1 Ekman pumpings in the ocean and atmosphere are computed as the curl of surface stress (dropping left-superscripts) The ocean and atmosphere surface-stresses are computed as a linear combination of lower-isopycnal atmosphere velocities (Hogg et al. 2014). This means that ocean Ekman pumping depends only on the atmospheric flow and does not explicitly depend on the ocean circulation. 2 The atmosphere diabatic entrainment a e 1 in (7) contains radiative heat flux terms that occur due to changes in the atmospheric mixed layer. When an unstable stratification occurs, the entrainment a e 1 is adjusted with a correction term (Hogg et al. 2014).

Boundary conditions
Partial-slip boundary conditions are applied on the lateral boundaries (Haidvogel et al. 1992) by where bc is the non-dimensional partial-slip boundary condition parameter (see Tables 1, 2) and Δ o x is the horizontal grid spacing. The subscript n denotes outward normal partial derivatives. We choose a partial-slip boundary condition parameter that is close to free-slip as this gives a more realistic representation of the nonlinear dynamics in the WBCs, see Deremble et al. (2011) for details.
There are also mass and momentum conservation constraints applied to the lateral boundaries (McWilliams 1977;Hogg et al. 2014).

Radiative heat fluxes
The entire model is driven by solar radiation which is fixed in time and depends on latitude. A linearised radiation scheme then redistributes heat throughout the atmosphere (Hogg et al. 2014).

Sensible and latent heat fluxes
Sensible and latent heat fluxes are parametrised by a single term where = 35 Wm −2 is the sensible and latent heat flux coefficient.

The ocean surface mixed layer
The SST evolution in the ocean surface mixed layer, denoted by o T m , is given by We also have o T 1 , o T 2 , o T 3 which are fixed, constant temperature of the isopycnal layers. We will present T m as a deviation from a temperature constant which is obtained through radiative balance (Hogg et al. 2014). For our model parameters, this constant temperature for the ocean mixed layer is 27.0 • C . Note that (12) also contains Laplacian and biharmonic diffusivities given by o K 2 and o K 4 , respectively (see Table 1). The forcing terms in (12) contain thermal forcing at the ocean surface a F 0 and entrainment heat flux out of the bottom of the ocean mixed layer o F e+ m .

The atmosphere surface mixed layer
The temperature evolution of the atmospheric mixed layer is given by a similar equation to (12) except the forcing terms are slightly different (Hogg et al. 2014). Heat flux at the top of the atmospheric mixed layer is given by a F m and surface forcing given by a F 0 .
The combining of the adiabatic and diabatic forcings is justified because, as it turns out in Sect. 4.4, the leading modes of variability of both components are controlled by the same modes of atmospheric variability. 2 Q-GCM does allow for ocean velocities to be included into the computation of surface stresses but we have chosen not to do so in order to simplify mechanical atmosphere-ocean interactions.

Partially-coupled climate
In order to attribute features of the modelled climate and pinpoint sources of variability to ocean-induced or coupled effects, we run a partially-coupled modelled climate alongside our fully-coupled climate simulation. This is done by changing a parameter X C which controls the coupling between the atmospheric mixed layer and the lowerisopycnal atmosphere. When X C = 1 , the model is fullycoupled. On the other hand, when X C = 0 , the atmosphere is decoupled from atmospheric mixed layer variability (see Appendix 2 for details). This will be our parameter choice for our partially-coupled climate. This coupling parameter is useful as any SST-driven variability will be suppressed while purely atmosphere-driven variability will remain. Note that mechanical interactions between the ocean and atmosphere are implemented identically to the fully coupled case.

The inter-gyre boundary
An inter-gyre boundary that partitions the ocean basin is computed by following the time-averaged pressure contour emanating from the western boundary. This approximates the position of the WBCE as it separates from the western boundary and extends into the open ocean. This inter-gyre boundary will be used to highlight the location of vorticity fluxes into the ocean gyres which arise as a result of the atmospheric variability.

Filtering data
Before any statistical techniques are applied to the data, a 2-year low-pass filter is applied to remove high frequency oscillations. This allows for any low-frequency signal to be picked up more easily in the subsequent analysis.

Empirical orthogonal functions (EOFs)
An empirical orthogonal function (EOF) analysis is performed by computing an eigendecomposition of a covariance matrix obtained from the oceanic/atmospheric state variables (von Storch and Zwiers 1999). This gives from most to least important, the patterns of variability (EOFs) that explain the most variance in the data. We will present the explained variance by each EOF as a percentage of the total variance. The filtered data set is then projected onto the EOFs to obtain the principal components (PCs), which give information about the temporal evolution of the EOFs. Analysing the power spectra will also indicate the most important time-scales at play.

Lagged singular value decompositions (SVDs)
Although EOFs are useful in decomposing single field variability, they lack information regarding modes of covariability in the modelled climate, e.g. ocean-atmosphere interactions. In order to pick out these modes of covariability, a singular value decomposition (SVD) analysis is performed on a lagged cross-covariance matrix constructed using two state variables (Bretherton et al. 1992). Similar to the EOF decomposition, the data is first 2-year low-pass filtered and both fields are subsampled onto an 80 km grid. This is to save on memory cost as cross-covariances matrices require considerably larger storage than covariance matrices. By maximising squared covariance fraction with respect to time-lag, the likely time-scale for any delayed response may also be obtained, if one exists. Analogous to the EOF decomposition, a SVD produces patterns that maximise the cross-covariance in the data. Note that as we are concerned with pairs of fields, two sets of patterns are produced in the decomposition-one for each field. These are called the singular vectors or SVD modes. Again, by projecting the original data sets onto these modes, we obtain two sets of temporal coefficients. These may be correlated to obtain a measure of the strength of the coupling between the pair of singular vectors. A description of the required SVDs will be given in Sect. 5.

Modelled climate
Reference solutions for the benchmark and partially-coupled modelled climates are computed for a 120-year simulation length with a 20-year spin-up period and data saved every 5 days. Full list of parameters are presented in Tables 1 and 2. The simulation was started from rest and in radiative balance. Note that barred variables • denote time-averaged variables.

Ocean circulation
The modelled wind-driven, double-gyre ocean circulation is presented in Fig. 1. The upper-isopycnal ocean is spun up by both adiabatic Ekman pumping and diabatic entrainment which creates the subpolar and subtropical gyres (see Fig. 3). Western intensification leads to the formation of powerful WBCs and downstream inertial recirculation zones and WBCE. The basin shift and addition of land orography has created a larger, but weaker, subtropical gyre and has deflected the jet-axis tilt slightly poleward. The jet-axis tilt is determined by the relative strengths of the recirculation zones with the stronger recirculation zone tending to pull the jet towards itself (Moro 1988;Hogg et al. 2005). The middle-isopycnal ocean is forced only by the diabatic entrainment term (7) which creates gyres of weaker circulation with PV of the opposite sign.

Atmosphere circulation
The modelled midlatitude, channel atmosphere circulation may be seen by looking at Fig. 2. Latitudinally varying solar radiation generates a PV gradient in the lower and middleisopycnal atmosphere layers, which leads to the formation of a highly turbulent, zonally inhomogeneous westerly jet. The zonal inhomogeneity of the jet is created by the ocean basin shift and addition of land orography. The westerly jet formed in the middle-isopycnal atmosphere dominates the circulation with its momentum supporting the lower-isopycnal layer via baroclinic instabilities (Hogg et al. 2003). The upper-isopycnal atmosphere is not forced directly by diabatic entrainments, and is only driven by eddies, so we will neglect its dynamics for the purposes of our study. Average speeds in the middle-isopycnal atmosphere are around 20 ms −1 , while they are around 6 ms −1 in the lower-isopycnal atmosphere. This difference in circulation strengths is responsible for the sign change in PV in the lower-isopycnal layer. The circulation in the middle-isopycnal atmosphere dominates despite receiving weaker forcing than the lower-isopycnal atmosphere because it is not spun down by frictional wind-stresses, and the effect of considerably lower reduced gravities at the layer-2,3 interface allows for enhanced PV generation via the corresponding interface displacements that occur.

Ocean mixed layer
Instantaneous and time-averaged fields of SST anomalies, wind-induced ocean forcing and surface stresses are shown in Figs. 3 and 14. The subtropical and subpolar gyres in the ocean mixed layer are clearly visible and separated by a sharp SST front across the WBCE (Fig. 3a, e). The WBCE acts as a partial inter-gyre barrier, separating warmer waters in the subtropical gyre from the cooler waters in the subpolar gyre. Inter-gyre exchanges of warm and cool waters by mesoscale activity is also visible. Surface stresses generate wind-curls that drive the ocean gyres in the form of both adiabatic Ekman pumping and diabatic entrainments. Although the sum of these forcings, i.e. the upper-ocean total entrainment o e total , has roughly equal vorticity flux into each gyre, the adiabatic Ekman pumping component dominates in the subpolar gyre (see Fig. 3f, h). These positive and negative vorticity fluxes over the subpolar and subtropical gyres, respectively, are fuelled by the atmosphere through frictional wind-stresses. However, the sign of PV anomalies in the lower-isopycnal layer is inconsistent with the sign of Ekman pumping since negative PV anomalies sit at lower latitudes and increase as they move poleward 3 (see Fig. 2b). Hence, the PV driving the ocean gyres must originate from the middle-isopycnal layer which is consistent with our findings so far (and by Hogg et al. 2003) that the lowerisopycnal atmosphere is driven from above through momentum transfers that take place via baroclinic instabilities.

Atmosphere mixed layer
and 3 a X∕4 ≤ x ≤ a X are omitted where ocean-atmosphere interactions are assumed to be weak. This is repeated for subsequent figures unless stated otherwise. a, e Instantaneous and time-averaged lowerisopycnal atmosphere dynamic pressures. b, f Instantaneous and time-averaged lower-isopycnal atmosphere PV anomalies. c, g Instantaneous and time-averaged middle-isopycnal atmosphere dynamic pressures. d, h Instantaneous and time-averaged middle-isopycnal atmosphere PV anomalies the atmospheric circulation except over the ocean basin where the WBCs and double-gyre circulation leave an imprint through the surface heat fluxes. Upward transfers of heat are given by positive heat fluxes, i.e. ocean warming atmosphere, and downward transfers of heat are given by negative heat fluxes, i.e. atmosphere warming ocean. We confirm that the sensible and latent heat fluxes give realistic patterns (Shaman et al. 2010;Kwon et al. 2010) which are strongest near the WBC separation point and WBCE, and hence, strongly controlled by the ocean gyre circulation. The time-average effect of these heat fluxes appear to be the warming and cooling of ASTs over the subtropical and subpolar WBCs, respectively. This consequently leads to a sharpening in meridional AST gradients over the ocean basin through the convergence of warm and cold waters by WBCs. Similar patterns of surface heat fluxes are obtained in Kravtsov et al. (2007) but ours show higher variability and time-averaged fluxes that extend well into the ocean interior due to the powerful WBCE. The position of these strong sensible and latent heat fluxes over the WBCs and WBCE coincides with the narrowing of meridional eddy heat fluxes in Fig. 4g. This eddy heat flux implies the convergence of ASTs due to eddies in the atmosphere mixed layer. Meridional eddy heat fluxes are active over the entire atmospheric jet, but it acts over a wider latitudinal band outside the ocean basin region, where synoptic scale turbulence destroys the low-level baroclinicity. The narrowing of this eddy heat flux latitudinal band by sharp SST gradients due to the WBCs is important in anchoring the westerly jet and increasing the low-level baroclinicity of the atmosphere (Nakamura et al. 2004(Nakamura et al. , 2008. Enhancement of time-averaged atmosphere diabatic entrainment due to increased land-sea contrast and the addition of land orography is clearly visible in Fig. 4d. Our modelled climate is qualitatively similar to those modelled by Hogg et al. (2006) and Kravtsov et al. (2007). However, some differences such as significantly increased turbulent activity in both the ocean and atmosphere are important to note. In the ocean, this leads to a strengthened WBCE via enhanced eddy backscatter in the inertial recirculation zones. This allows for greater amounts of heat transport by the jet as it is able to penetrate further into the ocean interior. This gives a greater longitudinal region for the ocean to thermally feedback onto the atmosphere, where previously, only regions close to the western boundary were able to supply large enough heat fluxes to the atmosphere such as to measurably impact upon its dynamics. Furthermore, the effects of eddies in the lower-isopycnal behaviour must also not be ignored. The convergence of heat by meridional heat fluxes over the WBCs and WBCE plays an important role in restoring the low-level baroclinicity. Mesoscale SST variability also appears well captured by the atmospheric mixed layer due to the high nominal resolution of the atmosphere.

Variability of modelled climate
We will now decompose the variability of the low-pass filtered ocean and atmosphere data. Furthermore, we must look at the variability of mixed layer variables that are responsible for mediating ocean-atmosphere interactions. These include forcings such as wind-induced Ekman pumping, ocean and atmosphere diabatic entrainments, as well as transfers of heat and momentum through wind-stresses and surface heat fluxes. Although similar analysis is performed in Hogg et al. (2005Hogg et al. ( , 2006, our modelled ocean and atmosphere circulation is significantly more turbulent so this analysis is necessary for the purposes of this study.

Upper-isopycnal ocean variability
The leading EOFs for the filtered upper-isopycnal ocean PV anomaly o q 1 from the benchmark modelled climate are given in Fig. 5. EOFs 1 and 2 ( Fig. 5a, b) explain 24.5% and 18.1% of the total variance in the data, respectively, with strong power at low frequencies (see Fig. 5c, d). EOF 1 consists of strong negative PV anomalies situated along the meandering WBCE. EOF 2 consists of similar, but more zonallysymmetric, negative PV anomalies along the WBCE. EOF 1 also appears to show a change in the subtropical gyre PV concentration, and EOF 2 very weakly so, but these are secondary to the strong variability of the WBCE. Comparing the power spectra of the corresponding PCs from the decoupled simulation shows that coupling of the atmosphere to its mixed layer has a strong effect on the ocean gyre variability (see Fig. 5c, d). Both EOFs have increased variability across almost all low frequencies through coupling. Such changes in variability are likely a consequence of thermal feedbacks between the ocean and atmosphere but it is unclear why this is so without further analysis of the lagged cross-covariability (see Sect. 5).
Although it is evident that these EOFs control WBCE variability, it is less clear how they affect the gyres by looking at the EOFs alone. In such situations it is better to overlay the EOFs onto the mean flow to show positive and negative phases of an oscillation generated by each EOF. However, we found that this was also not very insightful since the yielded PV anomaly fields still looked very similar. This is largely due to the fact that the PV anomaly fields are more homogenised throughout the ocean basin in comparison to streamfunction fields which are considerably stronger around the WBCE and recirculation zones. Thus, it is more useful to look at the leading EOFs of filtered upper-isopycnal transport streamfunction o Ψ 1 which gave EOFs that were highly correlated to the EOFs obtained from the filtered upper-isopycnal PV anomaly data. Positive and negative phases are defined by a ±1 oscillation on top of the time-averaged circulation where represents 1 s.d. of the corresponding PC.
We found that o q 1 EOF 1 was correlated to o Ψ 1 EOF 2 with correlation coefficient 0.85 and o q 1 EOF 2 was correlated to o Ψ 1 EOF 1 with correlation coefficient of also 0.85. The positive and negative phases of the upper-isopycnal ocean transport streamfunction EOFs are presented in Fig. 6 along with the original EOFs and explained variances. From this figure, it is now apparent that o q 1 EOF 1 affects the strength of the subtropical recirculation zone which coincides with a change in the subtropical gyre PV concentration. Changes in the subpolar gyre and recirculation zone appear very small or negligible. In the positive phase, the subtropical recirculation zone strengthens, leading to an increase in the jet meander as the jet is consequently deflected more strongly equatorward as it separates from the western boundary, and the size of the subtropical gyre shrinks. In the negative phase, the opposite happens. The subtropical recirculation zone weakens, the jet is deflected less as it separates from the western boundary, and the size of the subtropical gyre increases. The maximum volume transports achieved by the positive and negative phases for o Ψ 1 EOF 1 are 70 Sv and 65 Sv , respectively.
On the other hand, o q 1 EOF 2 (or o Ψ 1 EOF 1) appears to affect the meridional displacement of the WBCE, with the entire jet shifting poleward in its positive phase and equatorward in its negative phase (see Fig. 6c, d, e). The magnitude of this meridional displacement appears to be ∼ 100 km from its mean position and is similar to observed values (Sasaki et al. 2013). There is also a smaller but measurable change in the strength of the subtropical gyre, recirculation zones and WBCE in the different phases of o q 1 EOF 2 than what is observed in o q 1 EOF 1. The maximum volume transports achieved by the positive and negative phases for o Ψ EOF 2 are slightly less than EOF 1 with values of 67 Sv and 63 Sv , respectively, for phases at ±1 . Sasaki et al. (2013) also found this increase in jet-strength as the WBCE moved poleward. These are intrinsic modes of ocean gyre variability, with similar modes observed by Hogg et al. (2005), Hogg et al. (2006) and Berloff et al. (2007b). Our interests in subsequent sections lie in understanding how these modes of variability may couple with the overlying atmosphere through lagged responses and coupling.

Ocean mixed layer variability
Analysis of SST o T m and surface heat flux a F variability found similar EOFs to those seen in Fig. 5 and are presented in Fig. 6a, b, f, g. These EOFs are highly correlated to the EOFs in Fig. 5 and are similar, but not identical, in structure to each other. The main difference being that a F EOF 2 does not have the corresponding negative heat flux due to cooling of the subtropical gyre in o T m EOF 2 (see Fig. 6f, g). This is likely an effect of damped surface heat fluxes in the ocean interior due to the white-noise integrator effect (Hasselmann 1976;Frankignoul and Hasselmann 1977). The top row of panels in Fig. 6, o T m EOF 2, a F EOF 2 and o Ψ 1 EOF 1, are correlated to o q 1 EOF 2; the bottom row of panels in Fig. 6, o T m EOF 1, a F EOF 1 and o Ψ 1 EOF 2, are correlated to o q 1 EOF 1. All of the respective correlation coefficients are greater than 0.8 indicating that all the EOFs correspond to the same changes in the underlying ocean gyre circulation. Similar patterns due to o Ψ 1 EOF phase changes were found c Power spectra of PC 1 for coupled (X C = 1) and decoupled (X C = 0) simulation. d Power spectra of PC 2 for coupled (X C = 1) and decoupled (X C = 0) simulation for a F phase changes. It appears that in the positive phases of a F EOFs 1 and 2, the total sensible and latent heat flux due to the subtropical WBCs increase by 4% and 2%, respectively. In the negative phase, the sign of the anomaly flips and the total sensible and latent heat fluxes decrease by the same amount.

Atmosphere and diabatic entrainment variability
The leading EOFs of the filtered atmospheric variability, namely lower-and middle-isopycnal atmosphere PV and dynamic pressures, as well as their associated diabatic entrainments are presented in Fig. 7. These EOFs, which are all highly correlated to each other, 4 consist of a wavenumber-6 mode that affects the meander and zonal symmetry of the atmospheric jet with an east-west dipole anomaly sitting over the ocean basin. This mode is a synoptic scale Rossby wave that forms in the atmosphere through growing instabilities in the westerly jet. Similar modes of variability were found in Hogg et al. (2006) with a pair of EOFs identified as a travelling-wave mode pair. We also found similar modes shifted out of phase in the longitudinal direction but these were omitted since they appeared mostly unaffected by coupling of the atmosphere to its mixed layer and had significantly lower explained variances. For example, a q 1 EOF 1 (Fig. 7b) has an explained variance of 21.0% in the fully-coupled modelled climate (X C = 1) , but only has an explained variance of 10.8% in the equivalent EOF from the partially-coupled modelled climate (X C = 0) . This doubling of the explained variance is apparent in the power spectra of the respective PCs in Fig. 7g with increases in variability seen over the 5-16-year band range (increased variability over longer time-scales > 20-year are less certain). The increase in variability of this frequency band indicates the most dominant time-scales that are likely associated with ocean-atmosphere interactions. Similar changes in explained variance and power spectra were found for the corresponding EOFs in the middle-isopycnal layer (see Fig. 7h) when comparing the coupled and partially-coupled modelled climates.
Yet, no such change in explained variance was observed for the omitted phase-shifted EOFs. The question now arises why the variability of one EOF is amplified while the other is not. It is natural to expect that the time-varying diabatic entrainment would increase the atmospheric variability of these two leading EOFs with roughly equal measure as their structures are very similar. The fact that the EOFs in Fig. 7b, c, e, f alone appear affected is indicative of the influence of ocean gyre dynamics on this particular mode of variability. Due to these reasons, we will treat the EOFs in Fig. 7 as separate EOFs from their phase-shifted counterparts. The corresponding diabatic entrainment mode, or atmospheric forcing, that generates these modes of atmospheric variability are presented in Fig. 7a, d. In fact, they are the same EOF but sign-flipped since the lower-and middle-isopycnal layers receive proportional but opposite-signed forcings (7). These diabatic entrainments generate the strongest circulation response in the middle-isopycnal layer, with the PV response pattern advected eastwards of the diabatic entrainment forcing by the jet before reaching maximum strength. The corresponding dynamic pressure anomalies in the middle-isopycnal layer are of the correct sign with cyclonic circulations for positive PV anomalies and anticyclonic circulations for negative PV anomalies. Momentum transfers into the lower-isopycnal layer via baroclinic instabilities flip the sign of the circulation anomalies in the lower-isopycnal, despite receiving forcing of the opposite sign, such that the flow becomes more barotropic. The atmosphere diabatic entrainment EOF 1 is strengthened and more coherent over the ocean basin with enhanced forcing first appearing over the western boundary near the WBCE separation region. These then reach a maximum strength over the eastern boundary before they quickly decay. We Fig. 7 Leading EOFs of atmosphere diabatic entrainment a e 1 and filtered lower-and middle-isopycnal atmosphere PV. a a e 1 EOF 1. b a q 1 EOF 1. c a p 1 EOF 2. d − a e 1 EOF 1. e a q 2 EOF 1. f a p 2 EOF 1. g Corresponding power spectra of EOF in panel b for coupled (X C = 1) and decoupled (X C = 0) simulations. h Corresponding power spectra of EOF in panel e for coupled (X C = 1) and decoupled (X C = 0) simulations interpret this as the generation and growth of baroclinic instabilities over the ocean basin due to land-sea contrast and strong surface heat fluxes generated by the presence of WBCs and SST front (see Fig. 4f, g). The decay of this mode to the east of the ocean basin is associated with strong eddy mixing of heat by turbulence of the jet that destroys the increased baroclinicity created by the land-sea contrast. This tendency of the atmosphere to reduce meridional temperature gradients and low-level baroclinicity has been posed by Hoskins and Valdes (1990). Although it is clear that the atmospheric westerly jet is interacting with the SST front and WBCs, it is unclear through this analysis alone what role the SST variability plays. In summary, the strongest mode of atmospheric variability is controlled by the growth of Rossby wave modes through baroclinic instabilities that are enhanced by the presence of the SST front.
When the low-pass filtered, lower-isopycnal atmosphere pressure variability a p 1 was decomposed in our model, the leading EOF was a standing wave, meridional jet shift mode (see Fig. 8a) explaining a large 37.3% of the total variability. This mode is similar to the jet-shift mode found in Kravtsov and Robertson (2002) and Kravtsov et al. (2003) but with weak circulation anomalies sitting in-between strong circulation anomalies. These weaker circulation anomalies are regions that are largely unaffected by the anti-cyclonic circulations and increase the zonal asymmetry of the mode. Over the ocean basin, this mode corresponds to a strong anti-cyclonic circulation anomaly over the eastern basin and a much weaker, almost stagnant, anti-cyclonic anomaly over the western basin. The differences in structure with the equivalent EOF obtained from the partially-coupled modelled climate are notable with a more zonal structure and no such weak circulation anomalies, i.e. more similar to modes found by Kravtsov and Robertson (2002) and Kravtsov et al. (2003). Indeed the EOF in Fig. 8a appears similar in structure to a q 1 EOF 1 in Fig. 7b as it is of the same wavenumber. Since ocean gyre dynamics has already been posed as the likely culprit in increasing the variability of this particular EOF, it is also likely that it is responsible for the changes Fig. 8 Leading EOF of filtered lower-layer atmospheric dynamics pressure. a: a p 1 EOF 1 for benchmark modelled climate. b a p 1 EOF 1 for partially-coupled modelled climate. c Relative power spectra of the respective PCs in structure between Fig. 8a, b. 5 Such a change in zonal structure an atmospheric mode of variability through coupling has yet to be seen in such models with Hogg et al. (2006), Kravtsov et al. (2007) and Berloff et al. (2007a) only being able to find a modification of the preferred time-scales, rather than a change in structure of the EOF itself.

Oceanic forcing variability
Finally, we decompose the variability of the wind-induced ocean forcing. Namely we look at ocean Ekman pumping o w ek , ocean diabatic entrainment o e 1 , and in the case of the upper-isopycnal ocean, the upper-ocean total entrainment o e total . This is left until after the atmospheric variability is decomposed as it turns out that ocean forcing variability is controlled by the atmosphere at lag-zero. The leading EOFs and their associated ocean surface stresses are presented in Fig. 9. Unlike the variability of the diabatic atmosphere entrainment, which is correlated at lag-zero to the modes of lower atmosphere PV variability, ocean forcing does not correlate at lag-zero to modes of ocean gyre variability due to its large inertia. The leading EOF of the filtered ocean Ekman pumping, ocean diabatic entrainment and upper-ocean total entrainment explain 11-12% of the total variance (see Fig. 9a-c). These modes looks similar and are correlated (correlation coefficients > 0.7 ) to the leading EOF of a q 1 , a q 2 in Fig. 7b, e. These EOFs consists of an east-west positive-negative dipole vorticity flux over the ocean basin. The positive anomaly is situated over the WBCE and the negative anomaly is situated over the eastern basin. The centroid of both anomalies shift north and south of the inter-gyre boundary depending on whether we are looking at Ekman pumping (Fig. 9a) or diabatic entrainment (see Fig. 9b). 6 The shifting of these anomalies is due to the aforementioned misalignment of the SST front and the zero wind-curl line (see Sect. 3.3). The o e total EOFs are advected slightly further downstream of the atmospheric westerly jet compared to the a q 1 EOFs which is due to Ekman layer effects when computing surface stresses in the mixed layers. The signs of the anomalies is such that losses of PV in the lower-isopycnal atmosphere correspond to gains of PV in the ocean. However, the responsible PV flux must be generated by the middle-isopycnal atmosphere where the signs are correct, i.e. positive vorticity fluxes in the ocean must be due to positive PV anomalies that arise in the atmosphere otherwise enstrophy is not conserved. This pathway where circulation anomalies in the middle-isopycnal layer transfer momentum to the lower-isopycnal atmosphere and into the ocean is again consistent with our earlier findings (see Sect. 3.2). The responsible surface stress patterns (Fig. 9d, e) are EOF 2 of zonal surface stress o x and EOF 1 of meridional surface stress o y . The wind-stress anomalies give cyclonic surface stresses for positive vorticity fluxes and anti-cyclonic surface stresses for negative vorticity flux as expected. The leading EOF of zonal wind-stress, o x EOF 1 consists of a zonal shear pattern (omitted). The zonal shear appears strongest over the eastern basin with the sign of the shear corresponding to a cyclonic wind-stress pattern.
Analysis of ocean forcing variability has yet to be looked at in detail using this model. This is likely because ocean gyre variability has been, so far, only explained as intrinsic modes of variability (e.g. Hogg et al. 2006) or within oceanonly simulations where the atmosphere model is switched off (Hogg et al. 2005). However, in order to be able to understand the role of the ocean gyres within any coupled interaction. Modes of atmospheric and ocean forcing variability must be related to delayed modes of ocean gyre variability. 7 Indeed, Kravtsov et al. (2007) identified monopolar Ekman pumping forcing anomalies over the ocean basin to be part of a coupled mode associated with meridional shifts of the atmosphere westerly jet preceding shifts in the WBCE through nonlinear adjustments of the recirculation zones (Dewar 2003). The structure of the ocean forcing EOFs in this study are different since there are zonal asymmetries present which make it more difficult to diagnose its effects on the gyres. The zonal structure of these forcing modes is more similar to those discussed in Jin (1997). However, the ocean model used in their study consisted of a linear Rossby wave model which does not take into account the complex nonlinear dynamics present in the ocean gyres.

Disentangling causes of climate variability
Although we have made some progress in understanding the modelled climate variability, we still have some work to do in understanding how certain aspects of it interact. Namely, we are most interested in the ocean gyre response to forcing as well as how the associated surface heat flux anomaly affects the atmosphere diabatic entrainment. The ocean gyre response to forcing will be looked at by computing SVDs of the relevant data fields under different time-lags, with forcing leading response variables.

Ocean gyre response to wind-induced forcing
When the leading SVD modes were computed for upperocean total entrainment o e total leading upper-isopycnal ocean PV anomaly o q 1 , two different leading SVD modes were found for a wide band of time-lags ranging from 0 to 16-years (see Fig. 10a, b, d, e). Longer time-lags were found to decorrelate with known EOFs and became too noisy to interpret. The two sets of SVD modes were aggregated for different time-lags by averaging over them with equal weights, i.e. computing their ensemble average. The o e total component of both of these SVD modes (Fig. 10a, b) were found to have average correlation coefficients of 0.72 and 0.78 with o e total EOF 1, respectively. Indeed both of these patterns are spatially correlated with the corresponding EOF (positive anomaly in the western basin and negative anomaly in the eastern basin) but the SVD mode appears both noisy and sheared in the anti-cyclonic direction. This shearing of the forcing pattern is either part of the noise or an important feature of the signal such as an effect of coupling. 8 Although the responsible entrainment EOF is the same for the pair of presented SVD modes, where they differ is in the upperisopycnal PV anomaly response. In fact, the upper-isopycnal ocean gyre response of every computed leading SVD mode, except the 12-year lag, was found to be either one of the following two responses. The first mode, shown in Fig. 10d, is anti-correlated to o q 1 EOF 1 with a correlation coefficient of −0.96 while the second mode, shown in Fig. 10e, is correlated to o q 1 EOF 2 with a correlation coefficient of 0.85. Both of these ocean gyre responses were found to have very low correlations of < 0.25 (for any time-lags up to 16-year lags in both directions) indicating that the responses are likely governed by distinct mechanisms.
The maximum lagged coexplained variance of the computed SVD modes peaks at a 2-3-year lag (Fig. 10g) with the upper-isoypcnal ocean PV response found to be correlated to the meridional jet shift mode in Fig. 10. The zerolag and 1-year lag upper-isopycnal ocean PV response was found to be correlated to Fig. 10d indicating the response of the subtropical inertial recirculation zone and gyre acts at a faster time-scale. There is some uncertainty in predicting this dominant time-scale response which is likely due to high levels of mesoscale ocean turbulence. Despite a peak in maximum covariance at 2-year, we notice that both the lagged covariances and correlations of the temporal SVD coefficients persist for up to 16-year, much longer than the expected upper-isopycnal ocean gyre response time-scale (typically less than 5-year, see Hogg et al. 2005;Berloff et al. 2007a;Kravtsov et al. 2007). We attribute this to likely repeated feedbacks of the ocean gyres with the atmosphere. In addition, the sign of the response always remains consistent, with either a weakening or a poleward shift of the WBCE detected for this particular forcing pattern implying a positive feedback signal.
For the middle-isopycnal ocean, the PV anomaly response largely consists of o q 2 EOF 1 activated by the diabatic component of the upper-isopycnal ocean forcing which is given by + o e 1 EOF 1. The PV response is contained within the subtropical gyre as this is where diabatic entrainment dominates (see Figs. 3g,9b,g). The average correlations are 0.64 for o e 1 EOF 1 and 0.99 for o q 2 EOF 1. Note the sign change in vorticity fluxes is consistent with the forcing in the upper-isopycnal layer. The middle-isopycnal ocean PV response shows a negative PV anomaly protruding out from the western boundary extending out into the ocean basin and a similar positive PV anomaly extending out from the eastern boundary. This pattern was found to be correlated to the upper-isopycnal subtropical recirculation zone weakening in the upper-isopycnal layer (Fig. 10d). By checking the streamfunction response in the middle-isopycnal layer, the western and eastern circulation anomalies correspond to cyclonic and anti-cyclonic anomalies, respectively. Such a circulation response in the middle-isopycnal layer hence implies that these circulation anomalies do not arise as a result of diabatic entrainment forcing in the middle-isopycnal layer, but must be due to momentum transfers from the upper-isopycnal layer due to changes in the subtropical gyre circulation. The middle-isopycnal layer response to meridional shifts of the WBCE were not detected when looking at PV anomalies but were detected for transport streamfunction responses. The pattern of the response was similar to Fig. 10e and also occur due to momentum transfers from the upper-isopycnal layer. The time-scale of the middle-isopycnal ocean PV response is less clear than the upper-isopycnal (see Fig. 10). There appears to be a small peak in maximum covariance at the 1-year lag, perhaps associated with a fast barotropic response along the WBCE (e.g. Dewar 2003), but this peak is small and not much greater than the other values of computed maximum covariances (see Fig. 10).
In summary, the east-west positive-negative upperocean total entrainment mode induces a weakening of the subtropical recirculation zone and/or a poleward shift of the WBCE in the upper-isopycnal layer. The middle-isopycnal layer response seems largely to be controlled by momentum transfers that occur due to baroclinic instabilities through the much stronger upper-isopycnal circulation response. Furthermore, the ocean does not appear to respond at any preferred time-scale and instead acts over a wider interannual to interdecadal band time-scale. This is likely due to reduced predictability of the model induced by high levels of mesoscale ocean turbulence (Nonaka et al. 2016). The east-west dipole forcing mode is difficult to compare to most studies of wind-induced ocean gyre variability which mostly focus on monopole forcings (e.g., Dewar 2003;Kravtsov et al. 2007). Both the ocean forcing and response modes are complex in nature and further analysis of the processes involved is required. Since the forcing mode may be split into two separate anomalies, one in the western and the other in the eastern basin, we hypothesise that the two uncorrelated ocean gyre responses are a result of two distinct mechanisms responding independently to the two ocean forcing anomalies. Further dynamical interpretation of this ocean gyre response is outside the scope of this paper and we leave this to Kurashina and Berloff (2023). Instead, we focus our attention on the associated SST response and how this affects the atmospheric circulation.

SST-induced atmosphere diabatic entrainment
As the ocean gyre responses in the previous Subsection induce changes in the shape and location of the WBCE, it is natural to check how the associated SST front changes affect the atmospheric circulation. Since we know the associated sensible and latent heat flux response to the discussed modes of ocean gyre variability (see Fig. 6), we would now like to see what atmosphere diabatic entrainments are associated with them. Figure 11 shows the sensible and latent heat fluxes a F and the consequently generated atmosphere diabatic entrainments a e 1 over the ocean basin. Only lag-zero SVDs are checked since the atmospheric response time-scale is short in comparison to the ocean. The polarities of these SVD modes correspond to the aforementioned, wind-induced changes in the ocean gyre circulation. We have only computed these SVDs using atmosphere diabatic entrainment data over the ocean basin region where the effect of sensible and latent heat fluxes are strongest. Only the positive atmosphere diabatic entrainment situated over the eastern subtropical gyre is unaccounted for by the sensible and latent heat fluxes, but it is present in the associated SST EOF (see Fig. 6f). The damped surface heat fluxes over the subtropical gyre are likely due to the white-noise integrator effect (Hasselmann 1976). The a F modes in Fig. 11a, c are highly correlated to a F EOFs 1 and 2 with correlation coefficients Fig. 11 Leading lag-zero SVD modes of sensible and latent heat flux and diabatic atmosphere entrainment at lag-zero in the lower-isopycnal atmosphere. a, b SVD mode 1 with a F mode presented in a and a e 1 mode presented in b. c, d SVD mode 2 with a F mode presented in c and a e 1 mode presented in d 0.9 and 0.73, respectively. The a e 1 modes in Fig. 11b, d are correlated to a e 1 EOFs 3 and 1 with correlation coefficients of 0.87 and 0.64, respectively. The temporal coefficients of the SVD mode 1 and 2 themselves are also correlated highly with coefficients of 0.85 and 0.85, respectively.
The leading SVD mode pair of a F and a e 1 (Fig. 11a, b) is associated with the weakening of the subtropical recirculation zone and the expansion of the same gyre. The negative sensible and latent heat flux anomaly is associated with a deformation of the recirculation zone and a change in the shape of the WBCE and SST fronts, while the positive heat flux anomaly is associated with a warming of the eastern subtropical gyre as it expands. On the other hand, the second SVD mode pair (Fig. 11c, d) is dominated by strong, positive sensible and latent heat fluxes over the SST front as the WBCE shifts poleward of its average position. This heat flux generates a strong diabatic entrainment anomaly over the same region with a strong diabatic entrainment response of similar meridional scale in the form of an anomalous line source of vorticity (the sign of this depends on whether we are looking at the lower-or middle-isopycnal layer). Although there does appear to be some mixing of this temperature anomaly in the atmospheric mixed layer, the narrow meridional scale of diabatic entrainment is maintained well. This is likely due to the influence of meridional eddy heat fluxes that are anchored to the location of the SST front (Fig. 4g).

Inertial response of westerly jet to diabatic entrainment forcing over ocean basin
Next, we look at the response of the atmosphere to diabatic entrainments. Although we have already correlated EOFs of these fields to each other in Sect. 4.3, we now isolate impacts of SST variability on the atmosphere by spatially subsampling over the ocean basin. This will allow us to better determine the atmospheric response to the line source of vorticity identified in the previous paragraph.
To do this, we computed lag-zero SVDs of atmosphere diabatic entrainment a e 1 and atmosphere PV a q 1 , a q 2 . Similar SVDs with entrainment and dynamic pressure anomalies a p 1 , a p 2 are also computed to see the inertial effects on the atmospheric jet. In the case that influence of SST variability on the atmosphere is small or negligible, we would expect to see the computed SVDs will be similar or identical to the EOFs shown in Fig. 7, i.e. no difference in atmosphere response when subsampling data over the ocean basin. In addition, to ascertain any impact of atmospheric eddies on the response to diabatic entrainments, we also computed a supplementary coupled climate with a moderately reduced atmospheric resolution of 80 km (see Table 3 in Appendix 4). The ocean model parameters are identical to the benchmark modelled climate. The same sets of analysis was performed using data obtained from this supplementary modelled climate.
The leading SVD modes for the 20 km and 80 km lowerisopycnal atmosphere dynamic pressure responses to diabatic entrainments correlated to a e 1 EOF 1 are presented in Fig. 12. Firstly, the atmosphere diabatic entrainment modes (Fig. 12a, c) are highly correlated to equivalent EOFs to the one shown in Fig. 7a with correlation coefficients > 0.9 . We already know that this mode of diabatic entrainment, in the benchmark modelled climate, is well correlated to meridional shifts of the SST front and WBCE (e.g. Figure 11c, d). We found that the structure of the modes for the corresponding PV response in both atmospheric resolutions were similar to those already shown in Fig. 7, so they are omitted. However, we found that the dynamic pressure or inertial response of the 20 km atmosphere to this diabatic entrainment was significantly different in structure to the EOFs shown in Fig. 7. In particular, the lower-isopycnal pressure response shows an asymmetry in the strength of the anomalous response over Fig. 12 Leading lag-zero SVD modes of atmosphere diabatic entrainment and lower-isopycnal atmosphere dynamic pressure anomaly at lag-zero. Top panel row shows the SVD modes for 20 km atmosphere while bottom panel row shows the SVD modes for 80 km atmosphere. Left-hand side panels show a e 1 modes while right-hand side panels show a p 1 modes the western and eastern basins. 9 Over the western basin, there is a weaker cyclonic anomaly, while over the eastern basin, there is a stronger anti-cyclonic anomaly. This is in contrast to the EOFs in Fig. 7c, f which are computed over the entire atmosphere. The EOFs computed in Fig. 7 do not show this asymmetry in response of the westerly jet to diabatic entrainments because ocean-atmosphere interactions are given less importance, in a statistical sense, as the variability of the entire atmosphere is being considered.
Reduction of the atmosphere model resolution, and a consequent damping in eddy-induced behaviours, leads to a reduction in this asymmetric response (Fig. 12d), with the lower-isopycnal pressure response correlated to EOFs of the same structure as those in Fig. 7. On the other hand, the asymmetric inertial response of the 20 km westerly jet (Fig. 12b) appears to be occur as a result of SST variability on the atmospheric westerly jet with the inertial response now correlated to both a p 1 EOF 1 and 2 with correlation coefficients 0.70 and 0.62, respectively. Note that the correlation to a p 1 EOF 1 was undetectable in the EOF analysis in Sect. 4.3. This atmospheric response is akin to those seen in Taguchi et al. (2012) and Famooss Paolini et al. (2022) with the an anti-cycloncic pressure anomaly downstream of the jet. Furthermore, for poleward shifts of the WBCE, the atmospheric westerly jet experiences acceleration at higher latitudes and deceleration at lower latitudes which may have influences on atmospheric blocking (Famooss Paolini et al. 2022). The opposite may be seen for equatorward shifts of the WBCE. Indeed, this ocean feedback is so strong that it likely leads to the reduced zonal structure of the atmospheric westerly jet meridional shift mode seen in Fig. 8 which has a similar spatial structure to Fig. 12b. The structure of the pressure response may also explain the zonally sheared ocean forcing patterns observed in Fig. 10a-c.
We found that a e 1 EOF 3 (Fig. 11d), or changes in strength of the subtropical recirculation zone, did not appear to strongly invoke modes of atmospheric variability despite spatial subsampling over the ocean basin region. We found some evidence of coupling to the westerly jet meridional shift modes in Fig. 8 but correlation of any computed SVD modes to this EOF was low ( ∼ 0.4 ) so these modes were omitted. This does not necessarily mean that this SST anomaly has no measured effect on the atmospheric circulation, but rather, the response of the atmosphere to this particular anomaly is likely more strongly controlled by eddy-mean flow feedbacks. In other words, the response of the atmosphere to this SST anomaly is likely more dependent on the background flow state and is hence nonlinearly related to the strength and sign of the SST anomaly.
Finally, to better understand why improved resolution of atmospheric eddies increases the asymmetry in the strength of the lower-isopycnal atmospheric pressure response, we compare the meridional eddy heat fluxes of the 20 km and 80 km atmospheres that are conditionally time-averaged over years where the WBCE is in its poleward and equatorward states, i.e. greater than or less than ±1 of their measured mode variabilities. These meridional eddy heat fluxes, and their difference between the poleward and equatorward WBCE states are shown in Fig. 13. Reduction of the atmospheric resolution to 80 km severely inhibits the atmosphere's ability to resolve the underlying meridional shifts in the SST front (see Fig. 13c, f), despite the time-averaged state itself being fairly well resolved. The anomalous fluxes between the poleward and equatorward states of the SST front are both smaller in magnitude and their zonal extent. Furthermore, the meridional eddy heat flux profiles for the WBCE in the poleward and equatorward states are almost identical for the 80 km atmosphere. Since the ocean circulation is almost identical with the nominal 5 km resolution maintained in both climates, this difference is attributed, in part, to the lack of the atmospheric eddies' ability to resolve the small meridional scales associated with shifts in the SST front. Indeed, the magnitude of the meridional shift between the two states is around 200 km which is poorly resolved by the 80 km atmosphere. Consequently for the 80 km atmosphere, this leads to a 70 − 80% reduction in strength of the anomalous line sources of heat and vorticity in Fig. 11c, d. Around half of this reduction in strength is due to a corresponding weaker SST anomaly which is likely due to damped positive feedbacks in the 80 km atmosphere climate. The other half of the reduction is unaccounted for and is likely due to inhibited atmospheric eddies. The correlation of the associated temporal SVD coefficients also drops from 0.85 in the 20 km atmosphere to 0.59 in the 80 km atmosphere. This difference in strength of the anomalous line source of vorticity likely explains the two different inertial responses of the atmospheric jet with the increased zonal asymmetry acting as a proxy for the atmosphere's sensitivity to underlying shifts of the SST front. Indeed, if the SST front was fixed in its meridional position, which is similar to the response of the 80 km atmosphere, we would expect to see pressure anomalies of similar strength over the western and eastern ocean basin. In addition, the reduction in correlation between the surface heat flux anomaly and associated westerly jet response is indicative of a stronger dependence on the background flow state (Peng et al. 1997;Peng and Whitaker 1999;Seo et al. 2017).

Summary of ocean gyre influences on the climate variability
Our analysis has revealed the following positive feedback mechanism for enhanced low-frequency climate variability: 1. East-west positive-negative PV anomalies arise in the middle-isopycnal atmosphere over the ocean basin as part of a wavenumber-6 Rossby wave disturbance. 2. This propagates into the lower-isopycnal atmosphere via baroclinic instabilities, and surface stresses generate positive and negative vorticity fluxes over the western and eastern ocean basins, respectively. The growth of baroclinic instabilities is significantly strengthened over the ocean basin due to the presence of the SST front. 3. The ocean gyres respond through a weakening of the subtropical recirculation zone and a poleward meridional shift of the WBCE. Full dynamical interpretation of this ocean gyre response may be seen in Kurashina and Berloff (2023), but preliminary analysis indicates that these two responses are governed by distinct mechanisms. 4. SST anomalies are generated as a consequence of shifts in the inter-gyre boundary and deformations of the WBCE. This then leads to sensible and latent heat flux anomalies that generate strong diabatic entrainment anomalies in the atmosphere. The sensitivity of the atmosphere to underlying SST front shifts is highly dependent on resolution. 5. The diabatic entrainments trigger baroclinic instabilities that grow downstream of the westerly jet and eventually leads to a strong barotropic circulation over the eastern basin. This atmospheric response is shown to be modified by eddies and dependent on the atmosphere's ability to resolve meridional scales associated with the SST front shift. 6. This atmosphere response then reinvigorates the PV variability and generates a positive feedback onto the gyres.

Summary and discussion
Motivated by the increasing body of evidence for coupling of atmospheric jets to WBCEs (Frankignoul et al. 2011;Taguchi et al. 2012;Révelard et al. 2018) and the potential of atmospheric mesoscales to strengthen tropospheric forcing (Czaja et al. 2019), we investigated the low-frequency, midlatitude climate variability using an idealised, eddy-resolving coupled model. Although the model we have used in this study is highly idealised, e.g., lack of moisture, simple model geometry, lack of overturning circulation and vertical motions, proper capturing of complex eddy dynamics in the ocean and atmosphere have revealed relevant coupled interactions that enhance the low-frequency variability. We believe these models play an important role in understanding the midlatitude climate and that this study opens up more avenues of future work. Part I of this two-paper study focused on the effects of ocean-atmosphere coupling on the modelled climate and established the anatomy of the low-frequency climate variability. EOF analysis of the 2-year low-pass filtered, atmosphere PV variability showed two wavenumber-6 Rossby wave modes shifted out of phase in the longitudinal direction. These EOFs consisted of patterns with alternatingsigned circulation anomalies, similar to those found in Hogg et al. (2006). Despite these EOFs appearing at first to form a travelling-wave pair, the disparity in the total explained variance in temporal behaviour meant that the EOFs should be treated separately. The difference in explained variance is due to stronger coupling of one of the modes with the ocean gyre circulation. Ocean gyre coupling with the atmosphere also likely leads to the reduced zonal symmetry of the atmospheric jet meridional shift pattern in the benchmark modelled climate. Previously, only time-scale modifications of atmospheric modes have been found using this class of models (e.g., Hogg et al. 2006) so this finding is, to our knowledge, new and important.
EOF analysis of the filtered, upper-ocean PV variability data showed the two leading EOFs consisted of changes in strength of the subtropical inertial recirculation zone and meridional displacements of the WBCE. These EOFs were both found to be strengthened significantly by coupling over the entire low-frequency spectrum. Such modes of variability have been seen frequently in QG models (Hogg et al. 2005;Berloff et al. 2007b) as they appear from the intrinsic low-frequency variability of the ocean gyres.
The ocean was found to respond in two possible ways to the east-west positive-negative dipole forcing over the basin: (1) weakening of the subtropical recirculation zone and shrinking of the respective gyre at a relatively fast 0 − 1 -year response time-scale; (2) a poleward meridional shift of the WBCE at a relatively slow 2-3-year response time-scale. However, the estimation of this response is uncertain with the associated correlation coefficients remaining at fairly similar levels over the entire range of computed time-lags. This is likely due to high levels of mesoscale turbulence that reduce ocean gyre predictability (Nonaka et al. 2016) and repeated positive feedbacks present in the modelled climate. The sensible and latent heat flux response associated with the meridional shift of the WBCE generated atmosphere diabatic forcings that reinforced the initial dipole forcing, and thus, reinvigorate the ocean gyre modes. Again, this is consistent with our findings of increased, broadband ocean gyre variability through coupling. The nature of this positive feedback between the ocean and atmosphere bears some resemblance to the interannual coupled interaction found by Goodman and Marshall (1999), Kravtsov et al. (2006) and Farneti (2007). This coupled interaction involves the westward propagation of coupled Rossby waves and enhanced interannual variability due to positive feedbacks. The timescale of this interaction is found to be controlled by the basin crossing time of the ocean baroclinic Rossby wave. The coupled interaction in our model appears similar, with an initial ocean gyre response time-scale of 2-3-year which may be due to wave propagation, but positive feedbacks persist for up to much longer interdecadal time-scales. In addition, the atmospheric mode of the coupled interaction in our model appears to be a standing wave Rossby wave mode, rather than a phase-locked mode seen by Goodman and Marshall (1999), Kravtsov et al. (2006) and Farneti (2007). In fact, the coupled interaction observed in this paper more closely resembles the positive feedback found in more recent reanalysis data Taguchi et al. (2012) between the Subarctic Frontal Zone (SAFZ) and Pacific North American (PNA) pattern. This PNA pattern response to shifts in the SAFZ has also been verified in atmospheric GCMs (Tatebe et al. 2017;Okajima et al. 2018), albeit at a weaker magnitude, which is most likely due to the relatively low resolution of the models used. Similar responses have also been observed for Gulf Stream shifts in the North Atlantic using a high-resolution atmospheric GCM by Famooss Paolini et al. (2022).
Their findings show a circulation response that is closer to observations.
The asymmetry in the atmosphere response, where the westerly jet is displaced less over the western basin and more over the eastern basin, is strongest in the lower-isopycnal layer due to the reduced momentum of the lower-isopycnal jet. The final pressure anomaly in the lower-isopycnal atmosphere is a barotropic circulation anomaly over the eastern basin. This asymmetry in the response is so much so that the lower-isopycnal pressure response is now also correlated to the modified westerly jet meridional shift mode. This mode was only detected after computing SVDs using data subsampled over the ocean basin region only, rather than the entire atmosphere. The structure of this westerly jet meridional shift EOF is different in structure to modes seen in Kravtsov et al. (2007) because it is less zonally symmetric.
The response of the atmosphere to meridional shifts in the SST front was also found to be dependent on the atmospheric resolution. For resolutions where the meridional scale of the shift is well-resolved, shifts in the low-level baroclinicity translate to a weaker inertial deflection of the jet over the SST front and a stronger inertial deflection of the jet downstream. The asymmetric inertial deflection of the jet is likely due to downstream organisation of the eddyinduced, line source of vorticity generated by the SST front shift. When the meridional scale of the SST front shift is poorly resolved, it is now unable to detect the changes in the low-level baroclinicity as well. So although an entrainment mode still exists, this mode is no longer controlled by meridional shifts of the WBCE, and the anomalous line source of vorticity is significantly weakened. Note that although the diabatic entrainment is still strengthened by the presence of the SST front which restores the low-level baroclinicity, the control of the variability is now handed over to the atmosphere, likely through a stronger dependence of the atmospheric response to the background flow state. This lack of atmospheric sensitivity to meridional shifts in the SST front at 80 km resolution may be a contributing factor to why previous studies using this model struggled to capture many of these feedbacks of SST variability on the atmosphere (Hogg et al. 2006). Inability of atmospheric eddies to resolve shifts in the SST front may be a contributing reason to weaker-than-observed PNA responses observed by Tatebe et al. (2017), Okajima et al. (2018).
The results in this paper are also consistent to findings by Piazza et al. (2016), O'Reilly et al. (2017 where more realistic, sharper SST fronts lead to greater poleward excursions of the atmospheric jet. Greater poleward excursions of the jet have been explained by increased eddy heat fluxes due to a sharper SST front (O'Reilly et al. 2017). This then leads to a greater anomalous line source of vorticity which is advected downstream by the atmospheric jet.
Although Part I of this study has revealed the anatomy of delayed, wind-induced ocean gyre variability, we are yet to show the ocean processes driving the ocean gyre modes. We have proposed that certain patterns of ocean forcing activate modes of ocean gyre variability but exactly how they do this is still unclear and further analysis of the underlying dynamics is required. To do this, Part II of this study will require a series of ocean-only simulations run under fixed windstresses and thermal forcings to determine whether we are able to reproduce the effects of ocean entrainment observed in the coupled model. An analysis of ocean PV budgets will also be required to quantify the sources and sinks of PV within each of the gyres (see Kurashina and Berloff 2023, for details).
where F ↑ is instantaneous diabatic heating and F ↑ is the time-averaged diabatic heating obtained from a fully-coupled simulation.

Appendix 3: Modelled wind-stresses for benchmark modelled climate
The time-averaged wind-stress is largely zonal with strong eastward wind-stresses focused in a fairly narrow band that sits predominantly over the subpolar gyre (Fig. 14). The meridional wind-stress has time-averaged stresses of just over a quarter in magnitude of the zonal wind-stress with an average forcing in the poleward direction. This is associated with the addition of land orography west of the ocean basin which creates a mean flow that has a poleward wind-stress component. This figure has been retained for the purposes of Kurashina and Berloff (2023) where the time-averaged wind-stresses are used to define the forcing for an ocean-only configuration of Q-GCM.

Appendix 4: Reduced atmospheric resolution parameters
See Table 3.
Author contributions RK performed the calculations, prepared the figures and wrote the manuscript. All authors reviewed the manuscript.
Funding R.K. was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Mathematics of Planet Earth. P.B. was supported by the NERC grant NE/T002220/1 and by the Moscow Centre for Fundamental and Applied Mathematics (supported by the Agreement 075-15-2019-1624 with the Ministry of Education and Science of the Russian Federation).

Availability of data and materials
The model used in this study is available to download here.

Conflict of interest
The authors declare that they have no competing interests.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will