## Abstract

This paper focuses on potential vorticity (PV) budgets in the North Atlantic with an emphasis on the wind-driven subtropical gyre. Since PV is the key dynamical variable of the wind-driven circulation, these budgets are important to understand. PV is a conservative quantity on isopycnals and can only enter or exit through the boundaries, like the lateral topography or the surface. The latter fluxes are diagnosed and tested against the evolution of the PV content in an isopycnal layer. The former are computed using the Bernoulli function. The essential result is found for all the tested isopycnals, and the dominant feature of PV is recirculation, with very little added at the surface or the boundaries. Density coordinates are well suited to understanding PV circulation. A novel technique for computing the Bernoulli function is proposed. The Bernoulli function is governed by a simple elliptic equation and the solutions demonstrate the dominant contribution of PV advection.

## 1. Introduction

In the past 50 years, the quasigeostrophic (QG) approximation has been widely used to describe the ocean [see Berloff et al. (2007) for a recent example]. This approximation reveals the large-scale characteristics of the ocean circulation via the use of potential vorticity (PV). In the QG formalism, sources and sinks of PV are explicitly written into the equation of evolution. For example, the classical double gyre circulation is obtained by applying wind stress curl to the surface layer (Berloff et al. 2007). This formalism is appropriate to describe the large-scale evolution of layers that do not outcrop.

The form of PV appropriate to primitive equations in the Boussinesq approximation is the Ertel PV (e.g., Vallis 2006), written as

with ** ω** =

**f**+

**∇**×

**u**, the sum of the planetary and the relative vorticity,

*ρ*

_{0}a reference density, and

*σ*any thermodynamic variable. We here will choose potential density as

*σ*.

It is now widely recognized that the impermeability theorem applies to PV in the primitive equations (Haynes and McIntyre 1987). This theorem states that PV cannot cross isopycnal surfaces and is exactly conserved in an isopycnal layer (unlike mass). The only locations where PV can enter or exit a layer are at contacts with boundaries like the sea surface or the topography.

The evolution equation for PV is

where **J** is a generalized flux vector (to be defined in section 2). Consider an isopycnal layer. At the air–sea interface, the flux of PV through the surface is given by the vertical component of **J**: **J** ⋅ **k** = *J*_{z}. Similarly, if the layer intersects the topography, the flux of exiting or entering PV is given by **J ⋅ s**, where **s** is the unit vector normal to the topography.

The circulation of PV in the ocean interior is also dictated by the orientation of **J**. This vector is not necessarily parallel to an isopycnal surface; components of **J** normal to the layer imply that the surface will move so that PV does not cross that surface.

Ertel and QG PV differ in that in the primitive equations a range of density layers receive PV at their out-/incrops, whereas in the QG formalism only the upper layer is forced by the atmosphere. This picture is further complicated by the seasonal cycle, which forces a considerable annual migration of the outcrops. This migration allows layers to receive atmospheric forcing contributions from various physical locations.

To better understand the PV dynamics in the ocean, we ask the following questions for a representative subset of isopycnals participating in the wind-driven gyre: How much PV enters the ocean via air–sea interaction? What is the main pattern of PV circulation in the interior and how much PV is created or destroyed at the topography?

To address these questions, we consider the PV budget using the output of a numerical model and separate the PV circulation in the ocean interior from the PV flux at the boundaries. All the numerical computations are performed with a ° run of the North Atlantic with the Nucleus for European Modelling of the Ocean (NEMO; Madec 2008). This run is described in Treguier (2008) and Maze et al. (2013). The atmospheric forcing is the Drakkar forcing set, version 4 (Brodeau et al. 2010), which is essentially based on the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40; Uppala et al. 2005). The run spans 27 years (1980–2006) of which we analyze the last 10 years. The sampling of the output is 5 days.

Such a resolution (°) is necessary to highlight the mesoscale role in the PV dynamics. As shown in Polton and Marshall (2007), the eddy component of the PV dynamics might be of the same order of magnitude as the mean PV dynamics. Our study can be seen as an extension of the work by Marshall et al. (2001), who computed such PV budgets in a noneddying model.

The paper is organized as follows. Section 2 is dedicated to the description of PV flux at the surface. We review several mechanisms that are responsible for generating/extracting PV at the surface. Section 3 is centered on the ocean interior PV circulation. Conclusions are given in section 4. Circulation and surface flux are combined in a companion paper (Deremble and Dewar 2013, hereafter D13) that addresses the dynamics of mode water.

## 2. Surface PV flux

### a. General statements

To understand subtropical gyre PV fluxes, it is necessary to compute the vertical component of the vector **J** at the air–sea interface. The form of the generalized PV flux [see Eq. (2)] can be found by combining the momentum and density equations and is written as

with the Bernoulli function *B*,

being the sum of the pressure and potential and kinetic energy (see Schär 1993). This expression for **J** is strictly identical to

