2.1. Quasi-geostrophic equations

The ocean model we used was identical to that used by Berloff (Reference Berloff2015). We considered an idealised, eddy-resolving QG model with three stacked isopycnal layers. It was designed to mimic the wind-driven, midlatitude, double-gyre ocean circulation. The ocean basin was configured as a flat-bottomed square box with $(x,y)$ as zonal and meridional coordinates. The square basin $D(x,y)$ had length $2L=3840\ \textrm {km}$ and depths $H_k=250\ \textrm {m}$, $750\ \textrm {m}$, $3000\ \textrm {m}$ for $k=1,2,3$ with $-L\leqslant x\leqslant L$ and $-L\leqslant y \leqslant L$. The governing QG equations are given in terms of the streamfunction $\psi _i$ and PV anomalies $q_i$, which are defined by

(2.1a)$$\begin{gather} q_1=\nabla^2\psi_1-S_{11}(\psi_1-\psi_2), \end{gather}$$

(2.1b)$$\begin{gather}q_2=\nabla^2\psi_2-S_{21}(\psi_2-\psi_1)-S_{22}(\psi_2-\psi_3), \end{gather}$$

(2.1c)$$\begin{gather}q_3=\nabla^2\psi_3-S_{31}(\psi_3-\psi_2) , \end{gather}$$

where the stratification parameters $S_1, S_{21}, S_{22}, S_3$ are defined by the first and second baroclinic Rossby deformation radii, chosen to be $40.0\ \textrm {km}$ and $20.6\ \textrm {km}$, respectively. This then allows us to formulate the QG equations:

(2.2a)$$\begin{gather} \frac{\partial q_1}{\partial t}+J(\psi_1,q_1)+\beta \frac{\partial \psi_1}{\partial x}=\frac{1}{\rho_1 H_1}W+\nu \nabla^4\psi_1 , \end{gather}$$

(2.2b)$$\begin{gather}\frac{\partial q_2}{\partial t}+J(\psi_2,q_2)+\beta \frac{\partial \psi_2}{\partial x}=\nu \nabla^4\psi_2 , \end{gather}$$

(2.2c)$$\begin{gather}\frac{\partial q_3}{\partial t}+J(\psi_3,q_3)+\beta \frac{\partial \psi_3}{\partial x}={-}\gamma\nabla^2\psi_3+\nu \nabla^4\psi_3 , \end{gather}$$

where $\beta =2\times 10^{-11}\ \textrm {m}^{-1}\ \textrm {s}^{-1}$ is the planetary vorticity gradient, $\rho _1=1000\ \textrm {kg}\ \textrm {m}^{-3}$ is the upper-layer density, $\nu =20\ \textrm {m}^2\ \textrm {s}^{-1}$ is the eddy viscosity coefficient, $\gamma =4\times 10^{-8}\ \textrm {s}^{-1}$ is the bottom friction parameter and $J(\cdot ,\cdot )$ is the Jacobian operator.

The asymmetric wind forcing was fixed in time and defined as

