Low-frequency variability enhancement of the midlatitude climate in an eddy-resolving coupled ocean–atmosphere model—part II: ocean mechanisms

This paper investigates the spatial inhomogeneity of the time-averaged, quasigeostrophic, double-gyre circulation response to fixed, realistic, large-scale modes of wind-stress forcing. While the companion paper of this study focused on understanding the anatomy of low-frequency, midlatitude climate variability in an idealised, eddy-resolving coupled model, this paper looked at understanding the nature of the wind-induced ocean gyre response using an ocean-only configuration of the same model. Our analysis revealed two, time-averaged responses to an east–west dipole, wind-stress curl anomaly in the ocean basin. Firstly, wind-stress anomalies in the western ocean basin led to changes in relative strength of the inertial recirculation zones and jet-axis tilt. This is consistent with an advection-dominated, nonlinear adjustment of the ocean gyres to anomalous forcing. Secondly, wind-stress curl anomalies in the eastern ocean basin was found to induce a largely independent response involving meridional shifts of the western boundary current extension (WBCE). The effects of time-averaged advection in this region are weak and the discovery of westward-propagating Rossby waves along the WBCE revealed the response is more akin to a baroclinic Rossby wave adjustment.


Background
The formation of the linear, Sverdrup gyres with Munk boundary layer approximation under fixed wind-stresses is simple and well known. Given some prescribed doublegyre forcing, the ocean gyres create two large pools of opposite-signed potential vorticity (PV) with powerful western boundary currents (WBCs) that occur as a consequence of western intensification. The vorticity sources and sinks for the gyres are such that wind-induced vorticity flux, in the form of Ekman pumping and/or diabatic entrainments, is balanced by viscous boundary fluxes primarily along the western boundary. The shape of the gyres is determined by the structure of the wind-forcing, with the inter-gyre boundary largely aligning with the zero vorticity flux line.
The addition of nonlinear dynamics (e.g. Veronis 1966a, b;Böning 1986;Haidvogel and Holland 1978;Holland 1978;Holland and Rhines 1980) leads to adjustments of the western boundary layer (WBL) structure and the formation of the western boundary current extension (WBCE) and inertial recirculation zones. PV advection from the WBCs is now advected into the inertial recirculation zones, which occurs due to insufficient relative vorticity loss at the western boundary. This is caused by effects of complex nonlinear dynamics in the WBL (Cessi et al. 1987;Lozier and Riser 1989;Nakano et al. 2008;Kurashina et al. 2021). Within the inertial recirculation zones, processes due to mesoscale dynamics such as eddy backscatter and inter-gyre PV exchanges occur. The competition of these two processes, with eddy backscatter fortifying the WBCE (Berloff 2005(Berloff , 2016Shevchenko and Berloff 2016) and inter-gyre PV exchange weakening it (Berloff et al. 2007), have been shown to be crucial in driving decadal ocean gyre variability. These processes are usually triggered by anomalous wind-forcing, which leads to nonlinear adjustments of the ocean gyres (Dewar 1 3 2003). This involves a fast, barotropic eddy response along the WBCE as the inter-gyre boundary reacts to anomalous forcing. This is followed by a slower, baroclinic response involving deep PV redistribution in the gyres. The windforcing anomalies may arise in both ocean models driven by fixed winds, e.g. through misalignments of the zero vorticity flux line and inter-gyre boundary (Berloff et al. 2007), or through modes of large-scale atmosphere variability such as meridional shifts of atmospheric jets (Czaja and Marshall 2001). Crucially, it is the strong advective effects within the inertial recirculation zones that is shown to be important in driving this decadal ocean gyre variability.
There are many more specific examples of how the inertial recirculation zones respond to anomalous windforcing. For example, asymmetric wind-forcing has been shown to induce deviations from straight, west-to-east WBCEs due to differences in relative strength of inertial recirculations (Moro 1988;Hogg et al. 2005). This is because the vorticity flux imbalance creates a stronger recirculation zone, which leads the WBCE to 'wind-up' around itself (Moro 1988). Other effects of asymmetric wind-forcing include WBC overshooting (Harrison and Stalos 1982) and decreases in WBCE penetration (Rhines and Schopp 1991). Ekman pumping anomalies near the jet separation region have also been shown to destabilise the WBCE and weaken the ocean gyre circulation (Hogg et al. 2009).
Another proposed mechanism of decadal ocean gyre variability is through baroclinic Rossby wave adjustments (Sasaki and Schneider 2011;Sasaki et al. 2013). This mechanism occurs due to the formation of a basinscale, baroclinic Rossby wave in the eastern ocean basin. This Rossby wave then propagates westwards along the WBCE, narrowing in meridional scale as it does so in order to conserve PV along its trajectory. For example, an anti-cyclonic wind-stress curl in the eastern North Pacific has been linked to a poleward meridional shift in the Kuroshio extension (Sasaki et al. 2013). A strengthening of the WBCE was found to occur very shortly afterwards with an eastward propagating signal, indicating the importance of advective effects by the WBCs. Such an acceleration of the WBCE has been suggested to be triggered by the arrival of the Rossby wave signal at the western boundary, which leads to local restructuring of the WBL (Taguchi et al. 2005(Taguchi et al. , 2007. In these studies, it is implied that advective effects are secondary to the Rossby wave propagation and that effects of advection may only be triggered by the arrival of the wave at the western boundary. In order to motivate and give context to our subsequent analysis, we will now also give a brief recap of the main results in Part I of this study (see Kurashina and Berloff 2023, for full details).

Summary of wind-induced ocean PV variability from part I
Part I of this study (Kurashina and Berloff 2022), hereby referred to as KB22, focused on understanding the anatomy of midlatitude climate variability and coupling of the wind-driven ocean gyres to the atmospheric westerly jet. The 2-year low-pass filtered upper-isopycnal ocean PV variability was found to be dominated by two statistically independent modes or empirical orthogonal functions (EOFs) (see von Storch and Zwiers 1999). These modes were not found to form travelling wave-like behaviours and showed a lack of correlation at any time-lag. The first mode involved changes in strength of the subtropical recirculation zone, as well as changes in size of the respective gyre, while the second mode involved meridional shifts of the WBCE. The low-pass filtered atmosphere PV variability was found to be dominated by two wavenumber-6 modes with one of these modes coupling strongly with the ocean gyre dynamics. These modes of the atmosphere were found to transfer momentum into the lower atmosphere through relatively fast events of baroclinic instabilities. Momentum is then transferred into the ocean through wind-stress anomalies that arise through frictional effects. Modes of covariability between wind-induced ocean forcing and the ocean gyre response were found through lagged singular value decompositions (SVDs) (Bretherton et al. 1992). The ocean and atmosphere modes associated with one phase, or polarity, of this wind-induced ocean gyre variability are presented in Fig. 1 and involve the aforementioned modes of low-pass filtered ocean and atmosphere PV variability. Figure 1a, b show the windstress modes generated by the atmosphere PV variability, with cyclonic wind-stresses over the western ocean basin and anti-cyclonic wind-stresses over the eastern ocean basin. This then results in east-west dipole wind-stress curl anomalies in the upper-isopycnal ocean (Fig. 1c), with a corresponding weaker and opposite-signed diabatic entrainment forcing in the middle-isopycnal ocean (omitted).
Two delayed, responses were found in the upper-isopycnal ocean gyres to this forcing. The first was a weakening of the negative PV-signed subtropical recirculation zone (Fig. 1d); the second response of the ocean gyres was a poleward meridional shift of the WBCE (Fig. 1e). An equal and opposite response was found for a forcing anomaly of the opposite sign.
The middle-isopycnal ocean PV response was found to be contained within the subtropical gyre where diabatic entrainment dominates. This response involved the PV redistribution of vorticity flux anomalies in the subtropical gyre with a negative-signed PV anomaly protruding from the western boundary and a positive-signed PV anomaly protruding from the eastern basin (Fig. 1f). This mode was found to be correlated to the upper-isopycnal PV response in Fig. 1d. A barotropic response in the middle-isopycnal layer for Fig. 1e was found when checking the streamfunction response to the same mode of ocean forcing (omitted).
Although a robust statistical link was established between these modes of ocean forcing and response variables, full interpretation of the mechanisms and physical processes responsible in driving these changes in the ocean has yet to be looked at. This is because both the nature of the ocean forcing and response was complex. Firstly, the ocean forcing mode contains two, oppositesigned, vorticity flux anomalies situated in the western and eastern ocean basins. Secondly, there are two uncorrelated (at any time-lag) upper-isopycnal ocean gyre responses to this anomalous forcing. Clearly, these results are not easily explained by any single mechanism proposed by studies involving monopolar vorticity fluxes (e.g. Dewar 2003;Sasaki et al. 2013), and are thus possibly governed by more than one mechanism. The simplest explanation for these observations is that two distinct mechanisms are involved in governing the delayed, ocean gyre responses. One for the western basin vorticity flux anomaly and one for the eastern basin vorticity flux anomaly. Our main aim in this study is to investigate the validity of this hypothesis through the modelling of multiple, fixed-wind, doublegyre circulations.
In Sect. 2 the models and methods used in this paper are discussed. Section 3 shows the reference solutions to the modelled ocean circulation; Sect. 4 shows the time-averaged responses of the double-gyre circulation under different Modes obtained from analysis in KB22. a Zonal ocean surface-stress anomaly. b Meridional ocean surface-stress anomaly. c Ocean Ekman pumping. d Upper-isopycnal ocean PV response mode 1. e Upperisopycnal ocean PV response mode 2. f Middle-isopycnal ocean PV response 1 3 fixed wind-stress forcings. A discussion and summary of the main results are given in Sect. 5.

Notation
The notation we used in Part I of this study included leftsuperscripts of variables such as 'o' or 'a' for ocean and atmosphere, as well as right-subscript, 'm' or integer i, to indicate the mixed layer or i-th isopycnal layer, respectively. Layers were counted away from the ocean-atmosphere interface. In this paper, we will be restricting our focus on ocean variables, only so we will drop the left superscript and assume all variables are oceanic.

The quasi-geostrophic coupled model (Q-GCM) [ocean-only configuration]
The model we will be using for Part II of this study is the ocean-only configuration of Q-GCM (Hogg et al. 2014) run under the same ocean model parameters as KB22. Details of the coupled configuration of Q-GCM is covered in KB22 but we will give a brief summary of the ocean-only configuration below.
The ocean-only configuration of Q-GCM consists of an idealised, eddy-resolving, quasigeostrophic model designed to mimic the midlatitude wind-driven ocean circulation. It consists of a double-gyre box ocean, which is driven by prescribed, fixed mechanical and thermal forcing. The box ocean model also consists of an active mixed layer, which is embedded within the upper-isopycnal layer. The mechanical forcing is given by wind-stresses, which are then converted to Ekman pumping and mixed layer velocities. The thermal forcing consists of some prescribed surface heat flux which contains sensible, latent, radiative and solar heat fluxes. This forces the sea surface temperature (SST) evolution in the ocean mixed layer as well as indirectly influencing entrainment velocities. Since model forcings are fixed and the atmospheric model is switched off, any influences on the ocean circulation through ocean-atmosphere coupling or atmospheric variability are damped.

Model geometry
The box ocean model consists of 3 isopycnal layers with depths H 1 = 350 m, H 2 = 750 m and H 3 = 2900 m . The horizontal grid spacing for our benchmark ocean circulation is chosen to be 5 km for a basin with dimensions 4800 km × 4800 km , which gives a discretised grid of size 961 × 961 points representing cell vertices in the longitudinal and latitudinal directions, respectively.

The quasigeostrophic (QG) layers
The ocean circulation evolution is governed by the quasigeostrophic equations in each layer (e.g. Pedlosky 1987). These are written in terms of PV anomalies and dynamic pressure anomalies denoted in vector form as q = (q 1 , q 2 , q 3 ) T and p = (p 1 , p 2 , p 3 ) T , respectively. PV anomalies are defined as where f 0 = 9.37 × 10 −5 s −1 is the Coriolis parameter. Full PV is obtained by adding the term y ⋅ 1 onto (1) where = 1.75 × 10 −11 m −1 s −1 is the planetary vorticity gradient. We use (x, y) as our zonal and meridional coordinates where x is always measured increasing from the western boundary and y is measured increasing away from the central latitude y 0 . We also have 1 = (1, 1, 1) T and A is a matrix which defines interactions between layers through dynamic pressure anomalies where H i is the unperturbed thickness of the i-th layer and g ′ i are the reduced gravities. In the ocean, g � 1 = 0.0222 m s −2 , g � 2 = 0.0169 m s −2 correspond to Rossby deformation radii of 40.0 km and 20.6 km , respectively.
The QG equations may now be succinctly written as where A 2 = 50 m 2 s −2 and A 4 = 2.0 × 10 9 m 4 s −2 are the Laplacian and biharmonic diffusion coefficients, respectively. J(f, g) = f x g y − f y g x is the Jacobian operator. The forcing in the isopycnal layers is defined by the 3 × 4 matrix B and the 4 × 1 vector e. The forcing in (3) is defined by B which is a 3 × 4 matrix given by Next, we have the forcing vector given by The forcing terms consist of Ekman pumping in the first term, which are calculated using wind-stresses and diabatic entrainments in the second term. Upper-ocean total entrainment is defined as the combined adiabatic Ekman pumping and diabatic entrainment forcings Ekman pumping is computed as the curl of wind-stress In the ocean-only configuration of Q-GCM, these windstresses are fixed by some prescribed ( x , y ) of our choosing. For our benchmark solution, we set ( x , y ) = x cpl , y cpl (see Fig. 2a, b) which are the time-averaged wind-stresses obtained from the benchmark coupled climate of KB22. Note that due to the basin being shifted equatorward from the central latitude, there is an asymmetry in the wind-stress forcing which causes the maximum wind-stress to be located ∼ 500 km poleward of the central latitude (see Fig. 2). The diabatic entrainment consists of a scaled product of ocean Ekman pumping and the temperature difference between the ocean mixed layer and the upper-isopycnal layer. We define Δ k T = T k − T k+1 for k = m, 1, 2 as the temperature change across the mixed and isopycnal layers. This is fixed for k = 1, 2 as the temperatures of the isopycnal layers are fixed and variable for k = m due to the time evolution of SSTs. This makes e 1 nonlinear in since SST implicitly depends on wind-stress.
Partial-slip boundary conditions and mass conservation constraints are applied on the lateral boundaries (Haidvogel et al. 1992;McWilliams 1977). The choice of partial-slip boundary condition parameter is identical to KB22.

The ocean surface mixed layer
The temperature evolution in the ocean surface mixed layer, denoted by T m , is given by where T m denotes SST and (u m , v m ) are the mixed layer velocities. The Laplacian and biharmonic temperature diffusion coefficients are denoted by K 2 = 50 m 2 s −1 and K 4 = 2.0 × 10 9 m 4 s −1 , respectively. We will present mixed layer temperatures as anomalies from a constant temperature which is obtained through radiative balance (see Hogg et al. 2014, for details). The thermal forcing term is encapsulated by the term F which contains sensible, latent, radiative and solar heat fluxes. In the ocean-only configuration this is fixed by some prescribed, fixed thermal forcing. For our benchmark solution and subsequent experiments under different mechanical forcings in Sect. 4, we set F = F cpl (see Fig. 2d), which is the time-averaged diabatic heating obtained from the reference coupled model solutions (Kurashina and Berloff 2023). These two terms F and fully define the forcings in this model. The ocean mixed layer velocities (u m , v m ) are computed through the Ekman balance integrated over the mixed layer depth H m (Hogg et al. 2014).

Fig. 2 Prescribed ocean forcings used in the benchmark solution for
the ocean-only configuration of Q-GCM. Note that the Ekman pumping forcing in c is not directly passed into the model, but rather, it is the wind-stress fields (a, b) that are used. These prescribed forcings are obtained from the time-averaged benchmark solutions obtained from the coupled configuration of Q-GCM in KB22. a Zonal windstress. b Meridional wind-stress. c Ekman pumping. d Diabatic heating

Switching off convection
Given an unstable temperature stratification of the ocean, i.e. when T m falls below T 1 , a convective event is modelled by adding a correction term e 1 . However, we found that this led to unrealistically large levels of modelled convection of nearly two orders of magnitude greater than expected around the jet separation region, which was not observed in KB22 (see Fig. 12 in Appendix). In order to manage this undesirable effect, we switched off convection in our modelled double-gyre circulations by setting the entrainment correction term to zero (see Appendix for full details).

Potential vorticity budgets
In KB22, we related wind-induced vorticity flux anomalies with changes in shape and strength of the WBCE. In order to better determine the physical processes involved, a potential vorticity budget is required. These budgets will be timeaveraged budgets over a given time period for the simulated double-gyre circulation in statistical equilibrium. This will allow us to fully quantify any potential vorticity fluxes which arise due to different physical processes.

Inter-gyre boundary
In order to compute potential vorticity budgets as explained in the subsequent paragraphs, an inter-gyre boundary Γ which partitions the ocean basin is required. This is computed in the same manner as in Kurashina et al. (2021) where the time-average contour emanating from the western boundary is followed eastwards.

Source and sink terms
As we are interested in wind-induced phenomena, our focus will lie in the upper-and middle-isopycnal layers which are driven directly by the atmosphere. The governing PV equation in the upper-isopycnal layer is given by A similar equation may be written for the middle-isopycnal layer. Note that the advection of full PV q + y is considered here.
We will now break down the source and sink terms in (9). Similar terms may be obtained for the middle-isopycnal layer. The two source terms in our equation are given by Ekman pumping w ek and entrainment e 1 = − Δ m T 2Δ 1 T w ek . For a given inter-gyre boundary Γ , the gyre-integrated PV (9) fluxes, denoted G, are computed by integrating the source terms within the region encapsulated by the ocean boundaries and the inter-gyre boundary: These two terms are combined into a single wind-induced PV flux source term to give In the middle isopycnal layer, there is no adiabatic Ekman pumping term present and the sign of the diabatic entrainment flux is flipped (5).
Moving onto the sink terms, the viscous boundary fluxes are computed by first integrating the last two terms in (9). Integrating the Laplacian viscosity term over the gyre gives: Similarly, integrating the biharmonic viscosity term gives: where first and second terms in (14) and (17) represent the viscous boundary and inter-gyre viscous fluxes for the Laplacian and biharmonic viscosity terms, respectively. The inter-gyre viscous fluxes were found to be negligible and not considered any further. The viscous boundary fluxes are combined into a single term The inter-gyre flux of full PV q + y is computed by integrating the Jacobian term in (9) (10) Finally, the tendency term is computed by either integrating q 1 t over the gyre or equivalently, and more simply, by summing the above PV fluxes This completes the sources and sink terms in the PV budget. Multiple PV budgets will be computed in this study which each correspond to different Ekman pumping forcing regimes due to the atmospheric variability.

Defining ocean forcings for double-gyre experiments
Although our benchmark ocean circulation is computed using the time-averaged wind-stresses ( x cpl , y cpl ) obtained from the benchmark coupled climate in KB22, we will also compute other double-gyre circulations by adding on some wind-stress anomaly ( �x , �y ) onto ( x cpl , y cpl ) . Since the primary aim of this study is to understand the spatial inhomogeneity of the ocean gyre response to wind-stress anomalies found in KB22, we will base our forcings on the aforementioned anomalies shown in Fig. 1a, b. We will check the effect of both positive and negative polarities of these windstress anomalies for completeness and scale them according to the standard deviations of their respective principal components ( x , y ) (von Storch and Zwiers 1999). This is computed using data obtained from KB22. Further supplementary double-gyre circulations with modified wind-stresses derived from Fig. 1a, b will also be computed later on in the study (see Sect. 4.5).
Although the diabatic heating term F is able to influence the diabatic entrainment e 1 , we will assume this effect is relatively small and leave it unchanged in order to isolate the effects of wind-stress anomalies on the gyres. This removes adding unnecessary variables to the experiments that may over-complicate the dynamical picture.

Benchmark double-gyre circulation
We now present the benchmark double-gyre circulation using the model described in Sect. 2.2 with a simulation period of 120-years, including a 20-year spinup period. The chosen fixed-in-time mechanical wind-stress and diabatic thermal forcings are given in Fig. 2. The instantaneous and time-averaged PV anomalies are shown in Fig. 3, as well as their constituent relative vorticity and buoyancy components. Figure 4 shows the instantaneous and time-averaged transport streamfunction fields. In the upper-isopycnal layer, the ocean gyres consist of large pools of positive and negative PV anomaly in the subpolar and subtropical gyres with inter-gyre exchanges of PV anomalies controlled by mesoscale eddies which are shed off by the powerful WBCE. The middle-isopycnal layer is driven by diabatic entrainment of opposite sign to the upper-isopycnal layer. The lowerisopycnal layer is not directly forced by Ekman pumping or diabatic entrainments and is only driven by nonlinear dynamics due to eddy effects (Holland and Rhines 1980). Figure 3g-l show the time-averaged relative vorticity and buoyancy components of PV anomaly in each layer. Relative vorticity appears weak in the vast majority of the ocean basin except near the western and southern boundaries where they lead to the generation of sharp relative vorticity gradients. This induces large viscous boundary fluxes of PV which drains the gyres of enstrophy. There are also relative vorticity contours situated along the WBCE which is advected from the WBCs. This is because although the WBLs drain PV from the gyres, the effects of nonlinear dynamics in this region inhibits its ability to do so (Cessi et al. 1987;Lozier and Riser 1989;Kurashina et al. 2021). On the other hand, buoyancy appears dominant in the remainder of the gyres and along the WBCE and recirculation zones. Large pools of time-averaged buoyancies in these regions allows for the development of baroclinic instabilities which propagate eddy momentum into the lower layers through form stresses (McWilliams 2008).
The instantaneous and time-averaged SSTs are also presented in Fig. 5a, b. SSTs have a similar structure to Fig. 3a, consisting of a warmer subtropical gyre which is separated from the cooler subpolar gyre by the SST front. The main effect of this SST front is to shift the vorticity fluxes from dominating in the subpolar gyre, as is the case for Ekman pumping (Fig. 2c), to dominating in the subtropical gyre (Fig. 5c). This creates a roughly equal total vorticity flux for each gyre in the upper-isopycnal layer but creates a significant vorticity flux imbalance in the middle isopycnal layer which only receives diabatic forcing.
Qualitatively, this benchmark ocean circulation is similar to the benchmark ocean circulation in KB22. This is as expected as the parameter choices for the ocean model are identical for both. One difference appears to be that the WBCE and recirculation zones appear stronger in the oceanonly configuration of the model, with the time-averaged mass transport increasing by around 10% in each layer. This is expected since time-varying ocean forcing, which is the case in the coupled configuration of this model from KB22, has the effect of destabilising the WBCE, which causes it to break apart more often. For example, highly variable Ekman pumping anomalies near the WBC separation region have been shown to weaken the global ocean gyre circulation (Hogg et al. 2009). We also observe that there are sharper PV and SST gradients across the WBCE and SST fronts, respectively. This is also due to reduced inter-gyre exchanges of PV and heat across the more stable WBCE.
Time-averaged PV budgets of the reference double-gyre circulation (Table 1) shows wind-induced vorticity flux is balanced primarily by inter-gyre PV fluxes and secondarily by viscous boundary fluxes. Note the vorticity flux input sign into the gyres switches due to the sign change in diabatic entrainment. The vorticity inputs into the middle-isopycnal layer are around a quarter of that of the upper-isopycnal layer. This is since adiabatic Ekman pumping, which consists of around half of total vorticity input, is no longer present, as well as the middle-isopycnal layer being around twice as thick as the upper-isopycnal layer. The standard deviations of the gyre-integrated inter-gyre PV flux were found to be 1-2 orders of magnitude larger than the other sources and sinks of vorticity, meaning that inter-gyre PV exchanges control the variability of the ocean gyres (see Table 1). In this configuration of the model, the wind-induced vorticity flux variability is smaller than what would be expected in KB22 since the wind-stresses are fixed. Hence, the only sources of variability in the wind-induced vorticity flux is through shifts in the inter-gyre boundary and changes in SST distribution across the basin (5).

Experiment design
We now run two double-gyre circulation experiments with the different mechanical forcings given in Table 2. The model parameters are identical to those used in the benchmark simulated ocean circulation with the only variable being the mechanical forcing, or wind-stresses, while keeping the thermal forcing fixed. Simulation periods of 120-years, with a 20-year spinup period, are chosen to reduce uncertainty due to strong decadal variability of the gyres. In particular, we would like to show that eastern and western basin vorticity fluxes generate two distinct responses in the upper-isopycnal ocean. After checking the responses of the double-gyre circulation under these different wind-stress forcings, we will compute time-averaged, gyre-integrated PV budgets of the upper-and middle-isopycnal layers. The lower-isopycnal layer budgets are neglected as they are not directly forced by winds.
From these two experiments run under different mechanical forcings, we are interested in: (1) Are we able to reproduce similar ocean gyre responses, in a time-averaged sense, to those found in KB22? (2) What are the associated changes in the PV budgets relative to the benchmark circulation?; (3) What are the likely physical processes involved?

Time-averaged ocean forcings under prescribed
wind-stresses Figure 6 shows the time-averaged modelled Ekman pumping, diabatic and upper-ocean total entrainments. The modelled Ekman pumping forcing is similar in structure to that of Fig. 1c which is expected since the corresponding windstress and Ekman pumping modes were found to be highly correlated in KB22. Note that our modelled Ekman pumping is 15% weaker than what is observed in the coupled model, which is likely due to the lack of positive feedbacks. These east-west dipole Ekman pumping anomalies that sit largely over the subpolar gyre are then pushed equatorward and into a narrower latitudinal band by the SST distribution in the form of diabatic entrainments. Although the structure of this forcing is modified and strengthened by the SST field, the east-west dipole forcing pattern is still maintained. However, we note the diabatic forcing is a considerable 60% stronger than what is measured in the coupled model. This is likely due to the highly idealised formulation of diabatic heating in the ocean (see Sect. 2.2.3). Overall, the modelled ocean forcings show that the simple addition of wind-stress anomalies, without modification of diabatic heating, is able to successfully reproduce east-west dipole vorticity flux patterns in the ocean basin. Although the modelled ocean forcings are not entirely perfect, e.g.   The Control forcing is cpl in Fig. 2a, b. The wind-stress anomalies are given in Fig. 1a, b and are scaled by 1 standard deviation of their measured variability in KB22 Experiment no.

Wind-stress
Control (Reference solution) = cpl 1 = cpl + eof 2 = cpl − eof overly strong diabatic entrainment, our aim is not in perfectly reproducing the exact ocean forcings, but in understanding the response of the ocean gyres to vorticity flux anomalies over the western and eastern basins. Hence, these modelled ocean forcings will suffice for the needs of our study.

Time-averaged circulation response to wind-stress anomalies
The time-averaged response of the upper-isopycnal ocean circulation is shown in Fig. 7. Note that the ocean gyre circulation in Experiments 1 and 2 have weakened WBCEs and gyres in comparison to the Control with both WBCEs showing reductions in time-averaged mass transports. This is expected since the derived wind-stress anomalies act to increase the wind-forcing asymmetry into the gyres from a relatively straight east-west forcing (Fig. 2c) to more convoluted forcings (Fig. 6). Weakened WBCEs due to increased wind-forcing asymmetry have been noted by Rhines and Schopp (1991), particularly near the jet separation region (Hogg et al. 2009). The wind-stress forcing in Experiment 1 leads to a weakened subtropical recirculation zone, strengthened subpolar recirculation zone and a poleward change in the jet-axis tilt. The change in jet-axis tilt is consistent with findings by Moro (1988) and Hogg et al. (2005), where the dominant recirculation zone, in this case the subpolar recirculation, pulls the WBCE towards itself. This slightly counter-intuitive result is because the dynamics of the ocean gyres is controlled by the PV circulation, rather than inertia. The entire WBCE is also shifted poleward in comparison to the benchmark circulation except near the  Table 2). Top row shows forcing anomalies for Experiment 1 and bottom row shows forcing anomalies for Experiment 2. a, d Ekman pumping anomaly. b, e Diabatic entrainment anomaly (sign-flipped to show upper-isopycnal). c, f Upper-ocean total entrainment. Black dashed line indicates timeaveraged position of inter-gyre boundary for the respective modelled double-gyre circulation jet separation region. These changes in the circulation are largely consistent with findings in KB22 where changes in strength of the recirculation zones and a poleward shift in the WBCE were also observed. The anomalous circulation shown for the upper-and middle-isopycnal layers (Fig. 7b, f) are almost identical in structure and show an asymmetric tripole anomaly along the WBCE and recirculation zones. The largely barotropic structure of this anomaly indicates downwards transfers of momentum from the upper-isopycnal layer in the form of eddy form stresses (McWilliams 2008). However, the response of the middle-isopycnal layer is such that the subpolar recirculation zone now has a greater vertical shear, while the subtropical recirculation zone has a lesser vertical shear (in comparison to Fig. 4d, e) This effect will be shown to have an important consequence on the baroclinicity of the inertial recirculations and vertical momentum transfers later on. Similar patterns as modes of decadal ocean gyre variability have been noted by Hogg et al. (2005) Berloff et al. (2007, etc.) and have been described more simply as 'meridional shift' modes. We will describe this mode as a combination of relative recirculation zone strength change, leading to a change in jet-axis tilt, and a separate poleward meridional shift of the entire WBCE. In our situation, this more detailed decomposition of the asymmetric tripole circulation anomaly is necessary as preliminary analysis from KB22 indicates the changes associated with the relative recirculation zone strength and WBCE meridional shift are controlled by distinct mechanisms.
In a similar but opposite fashion, wind-stress forcings in Experiment 2 lead to a weakening of the subpolar recirculation zone, strengthening of the subtropical recirculation zone, a consequent equatorward shift in the jet-axis tilt, and a equatorward meridional shift of the WBCE in the upperisopycnal layer. The WBCE also appears to be now deflected  Fig. 4d and e, respectively. Panels (b, d, f, h) are saturated to show the weaker circulation anomalies in the inertial recirculations. a, b Upper-isopycnal time-averaged transport streamfunction for Experiment 1 and anomaly from Control. c, d Upper-isopycnal time-averaged transport streamfunction for Experiment 2 and anomaly from Control. e, f Middle-isopycnal time-averaged transport streamfunction for Experiment 1 and anomaly from Control. g, h Middle-isopycnal time-averaged transport streamfunction for Experiment 2 and anomaly from Control more strongly equatorward as it separates from the western boundary. In the middle-isopycnal layer, the response is such that the subpolar recirculation zone now has a lesser vertical shear and the subtropical recirculation has a greater vertical shear (in comparison to Fig. 4d, e).
There are certainly differences in the modelled ocean circulation response in comparison to KB22. This includes changes in strength of the subpolar recirculation zone and gyre, along with greater changes in the jet-axis tilt. However, this is not necessarily a problem since both east-west dipole vorticity flux anomalies are still well modelled and generic changes in the inertial recirculations and meridional shifts of the WBCE are still captured. This means that any findings in this paper, including the responsible mechanisms, are still of value and applicable to the findings of KB22. Figures 8 and 9 show effects of wind-stress anomalies on the PV anomaly distribution in the gyres. In this case, we have separated out PV anomalies into relative vorticity, which are strongest along the WBCs and WBCE, and buoyancy components, which are strong over the entire ocean, except near the western boundary. Figures 8b, d and 9b, d reveal the dominant effect on the upper-isopycnal PV is the generation of time-averaged, asymmetric, tripole anomalies over the WBCE and recirculation zones. The time-averaged buoyancy anomalies are roughly 2-3 times stronger than the accompanying relative vorticity anomalies over the same region. The central anomaly of these tripole patterns are largely due to shifts in the location of the inter-gyre boundary, while the opposite-signed anomalies that lie adjacent on either side are associated with circulation changes in the inertial recirculations. Figure 8d shows that there is increased baroclinicity in the subpolar recirculation zone, while there is decrease baroclinicity in the subtropical recirculation zone. Similar but opposite changes may be observed in Fig. 9d. The response of the inertial recirculation zones is a direct consequence to the local anomalous forcing that sits over the region, namely the western basin wind-stress curl anomaly. Changes in the baroclinicity of the inertial recirculations allow for them to more or less efficiently transfer momentum into the lower layers through baroclinic instabilities. The positive time-averaged buoyancy anomaly in the eastern basin (Fig. 9d) appears to be connected to the central buoyancy anomaly along the WBCE. This indicates a nonlocal ocean gyre response but needs to be confirmed with further analysis (see Sect. 4.5). Buoyancy anomalies in the middle-isopycnal layer show similar but opposite-signed anomalies along the recirculation zones and WBCE. There are also time-averaged buoyancy anomalies in the eastern basin on either side of the WBCE associated with the gyre recirculations of eastern basin vorticity fluxes. Such an anomaly is not visible in the upper-isopycnal layer which indicates the ocean response to eastern basin forcing is baroclinic in nature. This is important since if the eastern basin wind-stress curl in the upper-isopycnal was being advected by the gyre recirculations, then we should be able to observe this, like we have done so in the middle-isopycnal layer. We note that we do not observe PV anomalies similar to those in Fig. 1f. This is likely due to inaccuracies in the modelled diabatic entrainments from idealised formulations of thermal forcing in the model.

Time-averaged PV response to wind-stress anomalies
Outside of the WBCE regions in the open ocean, relative vorticity anomalies only appear significant very close to the jet separation region where local modification of the WBL occurs. For example in Fig. 8b, relative vorticity anomalies in the jet separation region are due to poleward deflections of the subtropical WBC. A similar pattern may be seen in Fig. 9b for poleward deflections of the subtropical WBC. Changes in the relative vorticity profile along the remainder of the western boundary are almost negligible, which agree with our findings that the measured changes in ocean gyre circulation are not triggered by processes in the WBL (e.g. Cessi et al. 1987), but more likely by reorganisations of the WBCE and inertial recirculations due to anomalous wind-forcing.

Existence of waves and importance of eastern basin forcing
Although a local ocean gyre response may be described for the wind-stress curl anomaly in the western basin, it is unclear from the two double-gyre experiments alone whether the wind-stress curls in the eastern basin are important. For example, we are still not sure whether the meridional shift of the WBCE is part of the inertial recirculation zone response or whether they are related to the eastern basin forcing.
We check for the presence of westward-propagating signals from the eastern basin by plotting time-longitude plots of PV anomalies (see Fig. 10). This is to see if the eastern basin forcing is able to generate a non-local, baroclinic Rossby wave adjustment (Sasaki and Schneider 2011;Sasaki et al. 2013). We change coordinate systems such as we follow the PV along the time-averaged position of the WBCEs, i.e. the inter-gyre boundaries, in each experiment. Both experiments show there are clearly westward-propagating Rossby waves along the WBCE with wave formation in the eastern basin. There is a noticeable increase in concentration of PV as the Rossby waves enter the WBCE region. Such an intensification is consistent with the narrowing of the meridional scale of the wave as it enters the WBCE region. This is because as the wave narrows in structure, it increases in PV concentration in order to conserve PV (Sasaki et al. 2013). Rossby wave phase speed is approximated using a basin crossing time of ∼ 1.5-years, which yields a propagation speed of c R = 10 cm s −1 . These speeds are around twice as fast as those predicted by Dewar (2003) and Sasaki et al. (2013). This is likely due to quasilinear interaction of the Rossby wave with the mean flow, which increases wave speed propagation (Dewar and Morris 2000).
Note that although it is clear that westward-propagating wave behaviour is visible along the WBCE, the signal is strongly influenced by the powerful effects of advection. Indeed, many of the propagating waves appear to be destroyed by effects of WBCE advection, which may be seen by the presence of Rossby deformation scale standing waves ( Fig. 10), particularly near the western boundary where the WBCE is strongest. These regions also show strong, intrinsic, interannual variability of the WBCE that occurs due to vigorous mesoscale activity. This means that although Rossby wave signal travels at a largely predictable speed, their arrival at the western boundary is not guaranteed. This leads to reduced predictability of wind-induced ocean gyre variability due to strong effects of time-averaged advection and intrinsic interannual variability of the WBCE (Nonaka et al. 2016(Nonaka et al. , 2020.
Finally, to check that eastern basin vorticity flux anomalies do indeed induce meridional shifts of the WBCE, we run another two additional sets of double-gyre simulations with modified forcings from Experiments 1 and 2. These modified wind-stress forcings are identical except that the eastern basin vorticity flux anomaly is weakened by setting the eastern basin zonal wind-stress anomalies in Fig. 1b to zero. Note that it is not possible to set meridional wind-stresses in the eastern basin to zero without creating a discontinuity in the wind-stress forcing, which may have undesirable effects on the modelled double-gyre circulation. Hence, we have elected to only set the zonal wind-stress anomaly to zero which does not cause this problem. These modified Experiments 1 and 2 will be referred to as supplementary Experiments 3 and 4 forcings, respectively. Figure 11a, b, d, e show the modelled time-averaged Ekman pumping and upper-ocean total entrainments under the reduced eastern basin wind-stress forcings. It is easy to see that setting the eastern basin zonal wind-stress anomaly in Fig. 1b to zero does indeed reduce both the time-averaged Ekman pumping and upper-ocean total entrainments in the region by ∼ 50% . Furthermore, western basin forcings appear almost identical to those shown in Fig. 6a, c, d, f. Hence, it is reasonable to assume that any observed changes between the modelled circulations in Figs. 11c, f and 7a, c, respectively, are purely due to reduced wind-stress curls in the eastern ocean basin.
Indeed, the reduced eastern basin wind-stress curl anomaly leads to an equatorward meridional shift of the inter-gyre boundary in Fig. 11c and a poleward meridional shift of the inter-gyre boundary in Fig. 11f. The strengths of the inertial recirculation zones and axis tilt of the WBCE remain identical, or very similar, to Fig. 7a, c. This agrees with our hypothesis that the western basin vorticity flux anomaly is responsible for changes in strength of the inertial recirculations. Furthermore, and more interestingly, this also implies that the eastern basin wind-stress anomaly is responsible for independently controlling meridional shifts of the inter-gyre boundary. The direction of the meridional shifts in Fig. 11b, c are consistent with the baroclinic Rossby wave adjustment mechanism proposed by Sasaki et al. (2013), with an anti-cyclonic wind-stress curl leading to a poleward meridional shift of Fig. 10 Presence of nonlinear, baroclinic Rossby waves along the WBCE. Data is taken over a 6-year period for the system in statistical equilibrium. a PV anomalies along the time-averaged position of the WBCE in Experiment 1. b PV anomalies along the time-averaged position of the WBCE in Experiment 2. Approximate trajectory of wave is plotted in both panels with a phase speed c R = 10 cm s −1 the WBCE. The magnitude of the shift observed each panels is around 75 km which is slightly less than the meridional shifts measured in KB22. This is expected since the eastern basin wind-curls are only reduced in magnitude and not set to zero.

Time-averaged response of PV sources and sinks in the ocean gyres
The time-averaged, gyre-integrated PV budgets for Experiment 1 in the upper-and middle-isopycnal layers are given in Table 3. In the upper-isopycnal layer, there is a 7-8% drop in the wind-induced vorticity fluxes for each gyre. This is expected since the forcing asymmetry and jet-axis tilt both increase due to changes in the relative recirculation zone strengths. To counteract this, the inter-gyre PV flux in the gyres falls by 20% and the viscous boundary fluxes falls by 1-4%. In this experiment, the inter-gyre PV flux handles the vast majority of the decrease in wind-induced vorticity flux, while the viscous boundary fluxes only account for a small portion. This is consistent with our findings that the western boundary layer is only adjusted near the WBC separation region, so changes in the viscous boundary flux will consequently be small. The changes in the PV budget were found to be statistically significant under a one-tailed Student's t-test with a 95% confidence interval. Changes in the PV tendency were not found to be statistically significant. Similar but smaller changes were observed for Experiment 2 in the upper-isopycnal layer (see Table 4). The windinduced vorticity fluxes dropped by 2-7%, inter-gyre PV fluxes fell by 5% and viscous boundary fluxes by 4-9%. In this situation however, changes in the inter-gyre PV flux are also no longer statistically significant under the same null hypothesis along with the change in tendency. PV budgets changes in the middle-isopycnal layer appears to largely mirror the upper isopycnal layer except that the forcing is smaller due to only diabatic entrainment being present, which is further weakened by the larger depth of the isopycnal layer.
Overall, our analysis shows that increased wind-forcing asymmetry decreases vorticity fluxes into the ocean gyres. The response of the ocean gyre circulation is consistent with a reduction in PV sinks across the gyres. This is achieved either through a drop in inter-gyre PV flux and/ or a drop in viscous boundary fluxes where local restructuring of the WBL near the WBC separation region is able to reduce normal relative vorticity gradients. We note that although the PV budgets are useful in capturing basin-integrated changes in the ocean gyres, more location-dependent circulation changes such as responses to western and eastern basin forcings are not well captured.

Summary of time-averaged ocean gyre responses to east-west vorticity flux anomaly dipoles
The ocean gyre response was found to be different according to the location of the vorticity flux anomaly in the ocean basin. Wind-stress curl anomalies in the western basin were associated with changes in relative strength of the inertial recirculation zones, which consequently affected the jet-axis tilt. For example, a positive windstress curl in the western basin of upper-isopycnal layer was found to weaken the subtropical recirculation zone and strengthen the subpolar recirculation zone. This led to a relatively stronger subpolar gyre which pulled the WBCE into a poleward tilt. This increases the baroclinicity of the subpolar recirculation zone and decreases the baroclinicity of the subtropical recirculation zone. This allows the subpolar recirculation zone to more efficiently transfer excess momentum into the lower layers due to increased anomalous PV flux. A similar but opposite pattern may be observed for the subtropical recirculation zone. The wind-stress curl anomalies in the eastern basin on the other hand were found to affect the meridional position of the WBCE. Our results with reduced wind-stress forcings in the eastern basin led to reductions in the magnitude of the WBCE meridional shift, while the recirculation zone strengths remained unchanged. The presence of westward-propagating waves from the eastern boundary which travel along the WBCE confirmed the existence of baroclinic Rossby waves likely responsible for the meridional shift. These findings imply the presence of two distinct mechanisms of decadal ocean gyre variability (Dewar 2003;Sasaki et al. 2013). Furthermore, the response of the ocean gyres seems to depend on the spatial location of the forcing.

Discussion
Part II of this study investigated the time-averaged response of the ocean gyres to wind-stress curl anomalies found in Part I of this two-paper study (see KB22). These modes of atmospheric variability arise as a consequence of large-scale, standing Rossby wave disturbance growth and are hence relevant to the midlatitude climate variability. Although analysis in KB22 found a statistical link between these modes of atmospheric variability and the delayed ocean gyre responses, it did not reveal enough about the nature of the dynamics and physical processes involved. Indeed, the complexity of the involved mode of atmospheric variability, which contains zonal asymmetries, and the two seemingly uncorrelated responses in the ocean gyres warranted further investigation.
The main finding of this paper is the spatially-inhomogeneous response to wind-stress curl anomalies in the ocean basin. That is, wind-stress curl anomalies in the western and eastern ocean basins generate distinct responses in the modelled double-gyre circulation. Although these two observed mechanisms of decadal ocean gyre variability are not new (see Dewar 2003;Hogg et al. 2005;Berloff et al. 2007;Sasaki and Schneider 2011;Sasaki et al. 2013, for examples), the finding that both of these mechanisms act largely independently in the same modelled circulation is, to our knowledge, a new one. 1 The findings from our reduced eastern basin wind-stress curl experiments and the existence of Rossby wave propagation along the WBCE is crucial in showing that these two aforementioned mechanisms are distinct. For example, if it was true that the eastern basin forcing was important in determining recirculation zone strengths, then this should have been observed when the eastern basin forcing was reduced. Furthermore, the existence of wave-like behaviours along the WBCE provided a means for the ocean gyres to respond to anomalous forcing, which acts largely independently of the advection-dominated inertial recirculation zones. Our observation of propagating waves is potentially problematic with findings by Dewar (2003) where adjustment of the ocean gyres to anomalous wind-induced vorticity fluxes was unlikely to be caused by Rossby wave propagation. However, we believe this is likely due to the relatively small longitudinal ocean basin dimension (less than half the size of ours), which generates a double-gyre circulation with a WBCE that nearly reaches the eastern boundary. This means that any wind-induced Rossby waves are quickly broken up by the turbulence and powerful timeaveraged advection of the WBCE. The existence of Rossby waves in our modelled double-gyre circulation is nevertheless still rather remarkable given the fact that our modelled WBCE is both stronger in its time-averaged transport and in its turbulent behaviour.
In summary, we believe the spatially inhomogeneous response of the ocean gyres to large-scale wind-stress anomalies depends on properties of the time-averaged circulation itself. In regions where the time-averaged circulation is strong, i.e. the inertial recirculations and WBCE, the response is dominated by advective effects and changes in the time-averaged circulation of those regions. Conversely, the presence of strong time-averaged circulation provides unfavourable conditions for the formation of wind-induced Rossby waves, which are quickly broken up by the powerful WBCE.
On the other hand, in regions where the time-averaged circulation is weak, i.e. in the eastern ocean basin, the conditions are now favourable for Rossby wave formation due to the largely stagnant flow in this region. Again, conversely, these regions also provide unfavourable conditions for advective processes to respond to anomalous forcing. Indeed, the only possible advective response is by gyre recirculations which are unlikely to be able to generate the measured WBCE response. This is because mixing and stirring effects dominate and lead to the quick loss of PV memory by fluid parcels (Rhines and Schopp 1991; Kurashina et al. 2021).
Although these modelled ocean gyre responses show some differences with KB22, our findings remain highly relevant and applicable. For example, the fact that our two upper-isopycnal ocean gyre responses are governed by distinct mechanisms is consistent with the lack of temporal correlations in the ocean gyre response modes from KB22. Furthermore, the reduced predictability of these ocean gyre responses from KB22 is well explained by the finding that both of these mechanisms involve regions of the ocean gyres where the effects of mesoscale turbulence are strong.
Appendix: Overactive convection in the subpolar gyre Figure 12 shows the effect of overactive convection on the time-averaged, upper-ocean total entrainment. Extremely large values of entrainments due to persistent convection events were found near the subpolar WBC separation region, where extreme lows in SSTs develop. Convection events are not necessarily a problem in this model if similar values were found in the coupled configuration of this model. However, we found zero events of ocean convection in KB22, leading us to conclude that such values of entrainment are aphysical. Such large forcings of diabatic entrainment near the jet separation region led to weakening of the ocean gyre circulation. These modelled circulations were akin to those computed by Hogg et al. (2009), where strong Ekman pumping anomalies near the WBC separation region destabilised the WBCE. Furthermore, vorticity flux anomalies due to convection were found to be of similar magnitude to those due to wind-stress anomalies in Sect. 4. This would have made it difficult to distinguish between circulation anomalies due to changes in the wind-stress field, and circulation anomalies due to overactive convection in the model. The fact that entrainments over the remainder of the gyres are at normal levels indicates that thermal forcing over the subpolar WBC is likely the cause of the issue. Thermal forcing for an ocean-only configuration of Q-GCM is taken from the time-averaged ocean diabatic heating obtained from a fully coupled model run. In our case, the time-averaged diabatic heating from KB22 was used. If the position of the WBCs misaligns with the thermal forcing, which only occupies a small region over the basin (see Fig. 2d), then extreme values of surface temperature may develop due to the lack of thermal regulation by the ocean diabatic heating. In the subtropical gyre, surface temperatures are also regulated by a fixed temperature Dirichlet boundary condition. However, the subpolar gyre has no such temperature regulation along the northern boundary. In order to remove the effects of convection, temperature is capped below at the fixed temperature of the upper-isopycnal layer T 1 but we choose to not model the adjustment in entrainment which accompanies this. For example, whenever T m falls below T 1 , i.e. the ocean is unstably stratified, there is usually an adjustment of the ocean diabatic entrainment term by adding on where Δt is the ocean time-step and Δ 1 T = T 1 − T 2 . At the end of the convection step, T m is set to T 1 . In our case, we have removed this entrainment 'correction' term and simply only implemented the temperature capping at the end of the step. This essentially means we have applied an artificial thermal forcing over regions where overly extreme values of low temperatures may develop. Importantly, the model still conserves PV which is important since our analysis relies on this for accurate PV budgets. We recommend a more robust diabatic forcing implementation for future studies using this particular model configuration, which allows for timevarying surface heat fluxes against prescribed atmospheric surface temperatures. Availability of data and materials Not applicable.

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 need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 12
Overactive convection in the subpolar WBC. Time-averaged upper-ocean total entrainment e total with contour plots saturated by a factor of 65. Data was averaged over a 40-year period for a system in statistical equilibrium