(see also Marshall et al. 2001). In Eq. (5), the vector **J** is the sum of an advective flux (first term on the right-hand side) and nonadvective contributions (the next three terms). The term *ω**Dσ*/*Dt* is the vorticity multiplied by the Eulerian derivative of potential density. The term **F** × **∇***σ* represents mechanical PV forcing (with **F** being the nonconservative forces of the momentum equation) and the last term comes from pressure-dependent effects in the equation of state (sometimes referred to as “thermobaricity”).

### b. Validation in a model

Our first step is to validate Eq. (2) and the impermeability theorem in a numerical model. We compute the time rate of change of the total PV inside an isopycnal layer and compare it to the PV flux at the surface. This first evaluation is done while neglecting the flux of PV at the topography (see also D13 for an analytical derivation). The PV fluxes at the solid boundaries are discussed in section 3.

The integrated PV in a layer *σ* is well approximated by the vertical component of the vorticity multiplied by the vertical stratification integrated over the layer (Vallis 2006),

where the triple integral spans the entire isopycnal layer, and the subscript is a reference to the given isopycnal layer. This simplifies to an integration over an isopycnal layer, and is accurately approximated for a thin layer as

where *σ*_{2} and *σ*_{1} are the upper and lower density limits of the layer. Notice that the thickness of the layer does not appear in the integrated PV (aside from the multiplicative constant *σ*_{2} − *σ*_{1}). It only consists of the sum of the planetary and relative vorticity integrated over the geographic extension of the layer.

Figure 1a is a time series of the integrated PV in the layer 26.0 ± 0.25 kg m^{−3}. Henceforth, we omit the density units when labeling an isopycnal layer. According to Fig. 1a, the integrated PV *Q*^{26} decreases from September to March and increases during the rest of the year. This is consistent with the seasonal evolution of the outcropping. Indeed, from Eq. (7), we know that the integrated PV is a rough measure of the surface extension of the isopycnal layer scaled by the Coriolis parameter. Figure 1a largely reflects the northern migration (in spring) and southern retreat (in autumn) of the outcrop.

Differentiating this time series gives an estimate of the rate of PV production/destruction. This time rate of change of the integrated PV is plotted together with the integral of −*J*_{z} [computed according to Eq. (3)] over the outcropping window (Fig. 1b). This operation consists of defining a spatial mask taking values of 1 if the surface density is between 25.75 and 26.25 and zero elsewhere. This mask is variable in time as it follows the outcropping window. For every output of the model, the field *J*_{z} is multiplied by this mask and the spatial integration of the result is performed. We end up with a time series of 73 points (5-day sample of one year) that is plotted in Fig. 1b.

The correlation coefficient between these two curves is 0.96. A perfect matching between these two curves does not occur partly due to the neglect of processes related to the topography and because of spatial and temporal discretization (see the discussion in the next paragraph). However, the agreement between the two strongly supports the validity of *J*_{z} as expressed in Eq. (3).

In Fig. 1a, we also reconstructed the time series of the integrated PV using −*J*_{z} (blue curve). The red and blue curves are consistent with the curves showing the differentiated time series: *J*_{z} closely reflects the variations of PV (correlation coefficient: 0.99). The curve of integrated *J*_{z} shows an anomalous trend of PV creation, which we discuss next.

### c. Global budget and aliasing issue

We now examine more closely the different components of *J*_{z}. The three contributions to *J*_{z} in Eq. (3) are **k** · **ω**σ_{t}, **k** · (**u**_{t} × **∇***σ*), and **k** · (**∇***B* × **∇***σ*). If we look at a daily snapshot of *J*_{z}, the dominant term is **k** · **ω***σ*_{t} (not shown). This term is in fact another way to look at the huge seasonal meridional migration of isopycnals (several thousands of kilometers in one year). But in a 1-yr period, we expect an isopycnal to return near to its initial position, so the average of this term should be small compared to the other terms. After a long average is taken, we expect that the remnant will be the correlation between the relative vorticity and the time derivative of the density. This statement is true at every location, even when an isopycnal window is used as conditional averaging.

To demonstrate this proposition, suppose that we have perfect knowledge of the density at the surface of the ocean at any instant in time *t*. Suppose also that the state of the ocean at the end of the year is identical to the state of the ocean at the beginning of the year. We want to measure the contribution of the term **k** · **ω***σ*_{t} to the surface PV flux. To a first approximation, this term is equal to *fσ*_{t} (we replace the vorticity by the planetary vorticity). We consider an isopycnal layer with the interval [*σ*_{1}, *σ*_{2}]. At a given location (*x*, *y*) and a given time *t*, if the surface density *σ*(*x*, *y*, *t*) ∉ [*σ*_{1}, *σ*_{2}] then *J*_{z}(*x*, *y*, *t*) = 0. Suppose now the isopycnal layer sweeps by this location (e.g., during the spring migration). As long as *σ* ∈ [*σ*_{1}, *σ*_{2}], there is a component of *J*_{z} that is equal to *fσ*_{t}. The summation (or integral) of all the instants in which *σ* ∈ [*σ*_{1}, *σ*_{2}] returns a global contribution equal to . Then, during the autumn migration, the isopycnal layer travels in the opposite direction and the contribution is now *f*(*σ*_{1} − *σ*_{2}). The sum of these two contributions must be exactly zero.

We have computed the term *fσ*_{t} in various locations and over 10 years and we find that this term is not equal to zero, or even small. In view of the above, we interpret this as an error of the analysis and explain it primarily as being due to the 5-day sampling of the model output. The spring migration of the outcrop is faster than the autumn retreat (not shown). At a given point (*x*, *y*), it is possible that one day in spring *σ*(*x*, *y*, *t*) < *σ*_{2} and five days later *σ*(*x*, *y*, *t*) > *σ*_{1}. In such case there will not be any contribution at this precise location during the spring migration, whereas the real contribution should be *f*(*σ*_{2} − *σ*_{1}) (with perfect sampling). During the return phase, which is slower, there is more chance that at this same location there will be a day for which *σ*_{1} < *σ*(*x*, *y*, *t*) < *σ*_{2}, which means that there will be a contribution on the order of *f*(*σ*_{1} − *σ*_{2}). The global contribution will appear as negative from the analysis when it is in fact zero.

The numerical evaluation of **k** · **ω***σ*_{t} is at leading order equal to *fσ*_{t}. So replacing the planetary vorticity by the total vorticity does not significantly change the result.

We expect the quantity **k** · (**u**_{t} × **∇***σ*) to be an order of magnitude smaller than **k** · (**∇**B × **∇**σ) (Marshall et al. 2001). This proposition is verified later. It remains that the leading-order time-averaged (written with a bar) vertical PV flux following an outcrop is approximately

[the last term of Eq. (3)]. The mean PV flux is given by the intersection of the Bernoulli contours and the potential density contours.

We compute the spatial- and temporal-mean PV flux [using Eq. (3)] on several outcropping windows using Lagrangian averaging. For example, the mean of the blue curve in Fig. 1b (extrapolated for the 10-yr time series) gives the total amount of PV extracted at the air–sea interface in isopycnal 26.0. We normalize the result by the size of the outcrop area: *dσ* = 0.5. The same operation is performed on several density classes from 24.0 to 27.0 and the result is plotted in Fig. 2.

In this figure, the solid line is computed with the full Eq. (3). The dotted line with squares is the component . We confirm that this quantity does not play a major role in modifying the global surface PV flux.

The dashed line with asterisks is . This curve is similar to (not shown), justifying the negligible role played by the relative vorticity in this term. As mentioned earlier, this curve should be very close to zero, but is in fact negative. As explained before, we attribute the negative values of this curve to the 5-day sampling that produces this negative bias.

Finally the dashed curve with crosses is . This curve takes always positive values with local maxima at 24.5 and 26.5. The shape of the curve is similar to the one obtained by Marshall et al. (2001) but our values are twice as large. There are several possible explanations for this. We use an eddy-resolving model (with possible larger values of **∇***B* and **∇***σ*). They use monthly averaged values for their computation that is again acting as a smoothing operator for **∇***B* and **∇***σ*.

Using the curve , we see that all density classes from 24 to 27 experience PV loss at the surface. This includes the mode water density range centered at 26.3, a result that differs from that in Olsina et al. (2013). We ascribe this difference to errors introduced by the mixed layer product they employed. Olsina et al. (2013) recognized this uncertainty as large; our result here supports that concern. Nonetheless, the order of magnitude of the PV flux is small. As a demonstration of this, consider the size of the flux expressed as an equivalent heat flux. Given an outcropping window of about 4000 km × 2000 km and using a mixed layer depth of 100 m, we obtain an equivalent *Q*_{net} = 10 W m^{−2} corresponding to the value of *J*_{z} = 4 kg m^{−1} s^{−2} in Fig. 2. In comparison, the annual-average heat flux in the North Atlantic is *O*(100 W m^{−2}). Nonetheless, if we assert that the ocean is approximately equilibrated in PV forcing (in the 10-yr time series we are using, there are no obvious trends in the subtropical gyre), then the long-term mean . From Eq. (2), the net flux at the surface has to be balanced at another boundary, which we expect intuitively will be centered in the western boundary currents.

### d. Discrepancy with the bulk formulation

When trying to evaluate the surface PV flux from data, it is more appealing to use Eq. (5) instead of Eq. (3). However Eq. (5) remains difficult to evaluate directly, requiring knowledge of many contributions to the momentum and buoyancy equations that are challenging to infer in both observations and models. As a result, observational studies (e.g., Czaja and Hausmann 2009; Deremble and Dewar 2012; Maze et al. 2013) have resorted to scaling analysis to estimate *J*_{z} (see Marshall and Nurser 1992). Specifically, the vertical component of both the buoyancy and the mechanical terms of the **J** vector are approximated by

and

with *Q*_{net} as the air–sea heat flux, *α* and *β* as the thermal and haline expansion coefficients, *h* as the depth of the mixed layer, *c*_{p} as the specific heat capacity, *S* as the surface salinity, (*E* − *P*) as the freshwater flux, ** τ** as the surface wind stress, and

*δ*

_{e}as the depth of the Ekman layer. This depth is usually approximated by

*δ*

_{e}= 0.4/

*f*, with , where 0.4 is an empirical constant, determined from observations (Wang and Huang 2004).

Another contribution to stress-driven PV creation measuring the extraction of PV due to mixing inferred by wind has been recently highlighted by Deremble and Dewar (2012). The scaling proposed for this effect is

This component models the entrainment of water at the base of the mixed layer. This component of the PV flux is always positive since the wind can only mix the upper layer and has been shown to be of the same order of magnitude as . The multiplicative factor (0.7) remains a matter of debate and is subject to adjustment depending on the sign of the buoyancy flux (Deremble and Dewar 2012).

These scaling laws [Eqs. (9)–(11)] are convenient as they are readily evaluated from observations external to the ocean and need little in situ data (only mixed layer depth). On the other hand, their accuracy has not been thoroughly assessed. We use model outputs to address this question, with the conclusion that these scaling laws are appropriate for the fall and winter season. Inaccuracies due to poorly observed shallow summer mixed layers tend to overwhelm the estimates away from the winter, leaving annual averages suspicious.

These scaling laws are compared to the time rate of change of the integrated PV [as was done in the previous subsection with the formulation of *J*_{z} given in Eq. (3)]. In Fig. 1c, the curve is compared to , and in Fig. 1d, is compared with . The correlation between the two curves is 0.79 and 0.84, respectively. In both cases, the flux is usefully estimated in winter whereas it is greatly and persistently overestimated in summer. The addition of helps to reduce the bias in summer but is not sufficient to bring about agreement with .

The primary reason for the difference is that is overestimated in summer because of the highly volatile nature of the thin summertime mixed layer and the nonlinear sensitivity of this quantity to it. We have tried different criteria to compute the mixed layer depth, but all these estimates give roughly the same value of (see also Olsina et al. 2013). It is also possible that heating does not only occur in the mixed layer; rather a fraction of it can be used to warm the layer below the mixed layer. A parameterization of this process would be to use a corrected value of the mixed layer thickness in Eq. (9).

We conclude from this comparison that bulk formulations of PV flux can be useful but must be viewed appropriately. Wintertime PV flux estimates are reliably portrayed in this manner, but summertime processes are too heavily dependent on the variable mixed layer to be trusted. Olsina et al. (2013) obtained annual-average PV flux estimates for the mode waters that suggested on average that the mode waters were receiving PV at the air–sea interface. It is likely that this result was biased due to an overestimated summer contribution.

### e. The spatial structure of surface PV forcing

Figure 3 shows estimates of the PV extraction/input by different processes. These plots are a 10-yr average of a given component of *J*_{z} in the layer *σ* = 26.0 ± 0.25. To compute this average, we compute *J*_{z} for each output of the model (5-day sampling) and we apply a mask taking values of one on the outcropping window and zero elsewhere. Each of these maps is summed and then divided by the number of maps. All the variables needed to compute *J*_{z} either come from the model output or from the forcing used in the model (Brodeau et al. 2010). These forcing terms (wind, *Q*_{net}) are sampled on the same 5-day basis.

To understand these maps of PV extraction/input, we need to make a clear distinction between the daily snapshot of *J*_{z} (not shown) and the time average of *J*_{z} (Fig. 3).

As mentioned earlier, instantaneous PV flux can differ significantly from the time-averaged PV flux. However, it is the daily flux pattern that leads to the idea that PV destruction occurs in winter and PV production occurs in summer. This is not true if we do the entire “summer” average (the period during which the outcrop is north of a certain latitude; roughly from March to September). Much of the summer production is removed in the following winter as the outcrop sweeps south (see also Fig. 1).

The first panel (Fig. 3a) is the component that has been validated in the previous section . In Fig. 3a, the region of the separated Gulf Stream (east of Cape Hatteras) is a PV input region. It corresponds to the wintertime for this outcrop. As we approach the northernmost outcrop (which occurs in summertime), *J*_{z} tends to be positive. This segregation between positive and negative values is seen in all isopycnal layers (not shown). For example, if we consider layer 24.0, which outcrops north of the Gulf Stream in summer (not shown), the blue zone in Fig. 3a is replaced by an orange zone.

We now analyze each individual contribution as given by the scaling laws and compare it to this reference map. Buoyant PV generation is mostly negative, consisting of a net restratifying PV flux. Diabatic heating fluxes PV into the ocean. Indeed, this term scales as *Q*_{net}/*h*, and because the mixed layer depth is much thinner in summer than in winter, the year average of this term will emphasize the summer value (negative).

As mentioned by Maze et al. (2013), mechanical PV generation remains small compared to even with an eddy-resolving dataset. Thomas (2005) suggested that this component can be large in the presence of a sharp density gradient like that in the Gulf Stream. However, persistent alignment of sharp surface density gradients and down-front wind is required to extract sizable PV in the mean. This does not appear to be the case in the Gulf Stream area where is negative (at least in the density range shown here). In very localized places (especially shelf water or the Grand Banks), can be large in the mean, showing values of the same order of magnitude as . This is due to the presence of an intense density gradient (stronger than the one observed in the Gulf Stream) because of the coexistence of shelf and deep water.

Note that can take only positive values (PV extracting). This term counterbalances since these two terms are of the same order of magnitude. Again, this term is maximum in summer as it scales with .

We now compare the sum of these three terms (Fig. 3e) with the full vertical PV flux *J*_{z} (Fig. 3a) computed using Eq. (8). This comparison tends not to support the classically used Eqs. (9)–(11). The estimates are different in many aspects.

The large-scale pattern is different. Whereas for the PV extraction/input occurs in the entire basin, for

*J*_{z}it is concentrated in the western part.The shallow water zones are completely different. This is probably due to the interaction between the mixed layer and the bottom topography leading to a completely different scaling than the one proposed in Eqs. (9)–(11) (see, e.g., the Grand Banks region or the continental shelf where the water is less than 100 m deep).

More fine structures occur in

*J*_{z}than . The filaments present in*J*_{z}seem to be the signature of the eddy field that is barely visible in .The overall magnitude of

*J*_{z}is smaller than . Moreover, the spatial mean of*J*_{z}is positive (PV extraction) and the spatial mean of is negative (PV input).

## 3. Interior circulation

Once PV is injected at the surface, it either circulates in the ocean interior or is re-extracted during the same annual cycle (or perhaps some years later). We now examine the main pathway of PV circulation and the different processes associated with its circulation.

### a. The isopycnal formalism

The impermeability theorem states that PV, defined as in Eq. (1), does not cross surfaces of constant *σ*. This leaves open many choices for *σ*. Historically, the most natural choice is potential density, to draw connections with theory and because of the dynamical importance of density. There are issues that arise with this choice, however, in a model employing a real equation of state. We adopt potential density in this paper to connect with the literature, but point out the consequences of this choice in the following analysis (see also Marshall et al. 2001).

Since an isopycnal surface is impermeable to PV, the latter can only circulate within this layer. The momentum equations in density coordinates are

with *H* = *Dσ*/*Dt* being the diabatic source and Φ = *gz* as the geopotential (Bleck 2002; Griffies 2004). Spatial derivatives are taken along constant *σ* surfaces, but we keep the same notation for the other terms. A subtle point is that at a given location, these equations literally might not exist for a given isopycnal. This occurs at a given point if, for example, the fluid column is all denser than the target isopycnal. If for the moment we restrict ourselves to locations where the isopycnal in question always is found, time averaging the equations removes the time derivatives, resulting in