(2.3) \begin{equation} W(x,y) = \left\{\begin{array}{@{}ll} -\dfrac{{\rm \pi}\tau_0A}{2L}\sin\left[\dfrac{{\rm \pi}(L+y)}{L+Bx}\right], & \textrm{if }y\leqslant Bx, \\ +\dfrac{{\rm \pi}\tau_0}{2LA}\sin\left[\dfrac{{\rm \pi}(y-Bx)}{L-Bx}\right], & \textrm{if }y> Bx. \end{array} \right. \end{equation}

where $A=0.9$, $B=0.2$ and $\tau _0=0.08\ \textrm {Nm}^{-2}$. The extra factor of $1/2$ in (2.3) and an order of magnitude correction in $\tau$ are both rectifications of typos by Berloff (Reference Berloff2015). Note that $A$ and $B$ are parameters that control the degree of asymmetry in the wind forcing. Partial-slip boundary conditions $\psi _{\boldsymbol {nn}}=\alpha \psi _{\boldsymbol {n}}$, where $\boldsymbol {n}$ is the normal-to-wall unit vector, were enforced on the lateral boundaries (Haidvogel et al. Reference Haidvogel, McWilliams and Gent1992). Mass conservation constraints were applied for each layer (McWilliams Reference McWilliams1977). The partial-slip boundary condition parameter $\alpha =120\ \textrm {km}$ was chosen such that it was close to a free-slip boundary condition but still allowed for momentum flux through the lateral boundaries (Deremble et al. Reference Deremble, Hogg, Berloff and Dewar2011). The equations were solved using the CABARET scheme (Karabasov, Berloff & Goloviznin Reference Karabasov, Berloff and Goloviznin2009) with a nominal resolution of $7.5\ \textrm {km}$, which corresponded to a uniform $513^2$ grid. The model was integrated for 60 model years with a 20-year spin-up period. Data were accumulated every model day.

We present typical snapshots and time-mean solutions of system (2.2) in figures 1 and 2, respectively. Wind acting on the upper isopycnal layer created the double-gyre flow regime through the imposed large-scale Ekman pumping anomalies. The two gyres were separated by the powerful eastward jet, which was supported by the subpolar and subtropical WBCs. Strong eddy activity was visible in the upper isopycnal layer around the WBC separation point and eastward jet, while the lower layers were dominated by weaker eddies. The middle isopycnal layer also showed a region of uniform PV, which was homogenised by eddies (Rhines & Young Reference Rhines and Young1982; Lozier & Riser Reference Lozier and Riser1989).

Figure 1.

Instantaneous snapshots of the statistically equilibrated reference solution of system (2.2). (a–c) Velocity streamfunctions in the upper, middle and lower isopycnal layers, respectively. (d–f) PV anomalies in the upper, middle and lower isopycnal layers, respectively. For presentation purposes, the colour-bar range for each plot is determined by the maximum and minimum values. This is repeated for all consequent plots unless stated otherwise. Values of Max for each panel are stated for comparison. (a) $\textrm {Max} = 7.3\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (b) $\textrm {Max} = 3.3\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (c) $\textrm {Max} = 2.4\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (d) $\textrm {Max} = 1.1\times 10^{-4}\ \textrm {s}^{-1}$. (e) $\textrm {Max} = 2.5\times 10^{-5}\ \textrm {s}^{-1}$. (f) $\textrm {Max} = 1.5\times 10^{-5}\ \textrm {s}^{-1}$.

Figure 2.

Time-mean solutions of system (2.2) in statistical equilibrium. (a–c) Velocity streamfunctions in the upper, middle and lower isopycnal layers, respectively. (d–f) PV anomalies in the upper, middle and lower isopycnal layers, respectively. Inter-gyre boundary, shown in (a), is defined by the time-mean contour emanating from the western boundary. Values of Max for each panel are stated for comparison. (a) $\textrm {Max} = 6.0\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (b) $\textrm {Max} = 2.3\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (c) $\textrm {Max} = 9.7\times 10^3\ \textrm {m}^2\ \textrm {s}^{-1}$. (d) $\textrm {Max} = 7.1\times 10^{-5}\ \textrm {s}^{-1}$. (e) $\textrm {Max} = 1.7\times 10^{-5}\ \textrm {s}^{-1}$. (f) $\textrm {Max} = 5.5\times 10^{-6}\ \textrm {s}^{-1}$.

2.2. Linearised dynamics

As we were interested in studying advection-induced anomalies, we also required solutions for the linearised dynamics. System (2.2) was linearised around the state of rest to give

(2.4a)$$\begin{gather} \frac{\partial q_1}{\partial t}+\beta \frac{\partial \psi_1}{\partial x}=\frac{1}{\rho_1 H_1}W+\nu \nabla^4\psi_1 , \end{gather}$$

(2.4b)$$\begin{gather}\frac{\partial q_2}{\partial t}+\beta \frac{\partial \psi_2}{\partial x}=\nu \nabla^4\psi_2 , \end{gather}$$

(2.4c)$$\begin{gather}\frac{\partial q_3}{\partial t}+\beta \frac{\partial \psi_3}{\partial x}={-}\gamma\nabla^2\psi_3+\nu \nabla^4\psi_3 . \end{gather}$$

The parameters for these equations, as well as the configuration of the numerical CABARET solver, were identical to system (2.2). These equations are also studied analytically in Appendix A by generalising the Munk boundary layer solution (Munk Reference Munk1949) for partial-slip boundary conditions.

Figure 3 shows that the Sverdrup (i.e. linear) gyres and WBCs in the upper isopycnal layer were present but differed from the nonlinear dynamics solutions. Furthermore, nonlinear phenomena, such as the eastward jet and recirculation zones, were missing. Contours plots for the lower layers were omitted as these only show the gravest basin modes propagating across the basin. This is because the eddy forcing that drives the lower layers (Holland Reference Holland1978) was no longer present.

Figure 3.

Contour plots of time-mean solutions for system (2.4) in the upper isopycnal layer. Values of Max for each panel are stated for comparison. (a) Velocity streamfunction, $\textrm {Max} = 7.0\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (b) PV anomaly, $\textrm {Max} = 3.7\times 10^{-4}\ \textrm {s}^{-1}$. Inter-gyre boundary in (a) is defined in the same manner as discussed in the caption of figure 2.

2.3. Lagrangian particles

In the subsequent sections, we use a Lagrangian particle analysis to identify inter-gyre particle pathways. Such techniques are useful owing to their ability to capture ensemble behaviour of fluid parcels, and hence, long-range material transports. For example, through the Lagrangian framework, we are able to study the long-range transport of PV, to which the Eulerian framework is not well suited.

Simulations were run by releasing an ensemble of $N$ particles into the upper isopycnal layer and solving

(2.5)\begin{equation} \frac{\textrm{d}\boldsymbol{x}_j(t)}{\textrm{d}t}=\boldsymbol{u}(t,\boldsymbol{x}_j(t)),\quad \textrm{for } j=1,\ldots,N , \end{equation}

where $\boldsymbol {x}_j(t)$ is the position of particle $j$ at time $t$ and $\boldsymbol {u}(t,\boldsymbol {x})=(u(t,\boldsymbol {x}), v(t,\boldsymbol {x}))$ is the velocity of the flow at position $\boldsymbol {x}$ and time $t$. System (2.5) was solved offline using the 4th-order Runga–Kutta scheme with a two-dimensional (2-D) cubic spatial interpolation and a one-dimensional (1-D) cubic temporal interpolation. Total PV $(Q=q+\beta y)$ was estimated along the trajectories by using a 2-D cubic interpolation (Total PV was used for Lagrangian particle analysis rather than PV anomalies because it is a materially-conserved quantity.).

3.1. Nonlinear western boundary layer

The differences in structure of the upper-layer WBL were particularly significant when inspecting figures 2 and 3. Figure 3(b) shows much stronger zonal PV gradients near the western boundary, which are accompanied by increased meridional velocities in the region. Typical time-mean meridional velocities for system (2.2) were $1\ \textrm {ms}^{-1}$, while they were up to $4\ \textrm {ms}^{-1}$ for system (2.4). This implied a cutting down of the high velocity and vorticity in the WBL through nonlinear effects. Indeed, Veronis (Reference Veronis1966a) also showed through a perturbation analysis that this restructuring of the WBL is likely owing to the meridional inhomogeneities in the structure of the wind forcing and ocean gyres. Furthermore, they showed that this restructuring is controlled by effects of PV advection. However, it is unclear what processes maintain the structure of the WBL after the assumptions of the stability analysis break down. A full analysis of the advection term in (2.2a) is required to understand this, which is outside the scope of this study. Reduction of WBC velocities also indicates that the nonlinear dynamics solutions will also show a reduction in viscous relative vorticity fluxes through the western boundary in comparison to the linear dynamics solutions (Spall Reference Spall2014, see also § 3.3). We now continue our analysis by defining the CGAs through a solution decomposition.

3.2. Solution decomposition

The CGAs, as well as other nonlinear circulation features, are defined by decomposing the nonlinear dynamics solutions into the linear dynamics (subscript ‘lin’), time-mean advection-induced nonlinear anomalies (subscript $\oplus$) and transient fluctuations (primed):

(3.1)\begin{equation} \psi_i = \bar{\psi}_{i,{lin}} + \bar{\psi}_{i,\oplus} + \psi_i' , \end{equation}

for layers $i=1,2,3$ (see Shevchenko & Berloff (Reference Shevchenko and Berloff2016) for identical decomposition). Note that the time-mean streamfunction is given by $\bar {\psi }_i=\bar {\psi }_{i,{lin}} + \bar {\psi }_{i,\oplus }$. This process was also repeated for PV anomalies. We chose this particular decomposition as we wanted to highlight the regions where the discrepancy between the nonlinear dynamics solutions, i.e. system (2.2), and linear dynamics solutions, i.e. system (2.4), was the greatest. This then allowed us to pinpoint the regions where nonlinear effects were particularly strong, which will guide our analysis in later sections. Comparisons of the nonlinear dynamics against the linear dynamics have been used in the past to study recirculation zones and eastward jets (e.g. Veronis Reference Veronis1966a; Harrison & Stalos Reference Harrison and Stalos1982), and we will be using the same technique for the CGAs.

The time-mean advection-induced anomalies (figure 4) revealed the recirculation zones on either side of the eastward jet, and the sharp zonal PV anomaly gradients near the western boundary that extended out into the gyre interiors. The CGAs were seen to be embedded in the subpolar (subtropical) gyres as opposite-signed, anti-cyclonic (cyclonic) anomalies. We stress here that because the CGAs were defined as a deviation from the linear dynamics solutions, they would not necessarily be visible in observations/models without this particular solution decomposition. However, because the decomposition was available to us, we will show that these anomalies are dynamically important and linked to mechanisms that control the double-gyre circulation.

Figure 4.

Contour plots of time-mean advection-induced anomalies in the upper isopycnal layer for the reference, statistically equilibrated flow regime. Values of Max for each panel are stated for comparison. (a) Velocity streamfunction, $\textrm {Max}=6.0\times 10^4\ \textrm {m}^2\ \textrm {s}^{-1}$. (b) PV anomaly, $\textrm {Max}/5=6.5\times 10^{-5}\ \textrm {s}^{-1}$. See figure 5 for zoomed-in plots of PV anomalies near the western boundary. Use of colour-bar range for (a) is the same as that in figure 1, but the colour-bar range for (b) is limited to 1/5 of the maximum value to better show the anomalies in the ocean interior, as they are weak compared with the anomalies in the WBL. Box surrounded by black dotted lines indicates the region near the western boundary plotted in figure 5.

Upon further inspection of these anomalies, we see that they were, in fact, two separate lobes (see figure 5c). The first lobe, which was by far the stronger anomaly, was situated within the viscous sublayer of the WBLs. We will refer to this lobe as the ‘WBL lobe’ for brevity. The second lobe, which was a much weaker but larger anomaly, was spread throughout the western half of each gyre interior. This lobe will be referred to as the ‘interior lobe’. Why the CGAs were separated into two separate lobes and the influences of nonlinearity in this regions on the global circulation is still not fully known.

Figure 5.

Contour plots of time-mean PV anomalies in the upper isopycnal layer near the western boundary. Figures are zoomed-in to show the region outlined in figure 4. Values of Max for each panel are stated for comparison. (a) Nonlinear dynamics solution, $\textrm {Max} = 7.1\times 10^{-5}\ \textrm {s}^{-1}$. (b) Linear dynamics solution, $\textrm {Max} = 3.7\times 10^{-4}\ \textrm {s}^{-1}$. (c) Time-mean advection-induced anomalies, $\textrm {Max} = 3.2\times 10^{-4}\ \textrm {s}^{-1}$. Areas that contain the subtropical WBL and interior lobes of the CGAs are marked by dotted lines in (c). Note that the interior lobe extends further outside of (c) into the ocean interior.

For the lower layers, the equivalent plots for figure 4 are shown in figure 2(b,c,e,f). This is because the lower layers were essentially stagnant for the linear dynamics solutions. Inspecting these plots, we observed that the recirculation zones were visible in the lower layers but no CGA-like patterns were present. Thus, we concluded that the CGAs were confined to the upper isopycnal layer, which was the only layer driven directly by wind-stress curl. This indicated that wind-stresses are important in the emergence of the CGAs.

Shevchenko & Berloff (Reference Shevchenko and Berloff2016) found that by strengthening the CGAs through an external forcing term, the eastward jet was weakened. However, they did not analyse in detail the dynamics associated with them. We continue our analysis of the CGAs by examining the relative-vorticity fluxes at the western boundary, where the time-mean advection-induced anomalies are strongest.

3.3. Excess PV buildup

To better understand how nonlinear effects adjust the large-scale circulation in the upper isopycnal layer, we calculated the PV budgets for both the linear and nonlinear dynamics solutions. A similar budget was computed by Berloff et al. (Reference Berloff, Hogg and Dewar2007b), albeit with different model configurations. The sources and sinks of PV in the upper isopycnal layer are considered below and PV fluxes are estimated and given in table 1.

Table 1.

PV budgets for linear and nonlinear dynamics solutions for the subtropical gyre in the upper isopycnal layer. PV sources are estimated using time-mean values with units $\textrm {m}^2\ \textrm {s}^{-2}$. Bracketed values are identical computations made for the subpolar gyre. Terms that account for ${<}0.05\,\%$ of the budget are deemed as negligible. We find that the PV budget is largely insensitive to small changes made in the inter-gyre boundary.

To compute a PV budget for each gyre, we need an inter-gyre boundary. This was calculated by following the time-mean contour emanating from the western boundary, which extended eastward across the basin. In the linear dynamics solutions, the inter-gyre boundary was a straight line dividing the two gyres (figure 3a). However, in the nonlinear dynamics solutions, the inter-gyre boundary was given by the time-mean eastward jet position (figure 2a). Note that the time-mean contour was particularly difficult to compute near the eastern boundary, where the eastward jet extension tore itself apart. However, we found that small changes to this inter-gyre boundary had negligible effects on the results from the subsequent PV budget.

The only source of PV for the gyres was through Ekman pumping anomalies generated by wind-curl. This was calculated by integrating (2.3) over each gyre:

(3.2)\begin{equation} G_W = \int_{{gyre}} W(x,y)\,\textrm{d} x\,\textrm{d} y . \end{equation}

We note here a phenomenon, which was more pronounced in the nonlinear dynamics solutions, and is referred to as the geometric wind effect (Berloff et al. Reference Berloff, Hogg and Dewar2007b). When the eastward jet extension was situated north of the zero wind-curl line, this led to a ${\sim }10\,\%$ attenuation of the wind-curl input for both gyres (table 1). A similar but weak effect occurred in the linear dynamics solutions, as the zero wind-curl line is given by the line of zero meridional velocity instead of zero total velocity (Rhines & Schopp Reference Rhines and Schopp1991). However, this only had a small effect on the wind-curl inputs for each gyre (${<}1\,\%$), so we neglected this effect on the linear dynamics.

The sink terms of the PV budget were computed by considering the viscous boundary and inter-gyre PV fluxes. Viscous boundary fluxes were calculated as

(3.3)\begin{equation} G_{VB} = \nu\overline{\int_{C\backslash\varGamma}\boldsymbol{\nabla}\left(\nabla^2\psi\right) \boldsymbol{\cdot}\boldsymbol{n}\,\mathrm{d}s}, \end{equation}

where $\boldsymbol {n}$ is the corresponding outward normal, $C$ is the contour running along the exterior of the gyre and $\varGamma$ is the inter-gyre boundary. The overline denotes the time averaging over the reference solution.

The inter-gyre viscous fluxes were negligible in the linear dynamics, where they accounted for ${<}0.05\,\%$ of the budget. They also remained extremely small in the nonlinear dynamics solutions, but we included them to close the PV budget. They were computed similarly to (3.3), as

(3.4)\begin{equation} G_{IV} = \nu\overline{\int_{\varGamma}\boldsymbol{\nabla}\left(\nabla^2\psi\right) \boldsymbol{\cdot}\boldsymbol{n}\,\mathrm{d}s} , \end{equation}

except that we integrated over the inter-gyre boundary, rather than along the basin boundary.

The inter-gyre PV fluxes were computed by projecting the PV anomaly fluxes onto unit normal vectors of the inter-gyre boundary:

(3.5)\begin{equation} G_{IF} ={-}\overline{\int_{\varGamma} q\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}\,\mathrm{d}s} . \end{equation}