The left-hand side is recognized as the flux form of PV in density coordinates and the mean Bernoulli function is identified as the streamfunction for the generalized PV flux. It is then a simple matter to show that

with *J*^{y} and *J*^{x} being the meridional and zonal components of the PV flux in density coordinates,

as written on the lhs of Eq. (13). We keep the same notation as in the previous section: even though this applies in a density coordinate system rather than a geopotential coordinate system, the interpretation remains the same.

The relation between the Bernoulli function and the mean circulation of PV has already been demonstrated by Schär (1993) and used by Maze and Marshall (2011). However, the formulation of Eq. (14) provides a useful framework to diagnose and understand the PV circulation.

Solving the elliptic Eq. (14) with the Neumann boundary conditions (∂*B*/∂*x* and ∂*B*/∂*y*) given in Eq. (13) yields the mean Bernoulli function (see the appendix for the method used to solve this equation). Of course, it is also possible to construct the Bernoulli function directly from the output of an OGCM. We perform this computation as well, but argue that using both procedures has an advantage. An issue with the Bernoulli is that it involves two enormous but opposite-signed inputs (the pressure and the quantity *ρgz*). The dynamical content in the Bernoulli is actually much smaller than either of these by a factor of 10 000 (McDougall 2003). Given the inherent errors in a numerical model, it is not obvious a priori that the Bernoulli constituted from the numerical data would be accurate. In contrast, the present approach based on the elliptic equation relates the Bernoulli solely to the dynamical state of the model. Of course, solving an elliptic equation is more difficult than simply creating the Bernoulli, so each technique offers some advantages. Last, we point out that the elliptic operator also allows us to compute the contributions to PV streamfunction from each of the contributing terms: we can explicitly compute the roles of the advection, diabatics, or equation of state [see Eq. (5)].

### b. PV circulation map

To validate this technique for finding the Bernoulli function, we plot in Fig. 4 the mean Bernoulli function for 1997 projected on layer 26.0. The northern boundary of the computational domain is given by the southernmost position of the outcrop during 1997 (red line). The two sets of contours (thick solid and thin dashed lines) show the two different methods used to compute the Bernoulli function. The thick lines are the result of the elliptic equation [Eq. (14)] and the thin lines are obtained with direct calculation [Eq. (4)]. A constant is added to the solution of the elliptic equation to match the Bernoulli contours obtained with the direct calculation (this unknown is due to the Neumann boundary conditions). The two estimates of the Bernoulli function agree to within a few percent, which argues for our purposes that either estimate will suffice.

As for the shape of the Bernoulli contours, we recognize the main features of the subtropical gyre: the Gulf Stream at the western boundary separates at Cape Hatteras and the return flow occurs at all longitudes.

As mentioned earlier, an advantage of solving the elliptic equation is to quantify the impact on the streamfunction of each term of the total PV flux. It appears that this field is almost exclusively due to the advection term on this surface. The effect of diabatic heating [*H* in Eq. (13)] is very small (four orders of magnitude lower). We expect that the ocean interior is largely adiabatic, hence this result agrees with our intuition. The new information here is that the ocean interior, as computed within the NEMO framework, also behaves in a largely adiabatic manner with regard to this important dynamical variable.