A proportion of this will consist of inter-gyre exchanges of PV through the action of mesoscale eddies, while the remainder consists of contributions arising from the time-mean flow. Note that PV fluxes owing to the time-mean flow were non-zero because the eastward jet extension did not reach the eastern boundary. The inter-gyre eddy PV fluxes largely presided over the eastward jet, where there was zero time-mean flow across it; while the inter-gyre time-mean PV fluxes occurred further east, where there was a ‘gap’ between the eastward jet and eastern boundary for flow to pass through.

Finally, we must consider the tendency term as, owing to the decadal variability of the ocean gyres in the nonlinear dynamics solutions, its time average is non-zero:

(3.6)\begin{equation} G_T ={-}\overline{\int_{{gyre}}\frac{\partial q}{\partial t}\,\textrm{d} x\,\textrm{d} y} . \end{equation}

Overall, we must have PV anomaly conserved over each gyre, i.e.

(3.7)\begin{equation} G_W+G_{VB}+G_{IV}+G_{IF}+G_T = 0 , \end{equation}

the results of which are shown in table 1. Note that the $\beta$-term contributions to the budget were zero through mass conservation.

Because the WBL lobes of the CGAs consisted of strong negative (positive) advection-induced PV anomalies in the subpolar (subtropical) WBL, we expected a reduction in the zonal relative vorticity gradients near the western boundary, which is confirmed in figure 6. Indeed, this led to ${>}80\,\%$ reduction in viscous boundary fluxes (table 1), which was balanced by the inter-gyre PV fluxes. The sign change in advection-induced PV anomalies at the western boundary did not increase the magnitude of relative vorticity in this region, but instead reduced it, as there was also a corresponding sign change in the relative vorticity near the western boundary. This then led to a reduction in viscous boundary fluxes, as they are governed by normal derivatives of the relative vorticity.