We do not explicitly compute the impact of “friction.” Indeed, this effect is expected to be weak. The confirmation of this expectation appears in that the Bernoulli functions computed using our two different approaches are very nearly identical. In principle, friction is an input to the directly constructed Bernoulli function. The absence of friction in our elliptic inversion, coupled with its coincidence with the direct Bernoulli function, implies that friction is a weak effect.

### c. Evolution of the circulation with depth

The projection of the Bernoulli contours on the consecutively deeper 26.0, 26.25, and 27.0 isopycnals appears in Fig. 5 and a surprisingly large impact of thermobaricity appears. For the relatively shallow layer (26.0), the mean circulation of PV is clockwise, in agreement with the global circulation pattern of the subtropical gyre. Using the decomposition of Eq. (13), we obtain that this pattern is mostly the result of the advection term. As we go deeper, the circulation weakens. For the mode water layer (26.25), net PV circulation is quite weak. Thermobaricity, however, is still weaker than advection (not shown). This result is further analyzed in D13. The last panel of Fig. 5 is quite different from the previous two. The circulation of PV is reversed (counterclockwise) and much more intense (the contour interval is the same in all panels). Again, using the decomposition of Eq. (13), this circulation results mostly from the thermobaric effect. In fact, because of the presence of mode water, the depth of this isopycnal surface varies from 700 (below the mode water) to 300 m elsewhere. This variable depth is responsible for a large in situ density gradient on this surface and creates a counterclockwise PV circulation on this density surface.

As mentioned by McDougall (1989), these plots can be counterintuitive because of the thermobaric effect, and the nature of this term depends on the choice of the variable *σ*. If *σ* were chosen as *ρ*, the thermobaric term would vanish, but then the source term *H* in the thermodynamic equation would contain a contribution from the material derivative of the pressure. In a previous study, Marshall et al. (2001) also mentioned this increasing role played by thermobaricity for deep isopycnal surfaces.

### d. Sinks and sources of PV

The Bernoulli function can also be used to diagnose the mean PV flux at the intersection of the density layer with the topography. The isopycnal flux of PV through a closed contour away from the boundaries is given by

with **n** being the unit vector outward pointing normal to the contour and Δ*σ* = *σ*_{2} − *σ*_{1} being the thickness of the layer. This factor is added as a result of the vertical integral of the PV flux over the thickness of an isopycnal layer. Since the mean PV flux is nondivergent in regions of the ocean equatorward of outcrops and topographic incrops [see Eq. (13)], the divergence theorem implies that is zero for any contour that closes within such regions in the fluid.

Consider now an arc from a closed contour bounded by points *P*_{1} and *P*_{2}. The flux through this contour is

We use this to diagnose the topographic PV generation. Consider the contour in Fig. 6, consisting of the topographic boundaries of the 26.0 isopycnal, the equator, and the position of the southernmost outcrop. For our numerical evaluation, we consider a layer of thickness Δ*σ* = 0.5. The net PV flux through this contour is zero since it is a closed contour. The values of the Bernoulli function on that contour are plotted in the second panel of Fig. 6 with the numbers matching the numbers of the first panel.

Starting at point 3 and moving counterclockwise to point 7, an arc that involves topographic interactions with South America and the cross-equatorial PV flux, the Bernoulli function is seen to be almost constant. Equivalently, this is almost a mean PV flux streamline, implying that the PV on this density surface is isolated from the Southern Hemisphere. From point 7 to 12, the Bernoulli function increases considerably. The map indicates that this is a net input of PV due to the atmosphere, as these streamlines all emanate from the southernmost outcrop line. The Bernoulli function increase along this arc is consistent with a southwestward flux of PV, in agreement with the overall sense of general circulation. We could essentially drop a line from point 12 to 4, with the Bernoulli difference between the two of 4 m^{2} s^{−2} being the net westward PV flux of the subtropical gyre, with about half entering the Caribbean Sea and a significant amount of the rest following the Antilles Current to the east of the Bahamas. Moving counterclockwise from point 15 to 4 is a mild Bernoulli increase of about 1 m^{2} s^{−2}. These are locations where the isopycnal grounds into the subsurface Gulf of Mexico, so we here see injection into the gyre of PV from topographic sources. The increase here is consistent with an eastward flux of positive PV.

However, clearly the major action in the figure occurs between points 14 and 15. The line between points 14 and 15 defines the East Coast of North America up to the point of Gulf Stream separation. One sees a growth of the Bernoulli function moving counterclockwise on the curve from Cape Hatteras to the tip of Florida that compares in size to that occurring over the interior of the open Atlantic. These Bernoullis, from about −1.5 to −6 m^{2} s^{−2}, have no counterparts along the surface outcrop and are thus unambiguously identified as of topographic origin. Equivalently, the sharp Bernoulli decrease from point 13 to 14 consists of recycling of interior PV (from Bernoullis 4 to −1.5 m^{2} s^{−2}) and the expulsion of topographically created PV (from Bernoullis −1.5 to −5.5 m^{2} s^{−2}) in roughly equal parts.