Figure 6.

Time-mean relative vorticity profiles near the western boundary of linear and nonlinear dynamics solutions, and advection-induced anomalies. The generalised Munk boundary layer solution (see Appendix A) is added to compare with the linear dynamics solution. Profiles have been meridionally averaged over (a) subtropical gyre, (b) subpolar gyre.

These results were consistent with those made by Cessi et al. (Reference Cessi, Ierley and Young1987) and Kiss (Reference Kiss2002), where there was also insufficient loss of PV. Cessi et al. (Reference Cessi, Ierley and Young1987) showed that the excess PV then drives the recirculation zones, and Kiss (Reference Kiss2002) argued that this effect may also be involved in the WBC separation. Our results differed from these studies, as we used a double-gyre, rather than a single-gyre model. This gave the model an extra mechanism to remove excess PV, which was through the inter-gyre PV fluxes downstream of the WBCs.

The extent to which the CGAs were consistent with a linear response to the geometric wind effect was checked by reducing the wind-stress amplitude $\tau _0$ by 20 % (We define the linear response to the geometric wind effect as the difference between the time-mean linear dynamics solutions for 80 % and 100 % wind-stress amplitudes.). This was to mimic the impact of the geometric wind effect on wind-curl input reduction for each gyre. We will see from the results of this experiment that even by overestimating the impact of the geometric wind effect, the linear response alone was not strong enough. Relative vorticities, rather than PV anomalies, were used to show the linear response, as they determine the viscous boundary fluxes. Figure 7(a) shows anti-cyclonic (cyclonic) anomalies situated within the subpolar (subtropical) WBL, similar to the WBL lobe of the CGAs seen in figure 5(c). Furthermore, figure 7(b) shows a weaker, second lobe similar to those observed in figure 4(b). However, the linear response to reduced wind-curl input in the ocean interior was basin-wide, rather than being confined to the western half of the gyres. This implied that the shape of the CGAs was consistent with a linear response to the geometric wind effect but this alone could not account for the nonlinear dynamics anomalies, which were significantly stronger. Indeed, the linear response only accounted for ${\sim }15$% of the viscous boundary flux reduction drop, which may also be seen by comparing the PV budgets of the linear and nonlinear dynamics solutions (table 1).