What is the fate of the topographically generated PV? The Bernoulli contours show that the topography is connected directly to the outcrop. The northern topographic input, north from Cape Hatteras to Newfoundland, forms its own circuit, directly from the generation site to the outcrop and is seemingly disconnected from the larger subtropical circulation. The U.S. East Coast PV injection takes a fast path into an area that in the mean experiences atmospheric interaction. This appears in Fig. 7 where a detailed view of the regional Bernoulli function appears superimposed on a map of the atmospheric PV injection. North of the outcrop, the PV flux is, strictly speaking, no longer nondivergent. Equivalently, the time-mean flux is the combination of two inputs, one from a streamfunction represented by the mean Bernoulli function and one from a velocity potential determined by atmospheric exchange. It is our experience that the vector sense of the PV flux is well approximated by lines of constant Bernoulli, so the sense of Fig. 7 is that the topographically generated PV enters a region of net extraction dominated by convection. In this sense, the topographically generated PV takes a fast path, entering at subsurface isopycnal incrops, moving north in the Gulf Stream to an atmospheric interaction that draws it out of the fluid. Viewed in terms of PV dynamics, part of the large region of heat loss in the separated Gulf Stream exists to balance the PV flux input to the gyre by the topography.

The regional net loss of PV appears to be stronger than the topographic injection. Overall balance on this isopycnal layer involves interactions of the recirculating fluid with regions of both injection and extraction of PV. Specifically, the recirculating PV near the Gulf Stream enters a region of net PV injection associated with the late winter early spring period. This injection is largely balanced by extraction just to the north. The recirculating fluid continues moving anticyclonically around the region of atmospheric interaction through a broader, but weaker, zone of PV injection. The fluid then exits the region of outcropping, anticyclonically winds around the gyre and shoots back into the atmospherically dominated PV forcing region.

### e. The role of eddies

It results from the comparison of the directly computed Bernoulli function [Eq. (4)] and the Bernoulli function given as the solution of the elliptic equation [Eq. (14)] that the PV flux in the ocean interior is dominated by advection down to about 26.25 (cf. Fig. 4). Deeper than this thermobaricity sets in as a dominant factor, but we for now consider a decomposition of the PV flux in the shallower layers to comment on the role of the eddies. Such decompositions have been used to parameterize the role of the eddies in coarse-resolution models (Greatbatch 1998). The technique normally employed is a thickness-weighted averaging, that is,

With this,

which decomposes the advective flux of PV in Eq. (13) into mean and eddy advection of PV (see also D13). The contribution from each to the Bernoulli function is obtained by solving an elliptic equation of the type of Eq. (14) where the right-hand-side forcing is built from one of the two parts instead of the full **J** vector. The same (reduced) advective component is used as boundary condition for the elliptic equation.

The solution so obtained is not unique because the boundary conditions can always be modified by a function that introduces no curl [see the discussion in Peterson and Greatbatch (2001)]. For that reason, we refer to these components as pseudo-Bernoulli contributions. Even though this ambiguity exists, we remark that the curvature of the field we compute will be unaffected by the boundary conditions, and hence the locations of extrema are unique. We will confine our interpretation of the pseudo-Bernoulli functions focusing on the extrema. We have also verified that the sum of all pseudo-Bernoulli functions obtained by solving the elliptic equation forced by the components of the flux corresponds to the total Bernoulli field with a maximum error of 10%. The full Bernoulli field was almost entirely recreated by the advective pseudo-Bernoulli functions and the thermobaric component.

In Fig. 8a we plot the pseudo-Bernoulli function due to and in Fig. 8b the pseudo-Bernoulli function due to . The main part of the PV circulation is due to the mean flow (note the different contour interval for the second panel). The eddies are mostly active in the Gulf Stream where they act to enhance the circulation set up by the mean flow.

This contribution of the eddies to the circulation of PV makes no statement about their flux divergence. This property is indeed a crucial component of the PV dynamics and is examined for the Eighteen Degree Water in D13.

To get a better understanding of the eddy contribution, we decompose this field into an along-gradient (of mean PV) component and an along-mean contour (or equivalently cross gradient) component

We evaluate *κ* and *ν* according to

Maps of *κ* and and *ν* appear in Fig. 9. These maps are very similar to those presented by Eden et al. (2007) for buoyancy fluxes. The downgradient property is mainly observed in the center of the subtropical gyre with a coefficient *κ* on the order of 3000 m^{2} s^{−1}. In the Caribbean Sea, *κ* is more intense. It is indeed a region of more intense eddy activity (Jouanno et al. 2008). The core of the Gulf Stream is dominated by an upgradient component (negative *κ*). Interestingly, this region coincides closely with the region of topographically created potential vorticity. The global spatial average of *κ* is 1000 m^{2} s^{−1}.

The cross-gradient component is a sizable component of the total field . It is at a maximum in the Gulf Stream core with values exceeding 10 000 m^{2} s^{−1} and corresponds to an enhanced circulation of PV due to the eddies.

This mean/eddy decomposition is used in D13 to explain the maintenance of subtropical mode water. It is shown there that while the mean flux of PV is directed outward of the mode water, it is balanced by an eddy flux of the same magnitude directed in the opposite way.

## 4. Conclusions

In this paper, we clarify the pathways of PV in the North Atlantic subtropical gyre. To address this topic we separate the in-/output flux from the interior circulation.

### a. Surface fluxes