Figure 7.

Contour plots of the time-mean, linear, weakened wind-curl response for relative vorticities in the upper isopycnal layer. Values of Max for each panel are stated for comparison (e.g. with figure 5c). (a) Zoomed-in WBL region, $\textrm {Max} = 8.5\times 10^{-5}\ \textrm {s}^{-1}$. (b) Ocean interior, $\textrm {Max}=1.9\times 10^{-7}\ \textrm {s}^{-1}$.

The separation of the two lobes of the CGAs (figure 5) also suggested that the mechanisms that drive these circulation features were unlikely to be the same. To confirm this, we produced scatter plots of PV anomalies against the velocity streamfunction (see figure 8). In the interior lobe, we observed a functional (nearly linear) relationship between PV anomaly and velocity streamfunction values, which implied the formation of Fofonoff-type gyres (Fofonoff Reference Fofonoff1954). However, in the WBL, the effects of friction made the above-mentioned relationship unfeasible.

Figure 8.

Scatter plots of PV anomalies against velocity streamfunction for grid points in the interior and WBL lobes of the CGAs: (a) subtropical gyre, (b) subpolar gyre.

Examining the inter-gyre PV fluxes in table 1, we see they took up a considerably larger proportion of the PV budget compared with that reported by Berloff et al. (Reference Berloff, Hogg and Dewar2007b). This was attributed to the much smaller eddy-viscosity parameter and a boundary condition parameter choice which created a more turbulent flow regime. Because higher inter-gyre PV fluxes were directly associated with stronger advection-induced anomalies in the WBL, we expected CGAs to only become observable in more nonlinear flow regimes. Along with the drop in viscous boundary fluxes, the inter-gyre PV flux increase accounted for the largest changes in the PV budget for the linear and nonlinear dynamics solutions.