We show that evaluating the PV input via air–sea interaction is not a trivial computation and must be addressed with care. In section 2, we argue that the scaling laws for the vertical PV flux at the surface (Fig. 3) are not accurate. They demonstrate skill in the wintertime, but yield inaccurate estimates in the summertime because of their sensitivity to mixed layer depth (Fig. 1). The complete formulation of flux is also hard to compute even from model outputs due to its dependence on nonconservative effects. The alternate formulation of the flux proposed in Eq. (3) is considerably simpler to use and returns generally more reliable estimates for the PV flux. We do remark, however, that computing the time-dependent flux accurately requires high-frequency observations of sea surface density to avoid aliasing. We compute mean PV budgets for several density classes (Fig. 2), using Eq. (8) to compute the air–sea PV flux. We show that for the density classes of the subtropical wind-driven gyre (from 24.0 to 27.0), the ocean loses PV at the air–sea interface. The annual mean remains very small compared to the daily PV flux.

### b. Ocean interior

The PV circulation in the interior is diagnosed using density coordinates, which reveals the importance of the Bernoulli function in describing the mean flux of PV [see Eqs. (13) and (14)]. This set of equations permits the computation of the Bernoulli function in an indirect way by solving an elliptic equation. We validate this computation by comparing the Bernoulli functions obtained using direct and indirect computation (Fig. 4). This novel technique allows us to directly connect the Bernoulli to its dynamical content and can independently quantify the amplitudes of the several contributions to PV flux.

The evolution of the Bernoulli contours with depth is plotted in Fig. 5. This figure reveals a surprising feature of the mean PV circulation in deep layers. We showed that the thermobaric effect is responsible for overwhelming the clockwise circulation produced by advection and setting up a strong counterclockwise PV circulation (Fig. 5). The mode water layer (26.25) appears to be a layer of very weak circulation. This point is further investigated in D13.

Figure 6 argues that the dominant behavior of the PV is to recycle in the interior. The PV route for this component of the circulation consists of atmospheric input in the central and eastern parts the North Atlantic and atmospheric extraction in the west after a circuit in the gyre. Farther west, against the U.S. East Coast, PV is input to the gyre from topographic sources, and then follows a short path to extraction where convective activity is driven by buoyancy exchange with the atmosphere (Fig. 7). The magnitude of the topographic production is comparable to that of the surface PV flux.

Finally, the role of eddies is investigated in terms of the proposition of Rhines and Young (1982) that they promote downgradient transport. We found mean values for the eddy PV diffusivity of *κ* = 1000 m^{2} s^{−1} (Fig. 9). The primary departure from downgradient flux occurred in the western boundary layer, in the zone affected by topographically generated PV flux. The cross-gradient component of the eddies is far from negligible in the Gulf Stream, although this contribution is not yet fully understood.

### c. Discussion and implication for future work

A full exploitation of Eqs. (13) and (14) requires an accurate method for the Helmholtz decomposition to decompose each contribution into a rotational part and divergent part. This decomposition is not unique and we are now developing a method to find the physically relevant rotational part and divergent part of a vector field (see also Fox-Kemper et al. 2003).

The interaction with the topography and the process by which PV is created is not yet fully understood. We are now trying to address this topic using high-resolution models. However, in contrast to surface PV flux that can be estimated from data, it is not currently possible to compare theoretical (or model) values of the PV flux at the topography with observed PV fluxes.

A direct application of this work can be found in D13. The results obtained here are used to explain the mode water dynamics. We show that the vertical component of the **J** vector does not play a direct role in mode water formation. The creation of low PV water mass is instead more directly related to the surface buoyancy flux.

## Acknowledgments

We thank the Drakkar team and especially A.-M. Treguier, J.-M. Molines, and T. Penduff for providing their dataset and the associated diagnostic tools. T. J. McDougall provided valuable input on the subtleties of the Bernoulli and two reviewers provided perceptive comments. The work reported here was sponsored by NSF Grant OCE-0960500.

### APPENDIX

#### The Elliptic Equation Solver

The resolution of an elliptic equation in an irregular domain is almost as straightforward as it is for a regular (rectangular) domain. We tackle this problem using direct computation of the type used by Gibou et al. (2002) rather than the immersed boundary condition described by Pares-Sierra and Vallis (1989).

Let us consider the Poisson equation

that we want to solve in the discretized domain Ω with Neumann boundary condition. Using a second-order differentiation scheme to solve this equation, one can write, away from the boundary,

with *e*1*u*_{i,j}, *e*1*υ*_{i,j}, …, being the grid steps of the curvilinear C grid at the location (*i*, *j*), as defined in Madec (2008).

Near a boundary, the derivative is simply prescribed: for example, suppose we are at a boundary such that is in the domain and is not. Then, the first term in Eq. (A2) is replaced by the boundary condition

where is the prescribed boundary condition. Using both Eqs. (A2) and (A3) permits us to write Eq. (A1) in matrix form,

with the matrix containing the boundary condition. Since is pentagonal, almost symmetric, it is easy to invert to get the solution **Ψ**. In our case, the field is first interpolated on a regular latitude–longitude grid. Then, we use the direct LU solver of matrix laboratory (MATLAB) to invert .

## REFERENCES

*Fundamentals of Ocean Climate Models.*Princeton University Press, 518 pp.

*Deep-Sea Res. II,*

**91,**128–138.

*Deep-Sea Res. II,*

**91,**84–95,

*Atmospheric and Oceanic Fluid Dynamics.*Cambridge University Press, 745 pp.