To understand the chain of events that led to both the viscous boundary flux reduction and increase in inter-gyre PV fluxes, we performed a numerical experiment, where we solved system (2.2), but used the time-mean solution to system (2.4) as the initial condition. What we would expect to see is the time-mean linear dynamics solution (figure 3) to converge towards the nonlinear dynamics solutions (figure 2), thus, leading to the adjustments in the PV budgets observed in table 1. The time scales associated with these adjustments can help us to identify the likely causal chains taking place. For example, it may be the case that viscous boundary flux reductions, induced by the nonlinear boundary layer, create a PV accumulation which must be rectified by the inter-gyre PV fluxes. Then, we would expect the time scale of this PV budget adjustment in our experiment to be of WBC advective time scales $T_{{wbc}} \sim L/U\approx 11\ \textrm {days}$, where $L$ is the basin length scale and $U=2\ \textrm {ms}^{-1}$ is the WBC velocity scale. Alternatively, the uptake in inter-gyre PV fluxes could reduce PV concentrations within each gyre, which then reduce the viscous boundary fluxes. If this were the case, then we would expect to see the PV budget adjustment to take place over gyre recirculation time scales $T_{{gyre}}\approx {O}(1\ \textrm {year})$. We found that the viscous boundary flux adjustment took place over an extremely short time frame $({\sim }20\ \textrm {days})$, which indicated that it is the WBC advective time scales that are important in the PV budget adjustment. Hence, it is likely that the nonlinear boundary layer is inducing PV accumulation within the viscous sublayer, which must be rectified downstream by the inter-gyre PV fluxes.

To summarise, the shapes of the CGAs are consistent with the linear, weakened wind-curl response created by the geometric wind effect. However, this cannot account fully for the CGAs, as the linear response by itself is too weak without considering nonlinear effects. In the interior lobe, we see the formation of Fofonoff-type gyres, which arise in nonlinear free-flow situations. Alternatively, in the WBL, nonlinear effects reduce zonal relative vorticity gradients near the western boundary, and this severely inhibits the ability of the WBCs to dissipate PV through viscous boundary fluxes. This creates a growing PV imbalance, which is only rectified downstream by the inter-gyre PV fluxes. We conclude that this rectification of the PV imbalance must be controlled by an inter-gyre PV exchange mechanism. Such inter-gyre mechanisms are not new, e.g. Yang (Reference Yang1996) and Coulliette & Wiggins (Reference Coulliette and Wiggins2001) have studied inter-gyre transports before. However, the link between insufficient PV loss in the viscous sublayers and inter-gyre exchanges of PV has not been made. A similar link between anomalous PV generated in the WBL driving the corresponding recirculation zone has been made before (Cessi et al. Reference Cessi, Ierley and Young1987), but this mechanism acts to strengthen the ocean gyres rather than to weaken them. In the next section, we will switch to a Lagrangian framework, which is more suited to identify the underlying mechanism as well as to confirm the hypotheses made in this section.

4.1. Experiment design

We released $N=400$ particles randomly positioned within the viscous sublayer of the upper isopycnal WBC in 60-day intervals. The viscous sublayer was chosen as this is the region where the PV accumulation has been found to occur (Lozier & Riser Reference Lozier and Riser1989; Kiss Reference Kiss2002). We define the viscous sublayer width to be the perpendicular distance from the western boundary at which the plane parallel to the western boundary gives zero viscous boundary fluxes for each gyre. This coincides with the zonal gradient of relative vorticity profiles changing signs, which we found to be ${\sim }15\ \textrm {km}$ (see figure 6). Particles were advected for a period of 8 years, with a total of 10 releases. This was repeated for the four sets of 10-year model runs computed for the reference flow regime in statistical equilibrium. The particle evolution time length of 8 years was chosen, because it allowed particles on average to make one complete gyre circuit, regardless of where a particle was seeded. Seedings were made up to $7.5\ \textrm {km}$ away from the western boundary, well within the viscous sublayers, which were about two grid cells wide. The total number of particles that permanently migrated from one gyre to the other after the evolution period were compared with the total number released in that gyre. Using the inter-gyre boundary by itself is unsuitable for this task as this produces false identifications when checking in which gyre the particle is positioned. This arises from the high variability of the eastward jet position, which will regularly diverge from its time-mean position. To remedy this, we introduced a buffer zone of $450\ \textrm {km}$ on either side of the inter-gyre boundary, which particles must cross for the particle to be flagged as ‘migrated’. Small adjustments of the width of the buffer zone did not significantly affect any of our results. This experiment was repeated with particles seeded randomly within the ocean interior, which acted as our control case. For Lagrangian particles that migrated over to the opposite gyre, their total PV was approximated using a 2-D cubic spatial interpolation at each time step. This experiment differed from that of Berloff et al. (Reference Berloff, McWilliams and Bracco2002), where a more general analysis of mixing/stirring processes within different regions of the ocean gyres was considered to develop stochastic parametrisations. The analysis in this study was more specific as we are attempting to link insufficient PV dissipation in the viscous sublayers with inter-gyre PV exchanges.

4.3. Inter-gyre PV exchange mechanism

The inter-gyre PV exchange is described by following an ensemble of Lagrangian particles released in the subpolar viscous sublayer (see figure 14 for a schematic).

Figure 14.

Schematic of inter-gyre PV exchange mechanism. Curved black arrow represents the eastward jet extension. Coloured arrows indicate orientation of mass fluxes associated with PV exchange mechanism, each colour and number represents a step of the mechanism described in § 4.3. (1) Fluid parcels in the viscous sublayer lose insufficient PV through the western boundary. (2) The majority of these fluid parcels then enter the recirculation zones where they backscatter. (3) Remaining fluid parcels migrate across the eastward jet extension where they close the PV budget. (4) Migrated particles get acclimatised to the background PV in the mid-ocean.

1. (i) The nonlinear western boundary layer reduces viscous boundary fluxes through effects of PV advection. This leads to an accumulation of PV within the viscous sublayers.

2. (ii) Particles now move downstream where they enter the eastward jet. The majority of these particles then enter the recirculation zone where they fortify the eddy backscatter.

3. (iii) The remainder of these particles migrate permanently to the opposite gyre, where they are able to rectify the PV imbalance created in (i).

4. (iv) Particles enter the Sverdrup-gyre circulation in the eastern half of the basin. By the time they reach the WBCs, they have already acclimatised to the background flow and the process is repeated.