# An “objective” definition of potential vorticity. Generalized evolution equation and application to the study of coastal upwelling instability

Yves Morel, Guillaume Morvan, Rachid Bershila, Lionel Renault, Jonathan

Gula, Francis Auclair

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yves Morel, Guillaume Morvan, Rachid Bershila, Lionel Renault, Jonathan Gula, et al.. An “objective” definition of potential vorticity. Generalized evolution equation and application to the study of coastal upwelling instability. Ocean Modelling, 2023, pp.102287. 10.1016/j.ocemod.2023.102287. hal-04274117

HAL Id: hal-04274117

https://hal.science/hal-04274117v1

Submitted on 7 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
1 An “objective” definition of potential vorticity.

2 Generalized evolution equation and application to the
study of coastal upwelling instability.

4 Yves Morel¹,²,*, Guillaume Morvan¹,*, Rachid Bershila¹,*, Lionel

5 Renault¹,*, Jonathan Gula²,*, Francis Auclair³,.*

6 Abstract

In this paper, we propose a form for potential vorticity (PV), rescaled using the Lorenz's rearranged density profile, the novelty being that we here take into account its time evolution. We argue this rescaled PV is more representative of the dynamics, in particular to evaluate the respective impact of mixing and friction on the generation of geostrophic circulation. The impact of mixing at global scale, which only modifies the global stratification at rest, is taken into account in the evolution equation of this “objective” definition of PV, in the sense that it scales the PV changes with respect to its effect on the circulation. Numerically, we show that all terms can be calculated coherently using a single computation cell.

We illustrate our purpose by studying the instability of coastal upwelling currents, using a numerical model at high resolution. The configuration is a periodic flat channel on the f-plane with vertical walls at the southern and northern boundaries. A constant wind is applied over a fluid at rest with an initial linear stratification. An upwelling current forms at the northern coast. After a few days, instabilities develop and vortices eventually emerge with surface intensified cyclones and subsurface anticyclones. We show that these instabilities and eddies are associated with (rescaled) PV anomalies,

<footnote>*Corresponding authors
Preprint submitted to Ocean Modelling
Email addresses: yves.morel@legos.obs-mip.fr (Yves Morel),
guillaume.morvan@legos.obs-mip.fr (Guillaume Morvan),
rachid.bershila@legos.obs-mip.fr (Rachid Bershila), lionel.renault@ird.fr
(Lionel Renault), jonathan.gula@univ-brest.fr (Jonathan Gula),
francis.auclair@aeroobs-mip.fr (Francis Auclair)

¹LEGOS, University of Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France
²LOPS, University of Brest, CNRS, IRD, Ifremer, IUEM, France
³LAero, University of Toulouse, CNRS, UPS, Toulouse, France</footnote>
triggered by mixing and friction.

We describe rescaled PV budgets in a layer bounded between the surface and an isopycnal level. Eulerian and Lagrangian diagnostics allow to analyze irreversible PV production terms, distinguishing the influence of advection, friction (associated with wind stress) and mixing. We find that friction plays the main role, generating negative PV anomalies, while mixing acts to dampen this negative PV production. The association of this negative PV anomaly with the outcropping front leads to the baroclinic destabilization of the upwelling front, creating subsurface anticyclonic vortices and surface intensified cyclonic vortices. Varying the strength of the wind forcing shows that mixing is the most sensitive process, with a net effect that is strongly reduced or even reversed with moderate to weak winds.

When the dynamics is fully turbulent, with filaments and vortices of small sizes, the PV production by mixing and friction is enhanced but the Lagrangian diagnostics are more difficult to analyze, since fluctuations at grid scale become significant and numerical effects – associated with imperfections of the numerical schemes – spoil the PV budget calculation.

Keywords: upwelling; instability; vortices; potential vorticity; mixing; wind stress.

## 1. Introduction

The ocean geostrophic circulation is strongly linked to the Potential Vorticity (PV) field. The most basic representation of the ocean dynamics is the quasigeostrophic (QG) model (Pedlosky, 1987; Cushman-Roisin and Beckers, 2011), based on the conservation of PV – in adiabatic evolution – and the PV inversion principle (the streamfunction, geostrophic velocity or vorticity fields can be inferred from the PV field and boundary conditions). The
QG model is based on a simplification of Ertel's general PV (Ertel, 1942). It
has been successfully used in numerous studies and helped to interpret many
aspects of ocean dynamics, from the ocean circulation at gyre scale (Rhines
and Young, 1982a,b; Luyten et al., 1983; Holland et al., 1984; Rhines, 1986;
Talley, 1988; Marshall and Nurser, 1992), to current instabilities (Charney
and Stern, 1962), geophysical turbulence (McWilliams, 1984), and isolated
vortices (McWilliams and Flierl, 1979; Sutyrin and Flierl, 1994; Morel and
McWilliams, 1997).

More recently, a few studies have analyzed the influence of non-conservative
effects on the evolution of Ertel PV. Haynes and McIntyre (1987) and Haynes
and McIntyre (1990) have shown that there are specific constraints on the
evolution of PV and that the net PV content in a layer bounded by two isopy-
cnic surfaces does not vary. The influence of non-conservative momentum
stress at boundaries on PV evolution has also been analyzed theoretically
(Thomas, 2005; Morel et al., 2006; Taylor and Ferrari, 2010; Benthuysen
and Thomas, 2012, 2013). Since then, several studies have investigated the
consequences of diabatic effects on the ocean dynamics from the prism of PV
modification, from basin scales (see for instance Hallberg and Rhines, 1996,
2000; Czaja and Hausmann, 2009) to meso and submesoscales<fnref n="4" /> (see for in-
stance Morel and McWilliams, 2001; Morel et al., 2006; Morel and Thomas,
2009; Rossi et al., 2010; Meunier et al., 2010; Thomas et al., 2013; Mole-
maker et al., 2015; Gula et al., 2015, 2016, 2019; Vic et al., 2015; Giordani
et al., 2017). However, the link between Ertel PV and the dynamics (vor-
ticity and velocity fields) is not straightforward, which makes the physical

<footnote n="4">⁴Mesoscale refers to horizontal length scales close to the internal radius of deformation (10 to 100 km) and submesoscale to scales below (1 to 10 km).</footnote>
analysis inconvenient.

Morel et al. (2019) (see also Assene et al., 2020; Delpech et al., 2020; Aguedjou et al., 2021; Napolitano et al., 2022; Ernst et al., 2023) proposed a rescaled PV, which is calculated taking into account a reference density profile “representative of the background stratification at rest”. The interest of the rescaled PV is that it is a generalization of the QG PV. It scales as a vorticity with a reference value at rest equal to the Coriolis parameter (Morel et al., 2019; Assene et al., 2020; Napolitano et al., 2022) and deviations of the rescaled PV from its background value at each latitude (also called PV anomalies) are the signature of the vortical geostrophic circulation and can be linked to the dynamics following the QG framework (Morel and McWilliams, 2001; Herbette et al., 2003, 2005; Morel and Thomas, 2009; Le Hénaff et al., 2012). The rescaled PV is conserved for each fluid particle in adiabatic dynamics, but its evolution under diabatic conditions follows similar constraints as the classical Ertel PV (Morel et al., 2019). The generated anomalies are therefore easier to link to the dynamics. Thanks to these properties, the rescaled PV is very useful to analyze the importance of adiabatic and diabatic processes on the generation and evolution of vortices (Assene et al., 2020; Delpech et al., 2020; Aguedjou et al., 2021; Napolitano et al., 2022; Ernst et al., 2023).

However, previous studies using rescaled PV have focused on the sub-surface layers and have assumed that the reference density profile does not evolve. Surface layers are subject to strong diabatic processes, so that the reference profile, which is associated with the background stratification and used to define the rescaled PV, may change with time. In addition, outcropping must be taken into account to evaluate the PV budget of the surface layer (bounded by the ocean surface and a deeper isopycnic level). Indeed,
outcropping of isopycnic levels at the surface is dynamically equivalent to a
positive PV anomaly (Bretherton, 1966; Held et al., 1995; Schneider et al.,
2003; Lapeyre et al., 2006; Lapeyre, 2017; Morel et al., 2019). We here pro-
pose a new form for the rescaled PV, based on Lorenz's rearranged profile
(Lorenz, 1955) and taking into account its time evolution in case of diabatic
mixing. We argue it is the objective form to link PV to the dynamics and
measure the respective influence of friction and mixing on the circulation.
The evolution equation of this new PV formulation is derived as well as the
layer PV budget taking into account surface outcropping. We then illus-
trate the use of this general rescaled PV by studying the development of
barotropic/baroclinic instabilities along coastal upwelling fronts in a simpli-
fied configuration.

Coastal upwellings are of particular interest here, since it is known
that the developing coastal currents, starting from rest, are subject to
barotropic/baroclinic instabilities that generate vortices (Roed and Shi, 1999;
Marchesiello et al., 2003; Capet et al., 2004, 2008a,b). Barotropic/baroclinic
instabilities can only develop if the PV structure has opposite sign gradients
along isopycnic levels (Charney and Stern, 1962; Ripa, 1991). In addition,
(cyclo)geostrophic vortices can only exist if their core consists of a local
PV anomaly. Thus, Morel et al. (2006) argued that, since the initial PV
structure is homogeneous (the ocean is initially at rest), a non-conservative
process must be invoked to generate PV anomalies and explain the instabil-
ity of upwelling currents and the generation of vortices. They also showed
that the friction associated with the wind stress acts differentially along an
isopycnic level, leading to a stress curl even with constant wind, which in this
case produces systematic negative isopycnal PV anomalies (see also Thomas,
2005; Morel and Thomas, 2009). The generated negative PV anomaly then
interacts with the surface outcropping, to produce instabilities and vortices
of both signs. However, Morel et al. (2006) used a layered model and did
not consider diapycnal mixing.

We here resume the study of Morel et al. (2006) in a configuration with
continuous stratification and taking into account mixing. The wind causes
mixing near the surface, so that a mixed layer develops. Since the wind
is constant, mixing occurs throughout the domain and modifies the global
stratification profile. Away from the coast, the evolution is 1D with an
Ekman spiral developing in a deepening mixed layer. The density profile and
PV structure evolve, but no isopycnal PV anomalies are generated, so no
instability vortices can be generated in this region. Near the coast, variations
in stratification and current structures, associated with the development
of the upwelling, locally modify the diapycnal mixing, so that isopycnal
anomalies can be created by mixing or friction. However, in order to assess
the dynamical significance of the generated anomalies, we need to calculate
the rescaled PV with a variable reference profile (representing the global
stratification at rest at any time).

The coastal upwelling configuration is thus a particularly interesting test
case for the rescaled PV evolution we propose here. It is also an interesting
test in terms of physics, and we will also compare the influence of the wind
stress intensity on the PV change and the characteristics of the generated
vortices. In section 2, we present the equations, the numerical model and the
configuration used in the study. The theoretical framework of the rescaled
PV is presented in section 3. We then analyze numerical simulations of
upwelling development and instabilities with the rescaled PV in section 4,
discussing the effect of friction and mixing with strong and moderate winds.
Our results are summarized and discussed in the final section.
## 2. Equations, numerical model and configuration

### 2.1. Equations

In this paper, we consider the Navier-Stokes equations with Boussinesq approximation:

$$ 
\begin{aligned}
& \frac{d}{dt} \vec{U} + \vec{f} \times \vec{U} = -\frac{\vec{\nabla} P}{\rho_0} + \vec{g} \frac{\rho}{\rho_0} + \vec{F} \\
& \vec{\nabla} \cdot \vec{U} = 0 \\
& \frac{d}{dt} \rho = \dot{\rho}
\end{aligned}
\qquad (1)
$$

where $\vec{U} = (u, v, w)$ is the velocity field, $\frac{d}{dt}\phi = \partial_t\phi + \vec{U} \cdot \vec{\nabla}\phi$, $\vec{\nabla} \cdot \vec{V} = \partial_x V_x + \partial_y V_y + \partial_z V_z$ is the divergence of vector field $\vec{V} = (V_x, V_y, V_z)$, $\vec{f} = (0, f_y, f_z)$ is the Coriolis vector ($f$, its vertical component, is the Coriolis parameter), $P$ is the pressure, $\rho$ is the potential density and $\vec{F} = (F_x, F_y, F_z)$ and $\dot{\rho}$ are terms associated with non-conservative processes for momentum (here the wind stress) and density fields (mixing). They are generally prescribed as diffusive terms $\partial_z(K_\phi \partial_z \phi)$, where $\phi$ is the velocity or the density field, and the diffusion coefficient $K_\phi$ is given by some parameterization.

### 2.2. Numerical model

The oceanic simulations were performed with the Coastal and Regional Ocean Community Model (CROCO, see Debreu et al., 2012), developed around the kernel of the Regional Oceanic Modeling System (ROMS, see Shchepetkin and McWilliams, 2005). CROCO is a free-surface, terrain- following coordinate model with split-explicit time stepping, used here with the Boussinesq and hydrostatic approximations. There are several options for the numerical schemes and here we have chosen a third-order, upstream- biased, dissipative advection scheme, with free slip conditions at lateral
boundaries, for horizontal advection of momentum and tracers, (Shchep-
etkin and McWilliams, 1998). There is no explicit lateral viscosity in the
model.

Vertical mixing of momentum and tracers is given by the K-profile pa-
rameterization (Large et al., 1994) with a critical Richardson number of 0.3.
This closure scheme gives vertical diffusion coefficients for momentum and
tracers as a function of wind stress and interior shear (other mixing pro-
cesses, associated with internal waves or double diffusion, can be taken into
account but are not considered here).

2.3. Upwelling configuration

The configuration is a periodic East-West channel with a flat bottom and
vertical walls at the Northern and Southern boundaries (Fig. 1). The chan-
nel depth is $H = 1000 \text{ m}$ and its width and length are similar $L_x = L_y = 250 \text{ km}$.

The stratification is initially uniform with a Brundt-Vaisala frequency $N = 5.2 \ 10^{-3} \text{s}^{-1}$.

The theory, presented in sections 3.1, 3.2 and 3.3, is valid in general con-
figurations taking into account the full Navier-Stokes equations (1), with
variable Coriolis parameter and winds, but to simplify the analysis, in the
numerical simulations presented here, the Coriolis parameter $f$ and the wind
are chosen constant. We fix $f = 7 \ 10^{-5} \text{s}^{-1}$ which gives a first internal radius
of deformation $Rd \simeq 23 \text{ km}$. The grid step is $\Delta x = 1 \text{ km}$. The vertical
sigma grid follows Song and Haidvogel (1994) with stretching parameters
$\theta_s = 5$, $\theta_b = 0$ and $hc = 30 \text{ m}$. We use 80 vertical sigma levels and the
vertical resolution ranges from $\Delta z \simeq 1 \text{ m}$ at the surface to $\Delta z \simeq 50 \text{ m}$ at
the bottom. There are 50 layers in the upper 150 m where the dynamics is
analyzed.

The ocean is initially at rest and we apply a constant West/East wind
196 $W$ (see Fig. 1). The West/East surface stress associated with the wind 197 is calculated using the bulk formulae $\tau_x = \rho_a C_d |W|W$ (and $\tau_y = 0$), 198 where $\rho_a \approx 1.225 \text{ kg.m}^{-3}$ is the air density and $C_d = 0.0012$ the turbulent 199 momentum transfer parameter. Two different wind velocities are considered, 200 strong with $W = 8 \text{ m/s}$, and moderate with $W = 2 \text{ m/s}$, corresponding 201 respectively to $\tau = 0.1 \text{ N/m}^2$ and $\tau = 0.006 \text{ N/m}^2$. Buoyancy (heat 202 and freshwater) fluxes are not considered here. A quadratic bottom drag, 203 with a drag coefficient $\text{rdrg2} = 2.10^{-3}$, is used, but it plays no role in 204 the transformation of PV in the upper layer. We focus on the Northern side 205 where the upwelling takes place. The downwelling occurring at the Southern 206 boundary is not studied (PV is also modified in the downwelling area, but 207 does not lead to instabilities. See Morel et al., 2006, for more details). For 208 each simulation, we record hourly outputs of all physical fields (velocity 209 $\vec{U} = (u, v, w)$ and density $\rho$) but also of the parameterized non-conservative 210 terms $(\vec{F} = (F_x, F_y, F_z)$ and $\dot{\rho}$).

211 **3. Potential Vorticity**

212 If the domain was infinite (without Southern and Northern boundaries), 213 the evolution would be 1D. Indeed, in this case an Ekman spiral devel- 214 ops and the upper part of the water column undergoes mixing (see Fig. 2 215 a). The stratification and Ertel PV are modified, with a homogenized up- 216 per layer overlying a pycnocline with stronger stratification than the initial 217 stratification. However, there are no horizontal variations of this structure, 218 and no geostrophic motion is associated with this homogeneous stratifica- 219 tion change, so the Ertel PV change is not significant in terms of geostrophic 220 dynamics.
<figure><img src="image_3.png" /><figcaption>Figure 1: Configuration characteristics, horizontal view (panel a) and vertical section (panel b).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Vertical section of the stratification evolution (panel b). In the open ocean the modifications are close to 1D (panel a), with an Ekman spiral developing in the mixed layer. A pycnocline forms but the stratification below is unchanged. Close to the vertical boundary, where the upwelling develops (panel c), the stratification and current structures are very different. Mixing is modified.</figcaption></figure>
In the current configuration, far from the boundaries, the evolution is close to 1D as depicted above. However, the stratification varies close to the boundary, where the upwelling develops. Isopycnic levels bend towards the surface and a vertically sheared geostrophic current develops. Both processes locally modify the mixing characteristics (Fig. 2 c). As discussed in Morel et al. (2006), even though the wind stress is constant, the momentum stress varies vertically, so that a stress curl is created along isopycnic levels bending towards the surface<fnref n="5" /> (see also section 4.1.2 and Fig. 10).

The stress curl and differential mixing create PV anomalies in the upwelling region and an associated dynamical signal, in particular leading to destabilization of the upwelling current (Morel et al., 2006). It is therefore interesting to isolate the PV anomalies generated by diabatic processes that are truly associated with geostrophic dynamics, and to measure the respective influence of mixing and friction on the observed dynamics. But such a measure depends on the definition of PV and the latter has to be carefully defined to get objective diagnostics. This is a delicate issue that we now discuss.

## 3.1. An objective definition of Potential Vorticity

As shown by Ertel (1942) (see also Muller, 2006), vorticity is not conserved, but a quantity combining vorticity and stratification can be defined that is conserved for each particle for adiabatic motions:

$$PV_{Ertel} = (\nabla \times \vec{U} + \vec{f}) \cdot \nabla \rho \quad (2)$$

<footnote n="5">⁵Note that the stress parameterization is generally sensitive to the vertical structure of mixing and thus also varies horizontally in the present case, even though the wind stress at the surface is constant.</footnote>
All dynamical fields (vorticity, velocity and stratification) can be calculated from PV under the simple assumption of (cyclo)geostrophic equilibrium and given boundary conditions. Most studies, using Navier-Stokes or Primitive Equations, invoking PV are based on Ertel PV, including recent studies analyzing the influence of mixing and/or friction. In fact, as already noticed by Ertel, PV is not uniquely defined in the sense that any form of the type

$$PV_{Gene} = (\vec{\nabla} \times \vec{U} + \vec{f}) \cdot \vec{\nabla} G(\rho) \quad (3)$$

-where G can be an arbitrary function- is still a Lagrangian tracer for adiabatic motions. So one may wonder which choice for G is the most objective to define PV. This is particularly important since in case mixing and friction are considered, the evolution equation for $PV_{Gene}$ is

$$\frac{d}{dt}PV_{Gene} = \vec{\nabla} \cdot [ (\vec{\nabla} \times \vec{F}) G(\rho) + (\vec{\nabla} \times \vec{U} + \vec{f}) \partial_{\rho}G(\rho) \dot{\rho} ] \quad (4)$$

This equation clearly shows that the strength of the PV fluxes associated with friction (first term on the right hand side) or mixing (second term) strongly depend on the choice of G. The form of Ertel PV (Eq. 2) with $G(\rho) = \rho$ has no objective justification apart from its simplicity. In fact, the interest of PV being its link with the circulation, the adequate form has to be quantitatively representative of the dynamics, that is to say directly invertible to estimate the vorticity field. This is not the case of Ertel PV.

The proper choice for G can be defined considering a purely barotropic circulation. Indeed, in this case the dynamics is independent of the vertical position z (oriented upward), and the background stratification is horizontally homogeneous: $\rho(x, y, z) = \rho^*(z)$ and is also the stratification at rest. To be representative of the dynamics, PV should correspond to the
barotropic vorticity in this case and be independent of $z$. The only choice is then that $G(\rho) = Z(\rho)$ where $Z$ is defined using the profile $\rho^*(z)$ so that $Z(\rho^*(z)) = z$ (the function $Z(\rho)$ corresponds to the vertical position of density $\rho$ along the reference profile $\rho^*$). In this case $PV = \partial_x v - \partial_y u + f$ is the expected absolute vertical vorticity with a value at rest corresponding to the Coriolis parameter. As argued in previous studies (Delpech et al., 2020; Asene et al., 2020; Aguedjou et al., 2021; Napolitano et al., 2022; Ernst et al., 2023), in the general case, $\rho^*$ should correspond to the stratification at rest, obtained from the Lorenz's (Lorenz, 1955) rearranged stratification (see also Nakamura, 1995; Winters and D'Asaro, 1996). Indeed, in this case only, the rescaled PV at rest is still given by the Coriolis parameter $PV_{rescaled}^{rest} = f$ and, at first order, the anomaly from this reference corresponds to the quasi-geostrophic PV, which is again directly linked to the geostrophic circulation. The link between the geostrophic circulation and rescaled PV is rigorous if the reference density profile chosen for rescaling correspond to the stratification at rest at each time, indicating the necessity to take into account time evolution of the reference profile. Indeed, Lorenz's rearranged profile can be modified by large scale diabatic processes, so that $\rho^* = \rho^*(z, t)$ also depends on time. For instance, considering again a barotropic circulation for which large scale ocean/atmosphere fluxes uniformly modify the stratification but not the dynamics (vertical mixing of momentum does not modify barotropic circulation), PV should remain equal to the absolute vorticity whatever the evolution of the background stratification. In addition, in general circumstances, mixing and ocean/atmosphere fluxes can destroy or create density classes in a domain, which can cause problems to evaluate rescaled PV if no time evolution is considered for the reference profile.

To achieve these properties and be representative of the dynamics at all
times, the objective form for the definition of PV is

$$PV_{\text{rescaled}} = \nabla \cdot ( (\nabla \times \vec{U} + \vec{f}) Z(\rho, t) ) \quad (5)$$

where $Z(\rho, t)$ is a function of both potential density and time and is defined using the time evolving rearranged Lorenz's profile $\rho^*(z, t)$: $Z$ is chosen so that at all times $Z(\rho^*(z, t), t) = z$. This is what is chosen in the present study to compare effects of mixing and friction on the dynamics, through their associated PV fluxes terms. In our numerical simulations, the calculation of the reference profile is based on previous work (Tseng and Ferziger, 2001) with an adaptation associated with the present configuration (see Appendix A).

Interestingly, previous studies, to our knowledge only using non-evolving Lorenz's profiles, have identified other interesting properties of the rescaled PV. The PV anomaly proposed by Morel and McWilliams (2001) (see also Herbette et al., 2003, 2004, 2005; Morel et al., 2006; Morel and Thomas, 2009) is similar to the rescaled PV for adiabatic evolution, but written in isopycnic coordinate (or for multi-layer shallow water models). This quantity was already shown to be a generalization of the QGPV for primitive equations (written in isopycnic coordinates) and was used instead of the classical form because it allows a direct inversion of PV to calculate the (cyclo)geostrophic circulation. The available potential vorticity, proposed by Wagner and Young (2015) (see also Early et al., 2021) for the filtering of internal gravity waves, is similar to the rescaled PV. Gravity waves have no PV signature along isopycnic surfaces, but in a Eulerian framework when using the classical form of Ertel PV, the signature of the pycnocline induces PV variations that are difficult to attribute to geostrophic dynamics or the displacement of the pycnocline by gravity waves. The use of a rescaled PV
Hereafter, PV will refer to the rescaled PV defined by Eq. 5 and whose evolution equation is Eq. 6 or 7.

For the upwelling simulation presented above, using the rescaled PV given by Eq. 5 ensures that the PV at rest is unchanged ($PV_{rescaled}^{rest} = f$) and that the isopycnic PV anomalies are representative of the geostrophic dynamics and in particular instabilities.

Finally, other forms are sometimes used for the right hand side of Eq. 6, but the present expression has strong similarities with the PV expression in Eq. 5, and is the preferred form for numerical calculations because all terms can be calculated using a single PV grid cell, which also simplifies coherency of PV budgets (see Morel et al., 2019, and Appendix C).

### 3.3. Average PV field in a layer

Since (rescaled) PV represents quasigeostrophic PV at first order, a vertical average of PV in a layer bounded by two isopycnals is informative about dynamical fields. For example, an anticyclonic vortex is associated with a negative PV anomaly, localized in a core within some layer, and the PV budget within the layer is related to the vorticity and stretching fields (see for instance McWilliams and Flierl, 1979; Hoskins et al., 1985; Morel and McWilliams, 1997; Le Hénaff et al., 2012). When the upper bound of the layer is the sea surface, which is generally not an isopycnic surface, an additional term associated with density variations at the surface -outcropping- has to be taken into account when evaluating the PV budget in relation with the dynamics. Indeed, density variation along the surface is equivalent to a Dirac delta sheet of PV that has to be taken into account (Bretherton, 1966; Schneider et al., 2003; Isern-Fontanet et al., 2006; Lapeyre et al., 2006). In this case, the correct calculation for the integrated PV is (Schneider et al.,
2003; Morel et al., 2019; Ernst et al., 2023)

$$ \overline{PV} = \frac{1}{h} \left( \int_{z=-h(\rho_{low})}^{z=0} PV \, dz - \left[ (\vec{\nabla} \times \vec{U} + \vec{f}) \cdot Z(\rho, t) \right]_{z=0} \right) \quad (8) $$

where $h = h(\rho_{low})$ is the depth of the isopycnic level that defines the lower boundary of the layer and $Z(\rho, t)$ is the vertical level at density $\rho$ of the reference profile $\rho^*(z, t)$. The additional term is calculated from velocity and density fields at the surface ($z = 0$). Note that for an upwelling, the region where deep isopycnic levels outcrop is equivalent to a positive PV anomaly ($Z(\rho, t)$ corresponds to the depth of density $\rho$ along the reference profile and is always negative), which can potentially generate cyclonic eddies (Bretherton, 1966; Legg and Marshall, 1993; Legg et al., 1996)$^6$.

## 4. Numerical results

### 4.1. Reference experiment

For the reference experiment, the wind stress is $\tau = 0.1 \, N/m^2$ and the model is run for 25 days with hourly outputs, starting from rest. Figure 3 shows the evolution of the integrated PV anomaly (using Eq. 8 and a layer bounded by the isopycnic level $\rho = 1025 \, kg/m^3$ and the surface), surface vorticity, a vertical section of PV anomaly and vorticity along a South-North transect. The PV anomaly is obtained from Eq. 5 and withdrawing the reference PV ($f$). The vertical sections are taken along the dashed line shown in the horizontal maps. It is a South-North section located at $x = 205 \, km$. Note that it crosses the core of a subsurface anticyclonic vortex at time $t = 21 \, days$ (Fig. 4), whose generation will be studied below.

<footnote>$^6$Physically, if there is no PV anomaly inside the water column but isopycnals outcrop at the surface, internal layers must be stretched and cyclonic vorticity is then created.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Evolution (days 1, 4, 16. Each day corresponds to a row) of horizontal maps of surface layer averaged PV anomaly $\overline{PV}$ (first column), surface vorticity (second column) and vertical sections of PV anomaly (third column) and vorticity (fourth column). In rows 1, 3 and 4, the solid black contours are associated with isopycnic levels.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Continued (days 19, 21 and 25).</figcaption></figure>
During the early stage, up to day 15, the evolution remains zonal, no instability is visible. From day 15 to 17, short scale instabilities with a maximum length scale of 10 km develop, but remain trapped near the coast. Then, from day 18 to 25, new instabilities develop with a dominant length scale starting at 20 km but evolving to about 40 km at day 25⁷. The evolution is similar to that found in Morel et al. (2006), with initial, short scale instabilities developing when the Richardson number becomes smaller than 0.25, corresponding to Kelvin-Helmholtz like instabilities<fnref n="8" /> and larger scale eddies developing from day 18 corresponding to geostrophic barotropic/baroclinic instabilities. At later stages, the merging of eddies also favors the formation of larger scale vortices.

Initially, at days 1 and 4, the evolution of the vertical section of PV anomaly and vorticity (columns 3 and 4) is limited to a region close to the northern wall: although mixing modifies the stratification and Ertel PV everywhere in the fluid, the rescaled PV mostly exhibits anomalies where it is dynamically significant, i.e. the generated PV can generate eddies. Indeed, at a later stage, when eddies develop (days 19, 21 and 25) there is a striking correspondence between the surface vorticity and the layer averaged PV. In particular, anticyclonic eddies are associated with negative PV anomalies and cyclonic eddies with positive ones. Note that the positive layer PV anomalies are associated with outcropping of isopycnic levels (red areas associated with positive PV anomalies and black contours, associated

<footnote n="7">⁷All length scales have been evaluated as the dominant wavelength of the x-axis Fourier transform of the surface PV in the upwelling area.</footnote>

<footnote n="8">⁸The initial growth of perturbations in KH instabilities is correctly reproduced in primitive equations, but their subsequent evolutions require non-hydrostatic dynamics.</footnote>
with surface density variations, in the first column) and that the order of magnitude of vorticity and PV are similar. Thus, as expected, the (rescaled) PV anomaly can be integrated to evaluate the geostrophic dynamics.

As discussed above, outcropping is always associated with equivalent positive PV anomaly (Bretherton, 1966; Schneider et al., 2003). If cyclonic vorticity seems to be mostly associated with the outcropping of isopycnic surfaces, both positive and negative interior PV anomalies are generated, as shown in the vertical section of PV anomaly (third column). They are necessarily associated with non-conservative processes. Thus, we now analyze the importance of diabatic mixing and friction in modifying the interior PV.

### 4.1.1. *Eulerian analysis*

The simplest diagnostics of PV production by mixing and friction are Eulerian maps of the terms in Eq. 7. As mentioned above, PV, friction and mixing tendency terms can be calculated coherently over a single numerical cell following Appendix C, but adiabatic Eulerian advection must be considered too.

Figure 5 shows snapshots of all terms appearing in the PV evolution equation (Eq. 7). We focus on the initial stages of the PV evolution when the dynamics remains 2D. The Eulerian PV time derivative at a given time $t$ is obtained from the difference between the PV fields at time $t + 1h$ and $t - 1h$ (using our hourly outputs). Other terms correspond to the tendency terms on the right hand side of Eq. 7: Advection, Effect of friction and diapycnal mixing. The good coherence between the sum of all terms (second column) and the Eulerian time derivative of PV (first column) confirms the accuracy of our computation. Note that for the present configuration high frequency outputs are necessary to obtain this accuracy. Indeed, using lower
frequency outputs (daily or even 12-hourly) can lead to significant discrepancies, especially when the geostrophic instability develops, since smaller scale structures appear and are rapidly advected, leading to rapid changes in dynamical fields. Morel et al. (2006) have argued that the effect of the wind stress associated with the friction tendency term (fourth row of Fig. 5) should produce negative PV anomalies (see also Thomas, 2005), as long as the stratification remains 2D and with constant wind. From day 15 on, the stratification becomes 3D, positive anomalies appear in the friction tendency term (see day 16, fourth row in Fig. 5), possibly generated by the wind stress acting on 3D structures (see for instance Morel and Thomas, 2009, for the influence of the wind stress on vortices). Note that the lateral friction could also play a role (even free slip boundary conditions can lead to the generation of PV anomalies, see Morel and McWilliams, 2001) but since horizontal viscosity is implicit, it is not calculated in the friction tendency term and cannot be linked to the positive anomalies seen in Fig. 5.

From the impermeability theorem (Haynes and McIntyre, 1987, 1990), one would expect mixing to produce negative PV anomalies in the upper part of the mixed layer – where the fluid is homogenized – and positive ones just below – where a pycnocline appears – (see also Morel and McWilliams, 2001). However, with the rescaled PV definition depending on time, this principle is modified. Indeed, in Eq. 7 the time evolution of the reference density (which represents the effect of mixing on a global scale) is discarded from the local diabatic density change, since it does not affect the geostrophic dynamics. Only the diapycnal flux anomaly with respect to the reference density evolution has a dynamical effect and is considered here.

This explains the particularly striking fact that over most of the mixed layer, the mixing tendency term remains zero away from the boundary (see
<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution (days 1, 4, 16. Each day corresponds to a row) of vertical sections, taken along a South-North section located at $x = 205$ km, of $\partial_t$PV (first column), the sum of all right hand side flux terms appearing in the Eulerian PV evolution Eq. 7 (second column), the advective PV flux (third column), the frictional term (fourth column) and the diapycnal mixing term (fifth column).</figcaption></figure>
Fig. 5 fifth column), even though mixing continuously modifies the density profile throughout the basin, with a pycnocline that penetrates deeper and deeper throughout the simulation. As discussed above, it is in fact the modification of stratification and velocity in the upwelling area near the coast that generates differential mixing and dynamically significant mixing tendency terms. Since isopycnic levels bend upward, as long as the mixing homogenizes a roughly similar depth portion of the water column, we can expect a similar behaviour as described above for the impermeability theorem (the isopycnic levels being mixed in the upwelling region are initially too deep to be affected by the mixing offshore, so the correction term for these isopycnic levels is zero). This is indeed what is seen close to the boundary in the upwelling area (fifth column in Fig. 5). However, the mixing associated with rapid offshore advection of denser water at the surface (generating convective mixing) leads to a strong deepening of the mixed layer in the upwelling area and the development of patterns that alter this simple rule, with a complex structure of the mixing tendency term.

In fact, it is difficult to evaluate the importance of each process and tendency term on the global evolution of PV and the generation of eddies from the Eulerian perspective. In particular, the advective term – which is not associated with PV generation – is of paramount importance. In addition, the tendency terms for mixing and friction often have similar magnitude and extension but with opposite signs, so that the net effect is difficult to evaluate. Note in particular the partial compensation between friction and mixing terms at days 16, when submesoscale structures appear, a process that has been observed and explained in Wenegrat et al. (2018). A better approach is to consider Lagrangian diagnostics, and this is what is presented in the next section.
Finally, Fig. 6 compares the previous results (upper panels) to Ertel PV and the associated traditional frictional and mixing productions terms at *t* = 1 *day* (lower panels). To obtain a scaling for Ertel PV that can be compared to the rescaled PV, we have here used the rescaled form but with the fixed reference profile given by the initial density profile. Indeed, since the initial state has a linear stratification, the non-evolving rescaled PV and tendency terms are proportional to Ertel PV calculations, and quantitatively comparable to the (time evolving) rescaled one. This also allows to illustrate the necessity to use a time evolving reference density profile.

As expected, the Ertel PV structure mostly represents the stratification.
In comparison with the rescaled PV, the mixing production term is primar-
ily marked by the homogenization in the mixed layer and the deepening
of the pycnocline over the whole basin, which is not representative of the
geostrophic circulation developing in the upwelling area. On the other hand,
at this stage, the frictional term for both formulations are pretty close (see
Fig. 5). Thus, measuring the respective influence of mixing and friction on
the development of the geostrophic circulation using this traditional form of
PV can be biased. In the present case negative PV production by mixing
production would be over evaluated.

4.1.2. *Lagrangian Analysis*

A Lagrangian approach is easier to interpret, since we can integrate
the global effect for particles forming the core of eddies emerging from the
instability of the upwelling current. However, it is numerically more de-
manding since it is combined with a tracking algorithm that has some posi-
tioning uncertainties associated with spatial and temporal resolution (again,
the hourly output is a minimum to achieve good positioning in the rapidly
<figure><img src="image_8.png" /><figcaption>Figure 6: Upper panels : rescaled PV (left panel), frictional term (middle panel) and diapycnal mixing term (right panel) at $t = 1 \ \text{day}$ (similar as in Fig. 3 and 5). Lower panels : Ertel PV (left panel), frictional term (middle panel) and diapycnal mixing term (right panel) at $t = 1 \ \text{day}$. Colorscales are given for PV and tendecy terms.</figcaption></figure>
evolving dynamical structure of the upwelling). The coherence between the
PV evolution of a particle and its evaluation from the tendency terms (here
only friction and mixing) following the particle in Eq. 6 is thus even more
difficult to achieve.

We first identify particles with negative PV anomalies in an anticyclonic
vortex whose center is located near x = 205 km, y = 235 km at day 21 (Fig.
4). 2000 particles were randomly seeded near this position, within a radius
of 10 km and with PV anomalies δPV < −10<sup>-5</sup>s<sup>-1</sup>. We backtracked them,
using time interpolated hourly velocity fields, to obtain their trajectories.
The tracking algorithm is the same as that used in Assene et al. (2020).
The position of the particles as well as the values of PV and the friction and
mixing tendency terms are calculated every 600 s, interpolating the hourly
outputs. For each particle, we then integrate temporally the mixing and
friction tendency terms obtained from the record of the non-conservative
terms (right hand side of Eq. 6) to reconstruct a PV evolution from these
diagnostics and we compare it to the PV variations of the particle directly
diagnosed from the physical fields. The correspondence is generally good,
but given the complexity of the tendency terms, the uncertainties associated
with the calculation of trajectories and the hourly sampling, there exist
some discrepancies. We then only keep the 1252 particles (representing 63%
of the initial 2000 particles) for which the correlation between both PV
evolution estimates is greater than 0.9. Figure 7 shows their initial (day 1,
in red) and final (day 21, in blue) positions. All particles forming the core
of the anticyclonic eddy were initially deeper than their final depth. They
also originate from a coastal band that extends about 350 km alongshore
(considering the periodicity of the domain) but only 25 km cross-shore, close
to the internal deformation radius. Figure 8 shows the PV evolution rate,
<figure><img src="image_9.png" /><figcaption>Figure 7: Initial ($t = 1$ day, in red) and final ($t = 21$ days, in blue) positions of selected particles constituting the anticyclonic eddy core. Panel (a) is an X/Y view, panel (b) is a side X/Z view and (c) a side Y/Z view.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

Figure 8: PV evolution rate for a single particle (with maximum correlation between direct PV diagnostic and reconstructed from tendency terms). The blue dashed line corresponds to the calculation from the model PV (calculated using an interpolation of the model PV on the particle’s trajectories) and the black plain line is the diagnostic associated with the sum of friction and mixing tendency terms (again interpolated on the particle’s trajectories).

every 600 s, calculated directly from physical fields and diagnosed from the sum of the friction and tendency terms for a particle with a typical evolution: the dashed blue curve is obtained from the PV record (calculating $\frac{d}{dt}PV$) and the plain black curve is the right hand side of Eq. 6, calculated from the non-conservative terms, along the particle’s trajectories. Note the good agreement between the two estimates before the circulation becomes 3D and exhibits small scale variability. This is generally the case for the retained particles, however, some discrepancies can be seen on the evolution rates, in particular, the model PV evolution rate exhibits abrupt changes that we do not see in the diagnosed PV (using tendency terms). But the variations are close enough to be exploited.

In Fig. 9 represents the latitudinal (upper panel) and vertical (middle
<figure><img src="image_11.png" /><figcaption>Figure 9: Latitudinal (upper panel) and vertical (middle panel) positions and PV evolution rate (lower panel) for the same particle as in Fig. 8. For the PV evolution rate, the green curve represents the friction tendency term and the red curve the mixing tendency term, the black plain line is the sum of both and is the same as in Fig. 8.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>b) Effect of wind stress along isopycnal</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>c) Upwelling</figcaption></figure>

<figure><img src="image_14.png" /></figure>

Figure 10: a) Mixing produces negative PV anomalies in region where the stratification is reduced and positive ones in regions where stratification becomes stronger. b) The effect of friction on PV is associated with the curl of the stress, but calculated along an isopycnic surface. The wind stress generally diminishes in amplitude from the surface to the base of the mixed layer. Thus even an homogeneous wind can be associates with a curl, and PV production, along an isopycnal whose vertical position varies. c) In the upwelling case, mixing deepens the mixed layer and thus increases PV at its base and decreases it in the homogenized region. Isopycnic surfaces bend upward and are thus subject to an increasing along-wind stress generating negative curl and negative PV anomalies. In this case, a Lagrangian particle first moves along an isopycnic surface when it is below the mixed layer, and its properties are unchanged. Eventually, it enters the mixed layer where it is subject to friction, which modifies its PV and produces negative anomalies, and mixing, which changes both its density and PV. Mixing generates first positive and then negative anomalies as the particle moves upward.
panel) position of the particle and the lower panel shows the respective contributions from friction (green) and mixing (red) to the global tendency term (black plain line, already plotted in Fig. 8). The behaviour observed here is the same for most particles: the particle first moves upward and Northward – following an isopycnic surface – but is initially not subject to diabatic forcings and the PV remains unchanged. PV is modified when the particle enters the base of the mixed layer. As explained in Fig. 10, it is first affected by mixing, producing positive anomalies (Fig. 10a and c). Soon after, friction acts too with a negative production rate (Fig. 10b and c). As the particle moves further upward into the mixed layer, it is entrained by the surface Southward currents and the mixing PV production strongly decreases (it can even become negative for some particle) but the negative production by friction is maintained. The global (friction+mixing) PV evolution thus evolves from positive to negative with longer and stronger negative production, so that the final PV anomaly of the particle is generally negative. After this strong modification episode when entering the mixed layer, the particle is subject to a phase where mixing and friction production equilibrate so that its PV value is maintained. After day 18, when the instability develops and produces small scale variability, marked by higher frequency changes in horizontal and vertical positions, PV evolves again but the physics of this phase is more difficult to explain since there are strong numerical uncertainties.

Figure 11 shows the mean evolution calculated from the 1252 selected particles with high correlation. The mean PV anomaly gradually decreases (plain black and dashed blue curves) and the final value is strongly negative (the PV anomaly is $\delta\overline{PV} \sim -6.10^{-5}s^{-1}$). Note that after day 18, discrepancies appear between the model PV (blue dashed curve) and the
542 PV diagnosed from the tendency terms (black plain curve). As explained
543 above, the development of the instability generates small, rapidly evolving
544 eddies. Implicit diabatic terms, associated with imperfections of the numer-
545 ical schemes at the grid scale and the limited (hourly) output frequency used
546 for our diagnostics, corrupt the PV budget. However, the development of the
547 PV structure responsible for the initial instability is correctly represented
548 before the emergence of small scale 3D structures, During this period, fric-
549 tion produces negative PV anomalies (green curves in Fig. 11). Although
550 it produces both positive and negative anomalies, the effect of mixing is
551 to dampen the effect of friction by producing a net positive anomaly (red
552 curves) at a rate that is – on average – similar in amplitude in this refer-
553 ence experiment. For particles entering the lateral boundary layer, lateral
554 friction can also affect the friction tendency term with an effect that can
555 be more difficult to predict as it depends on details of the velocity varia-
556 tions near the wall and the horizontal viscosity scheme used in the model
557 (D’Asaro, 1988; Morel and McWilliams, 2001; Akuetevi and Wirth, 2015;
558 Morel et al., 2019), but the general negative friction tendency term observed
559 here is mostly associated with the effect of the wind stress and its redistribu-
560 tion in the surface mixed layer. Testing with different anticyclonic vortices
561 or changing the selection criterion for retained particles does not change our
562 findings.

563 As can be seen from Fig. 4, cyclonic structures are constituted of smaller
564 vortices and filaments and are not as well defined as anticyclonic vortices. It
565 is also clear that the interior PV anomaly is mostly negative, although there
566 are traces of positive PV anomalies on vertical sections (third column in Fig.
567 3). We selected areas with positive interior PV anomalies and performed
568 a similar Lagrangian analysis (again backtracking 2000 particles from day
<figure><img src="image_15.png" /><figcaption>Figure 11: Mean PV evolution (upper panel) and evolution rate (lower panel) for all 1178 selected particles, constituting the anticyclonic eddy core (Fig. 7). For the upper panel, the blue dashed line corresponds to the PV anomaly of the model (calculated from the PV values of all selected particle at a given time from which we withdraw $f$), the green curve is obtained from the time integration of the friction tendency term and the red curve from the mixing tendency term, the black plain line is the sum of the integrated tendency terms and is to be compared to the blue dashed curve. The lower panel is similar but for the tendency terms (evolution rate).</figcaption></figure>
21). We found that the positive interior PV anomaly of particles constituting cyclonic eddies is determined during the last days of the evolution, i.e. after day 18, when the instability starts to develop. The criterion for particle selection had to be modified and the minimum correlation coefficient between model and diagnosed PV evolution was lowered to 0.7 to keep a few hundred particles. The mean PV anomaly of the particles constituting cyclonic eddies at day 21 was $\delta\overline{PV} \approx 1.5 \cdot 10^{-5}s^{-1}$, much lower than for anticyclonic structures (in absolute value). As mentioned above, the extension of the eddies is also much smaller, but it is also constituted of different subcores, i.e. disconnected parts associated with local maximum. In fact, eddies associated with positive vorticity or average PV anomalies are dominated by the surface outcropping signature, under which the positive interior PV anomaly eddies, generated during the development of the geostrophic instability, align. As seen before, mixing creates the positive PV anomalies at the base of the mixed layer, but during the late stage of the evolution, friction can also produce positive PV. Indeed, we have seen (see Fig. 10) that, during the upwelling development phase, the decrease of the stress from the surface and the upward bending of isopycnals is associated with a negative curl of stress along isopycnic levels, which explains the formation of negative PV anomalies by friction. If the isopycnal bends downward, which is the case when the instability develops and vortices emerge, the curl becomes positive and positive PV anomalies is then produced by friction (Morel et al., 2006).

## 4.2. Sensitivity to wind stress

We now evaluate the sensitivity of the PV change to the wind stress. As noted above, although we have considered a constant wind, the stress
along isopycnic levels is locally modified in the upwelling region, since it depends on the isopycnal position, velocity shear and local vertical viscosity coefficient (Morel et al., 2006, 2019).

We lowered the wind intensity with $W = 2 \text{ m/s}$ corresponding to a stress of $\tau \approx 0.006 \text{ N/m}^2$, 16 times smaller than before, and performed the same analysis on PV evolution. Figure 12 is similar to Fig. 3 and shows the evolution of vertically averaged upper-layer PV anomaly, surface vorticity, vertical section of PV anomaly, and vorticity along a South-North transect at $x = 205 \text{ km}$. Due to the lower forcing wind stress, mixing is less intense and the layer affected by mixing and wind stress is shallower. The upper layer over which PV is averaged is now bounded by the isopycnic level $\rho = 1024.6 \text{ kg/m}^3$ and the surface. The evolution is also slower than before with some peculiarities. First, Kelvin-Helmholtz like instabilities start at day 22 with several sequences until about day 37 when longer wavelength geostrophic instabilities develop, evolving into mesoscale structures after about 70 days. The layer PV and the surface vorticity structures are still strongly connected, but the patterns are different compared to the reference experiment with smaller and weaker anticyclonic vortices. The vertical sections of PV anomaly show little trace of positive PV anomalies in the interior during most of the evolution, and it appears only at the stage of fully developed turbulence (see the vertical section of PV at day 90, Fig. 12 last line, third row). Finally the isopycnic levels (black contours on the vertical sections) are less steep in the upwelling region, indicating that the fluid remains more stratified.

In Fig. 14 we plot the PV tendency terms associated with friction and mixing and their sum. Before the development of the geostrophic instability, PV anomalies generated by mixing show a different pattern for this lower
<figure><img src="image_16.png" /><figcaption>Figure 12: Evolution (days 10, 20, 35 for the experiment with a wind stress reduced to $\tau = 0.006 \, N/m^2$. Each day corresponding to a row) of horizontal maps of surface layer average PV anomaly (first column), surface vorticity (second column) and vertical sections of PV anomaly (third column) and vorticity (fourth column). The PV anomaly is obtained by withdrawing the reference PV, that is $f$. The vertical South-North sections are taken along the dashed line represented in the horizontal maps (at $x = 205 \, km$). Note that it crosses the core of an anticyclonic subsurface vortex at day 46. In rows 1, 3 and 4, the solid black contours are associated with isopycnic levels.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 13: Continued (days 46, 49 and 90).</figcaption></figure>
stress simulation, with sometimes a layer of negative PV production below
the positive one (row 1 and 2 of MixFlux, third column), which is inverted
compared to Fig. 5. The influence of the correction term, associated with
the evolution of the reference profile, is more important here: the presence of
the upwelling close to the coast generates a stronger stratification that is not
compensated by the shear, and mixing diminishes in the upwelling region
in comparison to the rearranged profile. This shows that nonlinearities and
details of the evolution, associated with the mixing closure scheme, can
lead to strong differences for the -rescaled- PV production and subsequent
dynamics.

As for the reference experiment, 2000 particles were seeded into the anti-cyclonic structure observed in the vertical section at day 46, whose center is located at x = 205 km and y = 242 km. These particles were backtracked to their initial positions while computing their PV and the PV tendency terms at their positions every 600 s. Using a similar selection criterion as for the reference experiment (correlation between the PV evolution and that reconstructed from the tendency terms above 0.9) about 1100 particles were retained. Figure 15 is similar to Fig. 11 and shows the evolution of the mean PV anomaly from the model (blue dashed line), the integration of the friction tendency term (green line), the mixing tendency term (red line) and the sum of both tendency terms (black line). There are striking differences with the reference experiment. First, both friction and mixing lead to negative PV production on average, but mixing has a modest effect; most of the variability is associated with friction, before the development of the geostrophic instability and vortex formation. The rate of PV production is also about 10 times weaker for friction or mixing. Another interesting aspect is that the analysis of individual particles (not shown) shows that after
<figure><img src="image_18.png" /></figure>

Figure 14: Evolution (days 10, 35, 46 and 90. Each day corresponding to a row) of vertical sections, taken along a South-North section located at x = 205 km, of the sum of diabatic tendency terms (first column), the frictional term (second column) and the diapycnal mixing term (third column) for the sensitivity test with a reduced wind stress ($\tau = 0.006 \text{ N/m}^2$).
<figure><img src="image_19.png" /><figcaption>Figure 15: Mean PV evolution (upper panel) and evolution rate (lower panel) of all 1083 selected particles forming the anticyclonic eddy core (Fig. 7) in the experiment with the wind stress reduced to $\tau = 0.006 \, N/m^2$. For the upper panel, the blue dashed line corresponds to the PV anomaly of the model (calculated from the PV values of all selected particles at a given time from which we remove $f$), the green curve is obtained from the time integration of the friction tendency term and the red curve from the mixing tendency term, the black solid line is the sum of the integrated tendency terms and is to be compared with the blue dashed curve. The lower panel is similar, but for the tendency terms (evolution rate).</figcaption></figure>
day 40, small scale processes spoil the PV budget for individual particles, but the mean evolution remains quite close to the estimated PV evolution (black and blue dashed curves in the upper panel of Fig. 15), whereas in the reference experiment, this leads to a bias.

## 5. Summary and Discussion

In the present paper, we advocate the use of a form of PV, rescaled by the time-evolving Lorenz's rearranged density profile (Lorenz, 1955; Naka-mura, 1995; Winters and D'Asaro, 1996). The proposed PV is the only form able to associate PV to the absolute vorticity for barotropic circulation. In general circumstances, it is also a generalization of the QG PV and is easily connected to the (cyclo)geostrophic circulation. We thus argue that this form of PV provides an objective evaluation of the effect of mixing and friction on the generation of PV. Its time evolution equation is derived and it is shown that global changes in stratification must be removed from the mixing production term, following isopycnal levels. The proposed form of PV and its evolution equation are valid in general (realistic) circumstances for the Navier-Stokes or Primitive equations. We believe this form has to be chosen when interpreting PV evolution in terms of dynamics (that is when using the inversion principle, relating PV to the circulation), especially to evaluate the respective influence of mixing and friction. Equation 6 shows that the objective rescaling of PV may lead to different interpretation depending on the characteristics of the water mass that is subject to diabatic changes: wind stress effect is rescaled by $Z(\rho)$ and mixing by $\partial_{\rho}Z(\rho)$. This can lead to strong differences with respect to the traditional Ertel PV form, for instance in region of ventilation, where a water mass is subject to heat
and momentum forcing at the surface before subducting and filling a deep layer.

Note that the evolution equation 6 is valid for other choices of the reference density profile, for instance using a "local"-but time evolving- reference density profile chosen at a fixed position. In this case, the interpretation of the PV evolution is different, relative to the PV and circulation fields at the chosen location. In the numerical configuration tested here, a choice of a reference profile chosen in the middle of the basin yielded similar results (not shown). In previous studies (Delpech et al., 2020; Assene et al., 2020; Aguedjou et al., 2021; Napolitano et al., 2022; Ernst et al., 2023) we chose a profile in the area of interest and fixed in time, arguing that, since we concentrated on the dynamics of deep layers, the reference stratification does not evolve much. Thus simpler choices for the reference profile are possible, but the ideal and rigorous choice remains Lorenz's rearranged profile.

The calculation of Lorenz's profile is generally straightforward and different methodologies have been proposed (Winters and D'Asaro, 1996; Tseng and Ferziger, 2001; Tailleux, 2013b). However, specific configurations may require adaptations, such as the periodic channel one used here. Other possible difficulties may arise from realistic configurations over a restricted region forced at its boundaries by fluxes from a larger domain. If the water mass characteristics entering the region varies, the interpretation of the PV evolution based on the rearranged profile in the region can be more difficult. For instance, in the present configuration, we mimicked open boundaries by restricting the calculation of the rearranged profile to the Northern half of the basin. We found this can lead to some biases in PV production by mixing. Indeed, because of transverse circulation associated with the upwelling development leads to light layers depletion and deep layers inflation, the
reference profile evolves even when there is no mixing. Interpretation of the
PV evolution must then be made with care.

We then numerically evaluate the use of this PV formulation using a
simplified configuration of an upwelling developing in a periodic channel. In
the reference upwelling simulation, with a strong mean wind, both diapycnal
mixing and friction influence the PV evolution of particles. Outcropping of
isopycnals at the surface provides the main reservoir for cyclonic eddy gen-
eration (Bretherton, 1966; Morel et al., 2006). Interior diabatic effects asso-
ciated with friction or mixing also create positive PV anomaly eddies, but
the latter are generated after the onset of the instability and remain weak.
They play little role in the early stages of the instability. On the contrary,
anticyclonic vortices are associated with the generation of interior negative
PV anomalies, which require non-conservative effects. Prior to the onset
of the geostrophic instability, the friction associated with the wind stress
generates negative PV anomalies and the interior PV is predominantly neg-
ative. Mixing generates both positive and negative PV anomalies, with a net
effect toward positive, only partially compensating the negative generation
by friction. This explains the initial PV structure of the upwelling and its
destabilization at geostrophic scales.

As the instability develops, the stratification and velocity fields become
complex and modulate the effect of mixing and friction at small scales,
leading to strong changes in the PV structure. At this stage, implicit dissi-
pation associated with numerical schemes becomes influential and spoils the
PV budget when using the explicit friction and mixing terms.

For a lower intensity of the wind stress, the effect of friction remains
similar with a negative production, but at a lower rate. On the other hand,
the effect of mixing is strongly modified: higher stratification associated with
the upwelling development under a weaker wind leads to less intense mixing in the upwelling region than in the open ocean and a net negative production rate for the rescaled PV. The global PV production due to mixing is an order of magnitude smaller than that due to friction. This shows the sensitivity of PV modification by non-conservative processes and their parameterizations.

A downwelling develops along the southern boundary and PV is also modified in this area, but the current remains stable in our configuration, so we did not focus on this area. In more realistic configurations, capes or topographic variations along the coast may generate localized eddies that are necessarily associated with the generation of PV anomalies. It is interesting to evaluate the influence of non-conservative processes in such a situation. Discrepancies between the characteristics of such eddies in models and in nature could well be associated with erroneous PV production associated with closure schemes.

More generally, the present study emphasizes that parameterizations play an important role in the determination of PV, the subsequent characteristics of the eddies and the structuring of the large scale circulation (which is partly determined by the redistribution of PV by the eddies). Testing different closure schemes, or parameters within a closure scheme, from the perspective of PV is therefore very informative. For example, recent studies have shown that accounting for the modulation of the wind stress by the surface ocean current feedback to the atmosphere drastically improves the surface eddy properties (Renault et al., 2016b, 2020). In upwelling systems, the wind drop-off near the coast (Capet et al., 2004; Renault et al., 2016a) or wind acceleration near capes modulates the surface stress, which, combined with the vertical redistribution of momentum by the turbulent fluxes, can lead to more complex effects of the wind stress (Kessouri et al.,
2022). It is also possible to refine the analysis and to distinguish specific non-conservative effects on PV evolution (e.g. separating vertical momentum fluxes from horizontal ones, or separating different diapycnal mixing processes). All of these mechanisms, and many more, could benefit from an analysis of PV generation associated with friction, mixing and their dependence on parameterization choices. The tools we have presented here can be very helpful for this perspective.

We have seen that when the dynamics becomes fully turbulent, there are discrepancies between the PV evolution and its reconstruction from the tendency terms. Temporal sub-sampling is one issue. It can be addressed by on-line Lagrangian diagnostics of PV, but the strategy for the initial particle positioning may be difficult, and backtracking of particles remains necessary to analyze the generation of PV for a particular region or eddy. We have also cited the imperfections of numerical schemes, which lead to non-conservative terms that can affect the PV evolution at the grid scale, when the dynamics generates small scale structures (filaments or other submesoscale eddies). If these uncertainties seem to cancel out when averaging many particles in the moderate wind experiment, we have seen that there is a bias possibly associated with this effect in the strong wind -reference- experiment. This can be problematic in achieving correct properties for eddies in numerical models, but again, the tools we have proposed can help to identify numerical effects and test the influence of improved numerical schemes.

Finally, Lorenz's profile corresponds to the minimum energy state which is linked to the available potential energy (APE). In some previous studies, we discussed a link between the PV structure of vortices or jets and their energy (Morel and McWilliams, 2001; Morel et al., 2019). The fact that, to get a rescaled PV related to the circulation, Lorenz's rearranged profile
is the ideal choice is thus probably not chance. The relationship between
rescaled PV and APE in general is thus worth investigating and some studies
dedicated to APE can be of great interest to substantiate the definition of the
ideal PV form and also possibly to identify other sources of PV modification.
For instance, Scotti and Passaggia (2019) have shown that the choice of
the Lorenz's rearranged profile to define APE can be justified because in
this case the effective energy does not depend on the time evolution of the
restratified flow, and argument that has strong similarities with the one used
here for the necessity of taking into account time evolution of the reference
profile for the rescaled PV. Finally, depending on the form of the equation
of state, the calculation of the Lorenz's profile is not straightforward, and
thermobaric effects can also make PV analysis more complicated (Straub,
1999). Recent studies in this field (Tailleux, 2013a,b; Saenz et al., 2015;
Tailleux, 2018) offer an interesting perspective for the use of PV analysis in
realistic configurations and are worth pursuing.

797 Acknowledgments

The authors thank contributions from their institutions : CNRS (French
National Center for Scientific Research), IRD (French National Research
Institute for Sustainable Development), University of Brest and University
of Toulouse. They also thank the 2 anonymous reviewers and Rémi Tailleux,
the editor in charge at Ocean Modelling, for their constructive remarks that
lead to a much better and more rigorous manuscript.
804

# Appendix A. Numerical calculation of the reference density profile

805

806

Calculating Lorenz's rearranged density profile can be efficiently done by computing the pdf of the density field (Tseng and Ferziger, 2001) or other methods proposed in the literature (see for instance Tailleux, 2013b). The pdf of density is easily calculated by defining fixed density classes and using the density values in existing numerical cells of the configuration. The obtained distribution should not suffer from biases when sigma levels and density fields vary over the domain because of a varying bottom topography and/or a complex -realistic- circulation is considered.

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

However, in the present configuration, we start from rest, with flat isopycnal over a flat topography. Discarding mixing, so that the initial linear reference profile should be maintained, the stratification evolution consists in bending isopycnals in narrow regions near the domain boundaries. In such a circumstance, the density pdf might be strongly biased, with large classes associated with density values over the middle of the domain and that are almost unchanged with respect to the initial state, and small classes associated with the slight variation near the boundaries. Such a pdf results in staircase like stratification instead of the linear one, which is very problematic for the evaluation of the rescaled PV. To avoid this problem, the existing methodologies to calculate the rearranged profile was adapted as follows:

825

826

827

828

1. Define the variation of volume and surface with depth, from the surface to the bottom of the domain, as a function of the vertical position $z$, with high enough discretization ($\Delta z \simeq 0.1 m$). For our flat bottom configuration, this is trivial.

829

2. Define intervals of density (density classes) covering the whole range
of densities present in the domain. The density range of each interval should be very small because we want to represent very small variations in the mixed layer. Here, to achieve good comparisons between direct calculation of the PV evolution rates and the one reconstructed from the tendency terms, we had to chose $\Delta\rho = 10^{-5} kg/m^{-3}$.

3. Calculate the density pdf: affect each density cell of the 3D domain to a density class and affect its volume to the class. This is the stage when adaptations are necessary to avoid biases for the present configuration. For each horizontal position, we interpolate the density profile on a virtual vertical grid with higher resolution. This allows to refine the pdf but is more computationally demanding. To achieve good results, the vertical step we had to choose for the present configuration is $\Delta z = 0.1 m$.

4. Calculate the net volume and mean density of each class and construct the initial rearranged 1D density profile filling the ocean from top and using the 1D volume and surface vectors determined at stage 1.

5. Interpolate the obtained profile using a regular density spacing (again a density step $\Delta\rho = 10^{-5} kg/m^{-3}$ is necessary here) to define the rescaling function at time $t Z(\rho, t)$.

## Appendix B. Generalized PV evolution equation

The traditional Ertel PV is defined as

$$PV_{Ertel} = -(\vec{\nabla} \times \vec{U} + \vec{f}) \cdot \vec{\nabla} \rho \qquad (B.1)$$

where $\vec{U}$ is the velocity field, $\vec{f}$ is the Earth rotation vector, whose projection on the local vertical axis defines the Coriolis parameter $f$, and $\rho$ is the
potential density. As already mentioned by Ertel (1942) (see also Muller, 2006; Morel et al., 2019), other forms for the potential vorticity are possible, with the same adiabatic conservation properties, obtained with

$$ 
\begin{aligned}
PV_{\text{rescaled}} &= (\vec{\nabla} \times \vec{U} + \vec{f}).\vec{\nabla}Z(\rho) \\
&= \vec{\nabla}.((\vec{\nabla} \times \vec{U} + \vec{f}) Z(\rho))
\end{aligned}
\quad (B.2) $$

where $Z(\rho)$ can be any function of potential density alone.

For geophysical flows, isopycnic variations of PV can be inverted to determine the geostrophic currents and associated stratification (Hoskins et al., 1985). But the Ertel PV is not the best PV expression to do so, in particular, for small geostrophic perturbations relative to a state at rest, quasi-geostrophic PV (Pedlosky, 1987; Muller, 2006; Cushman-Roisin and Beckers, 2011) should be recovered. Indeed the quasigeostrophic PV is given by the sum of the Coriolis parameter and the 3D Laplacian of the pressure variations associated with the small perturbation, from which both the geostrophic current and stratification can be obtained. With this straightforward link between PV and dynamical fields, the quasigeostrophic physics, associated with conservation and inversion of quasigeostrophic PV, provides important insights and is a fundamental approach for the understanding of the ocean dynamics from mesoscale to large scale (Pedlosky, 1987; Cushman-Roisin and Beckers, 2011). It is also at the base of the "PV thinking" approach (see in particular Hoskins et al., 1985, among many other) based on the analysis of the evolution of the PV field.

As discussed in the introduction, to overcome some problems associated with $PV_{Ertel}$ and simplify its interpretation, Morel et al. (2019) proposed $Z(\rho^*) = z$ (where $z$ is the vertical coordinate oriented upward), a rescaling function defined using a reference density profile $\rho^*(z)$ representing the
stratification at rest (see also Delpech et al., 2020; Assene et al., 2020). It is rigorously obtained using Lorenz (1955) diabatic rearrangement, where each fluid particle is classified according to its potential density and the domain is gradually uniformly filled following this stable rearrangement (see also Nakamura, 1995; Winters and D'Asaro, 1996). For such a choice, the rescaled PV at rest is $PV_{rescaled}^{rest} = f$. The calculation of the reference profile can be cumbersome, and in general taking into account a profile which is typical of the stratification of the studied region can be taken as a reference to rescale the PV.

In the present paper, to study the vortical dynamics of the upper layer in a context of rapid diabatic evolution associated with wind stress and mixing, the evolution of the reference stratification has to be taken into account.

We thus define a general rescaled PV field as

$$PV_{rescaled} = \vec{\nabla}.((\vec{\nabla} \times \vec{U} + \vec{f}) Z(\rho, t)) \quad (B.3)$$

where $Z(\rho, t)$ is now a function of both potential density and time and associated with a reference profile $\rho^*(z, t)$ that is time dependent. At any time, $Z(\rho^*(z, t), t) = z$, so that the rescaled PV at rest is still given by the Coriolis parameter $PV_{rescaled}^{rest} = f$.

The evolution equation for $PV_{rescaled}$ is obtained from Eq. 1 and, following Muller (2006), we get (see also Morel et al., 2019, for the specific divergence form we chose here)

$$\frac{d}{dt}PV_{rescaled} = \vec{\nabla}.((\vec{\nabla} \times \vec{F}) Z(\rho, t) + (\vec{\nabla} \times \vec{U} + \vec{f}) \frac{d}{dt}[Z(\rho, t)]) \quad (B.4)$$

The last term can be rewritten using

$$\frac{d}{dt}Z(\rho, t) = \partial_{\rho}Z(\rho, t) \dot{\rho} + \partial_{t}Z(\rho, t) \quad (B.5)$$
We now take into account the specific form chosen for the function $Z(\rho, t)$, which verifies $Z(\rho^*(z, t), t) = z$. Differentiating the latter with respect to time yields

$$ \partial_{\rho} Z(\rho^*(z, t), t) \partial_t \rho^*(z, t) + \partial_t Z(\rho^*(z, t), t) = 0 \quad (B.6) $$

or equivalently

$$ \partial_t Z(\rho^*, t) = -\partial_\rho Z(\rho^*, t) \partial_t \rho^* \quad (B.7) $$

from which we get

$$ \partial_t Z(\rho, t) = -\partial_\rho Z(\rho, t) \partial_t \rho^*(Z(\rho, t), t) \quad (B.8) $$

Using Eq. B.8 into B.5 and finally replacing terms in B.4 we get

$$ \frac{d}{dt} PV_{\text{rescaled}} = \vec{\nabla}.((\vec{\nabla} \times \vec{F}) Z(\rho, t) + (\vec{\nabla} \times \vec{U} + \vec{f}) \partial_{\rho} Z(\rho, t) (\dot{\rho} - \partial_t \rho^*(Z(\rho, t), t))) \quad (B.9) $$

or in Eulerian form

$$ \begin{aligned} \partial_t PV_{\text{rescaled}} = \vec{\nabla}.(&\vec{U} PV_{\text{rescaled}} \\ &+ (\vec{\nabla} \times \vec{F}) Z(\rho, t) \\ &+ (\vec{\nabla} \times \vec{U} + \vec{f}) \partial_{\rho} Z(\rho, t) (\dot{\rho} - \partial_t \rho^*|_{\rho,t})) \end{aligned} \quad (B.10) $$

<footnote n="887">887 where the first term of the right hand side divergence is associated with adiabatic advection, the second term with friction and the third one with diapycnal mixing. For the latter term, a correction is made for the generalized rescaled PV we propose here: the diabatic mixing term has to be corrected and the evolution of the reference profile $\partial_t \rho^*|_{\rho,t} = \partial_t \rho^*(Z(\rho, t), t)$ has to be withdrawn as it has no consequence on the generation of "dynamical" PV anomalies. Note it has to be evaluated using the reference profile</footnote>
time evolution at a vertical level corresponding to the density value in the physical domain (not the elevation).

## Appendix C. Calculation of PV and tendency terms

Morel et al. (2019) showed that the divergence form of PV (Eq. 5) preserves budgets and drastically simplifies the numerical calculation of PV, which can be expressed with a compact scheme, i.e. using a single PV grid cell on a 3D Arakawa C-grid (see Fig. C.16 and Arakawa and Lamb, 1977). As mentioned above, if we note the similarity between the divergence expression of PV (Eq. 5) and the diabatic PV tendency terms of its evolution equation (Eq. 6), we can see that the same scheme can be used to calculate the friction and mixing tendency terms.

Indeed, all terms can be written as

$$Div^{PV} = \nabla \cdot (\nabla \times \vec{F}^{vel} Z^{\rho}) \quad (C.1)$$

where $Div^{PV}$ is the divergence calculated at PV points of the staggered C-grid (see Fig. C.16), $Z^{\rho}$ is a function of density calculated at density points and $\nabla \times \vec{F}^{vel}$ is the curl of a 3D vector $\vec{F}^{vel} = (F^{U}, F^{V}, F^{W})$ whose components are located at $(U, V, W)$ points.

Integrating Eq. C.1 over the PV cell, whose corners are located at density points (see Fig. C.16), we get

$$\delta V^{PV} Div^{PV} = \iint_{S} Z^{\rho} \nabla \times \vec{F}^{vel} \cdot d\vec{S} \quad (C.2)$$

where $\delta V^{PV}$ is the volume of the PV cell and the right hand side integral is the flux of $Z^{\rho} \nabla \times \vec{F}^{vel}$ though all sides of the cell. Using the Stokes circulation theorem, the calculation of the latter term is simplified for the
<figure><img src="image_20.png" /><figcaption>Figure C.16: Elementary cell, for a 3D C-grid, used for the calculation of PV and the tendency terms. We consider Cartesian coordinates (x, y, z) associated with indices (i, j, k). Note that the calculation of the flux through each side of the cell is simplified for the C-grid because we can use the Stokes circulation theorem.</figcaption></figure>
C-grid where velocity points are located at the center of edges parallel to the velocity component (see Fig. C.16). For instance, the flux through the side given in Fig. C.16 (right panel) is

$$ \left[ \iint_S Z^\rho \vec{\nabla} \times \vec{F}^{\text{vel}} \cdot d\vec{S} \right]_{i,j,k}^x = \overline{Z^\rho}_{i,j,k} \cdot \left( F_{i,j-1,k}^W \delta z - F_{i,j,k}^W \delta z + F_{i,j,k}^V \delta y - F_{i,j,k-1}^V \delta y \right) \quad (C.3) $$

with

$$ \overline{Z^\rho}_{i,j,k} = \frac{Z_{i,j,k}^\rho + Z_{i,j,k-1}^\rho + Z_{i,j-1,k}^\rho + Z_{i,j-1,k-1}^\rho}{4} \quad (C.4) $$

<footnote n="909">909 and so on for the fluxes through other sides.</footnote>
910 Bibliography

911 Aguedjou, H. M. A., Chaigneau, A., Dadou, I., Morel, Y., Pegliasco, C., Da-
912 Allada, C. Y., Baloïtcha, E., 2021. What can we learn from observed tem-
913 perature and salinity isopycnal anomalies at eddy generation sites? ap-
914 plication in the tropical atlantic ocean. Journal of Geophysical Research:
915 Oceans 126 (11), e2021JC017630.

916 Akuetevi, C. Q. C., Wirth, A., 2015. Dynamics of turbulent western-
917 boundary currents at low latitude in a shallow-water model. Ocean Science
918 11 (3), 471–481.

919 Arakawa, A., Lamb, V., 1977. Computational design of the basic dynamical
920 processes of the ucla general circulation model. Methods Comput. Phys.
921 17, 174–267.

922 Assene, F., Morel, Y., Delpech, A., Aguedjou, M., Jouanno, J., Cravatte, S.,
923 Marin, F., Ménesguen, C., Chaigneau, A., Dadou, I., Alory, G., Holmes,
924 R., Bourlès, B., Koch-Larrouy, A., Aug 2020. From mixing to the large
925 scale circulation: How the inverse cascade is involved in the formation of
926 the subsurface currents in the gulf of guinea. Fluids 5 (3), 147.

927 Benthuysen, J., Thomas, L., 2012. Friction and diapycnal mixing at a slope:
928 Boundary control of potential vorticity. J. Phys. Oceanogr. 42, 1509–1523.

929 Benthuysen, J., Thomas, L., 2013. Nonlinear stratified spindown over a
930 slope. J. Fluid Mech. 726, 371–403.

931 Bretherton, F., 1966. Critical layer instability in baroclinic flows. Q.J.R.
932 Meteorol. Soc. 92 (2), 325–334.
933 Capet, X., Marchesiello, P., McWilliams, J. C., 2004. Upwelling response to
934 coastal wind profiles. Geophysical Research Letters 31 (13).

935 Capet, X., McWilliams, J. C., Molemaker, M. J., Shchepetkin, A. F.,
936 2008a. Mesoscale to submesoscale transition in the California Current sys-
937 tem. part i: Flow structure, eddy flux, and observational tests. J. Phys.
938 Oceanogr. 38 (1), 29-43.

939 Capet, X., McWilliams, J. C., Molemaker, M. J., Shchepetkin, A. F., 2008b. Mesoscale to submesoscale transition in the California Current system. part ii: Frontal processes. J. Phys. Oceanogr. 38 (1), 44-64.

942 Charney, J., Stern, M., 1962. On the stability of internal baroclinic jets in
943 a rotating atmosphere. J. Atmos. Sci. 19 (2), 159-172.

944 Cushman-Roisin, B., Beckers, J.-M., 2011. Introduction to Geophysical
945 Fluid Dynamics. Academic press, 875 pp.

946 Czaja, A., Hausmann, U., 2009. Observations of entry and exit of potential
947 vorticity at the sea surface. J. Phys. Oceanogr. 39, 2280-2294.

948 D'Asaro, E., 1988. Generation of submesoscale vortices: A new mechanism.
949 J. Geophys. Res.: Oceans. 93-C6, 2156-2202.

950 Debreu, L., Marchesiello, P., Penven, P., Cambon, G., 2012. Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean Modelling 49, 1-21.

953 Delpech, A., Cravatte, S., Marin, F., Morel, Y., Gronchi, E., Kestenare, E., 2020. Observed tracer fields structuration by middepth zonal jets in the tropical pacific. Journal of Physical Oceanography 50 (2), 281-304.
956 Early, J. J., Lelong, M. P., Sundermeyer, M. A., 2021. A generalized wave-
957 vortex decomposition for rotating boussinesq flows with arbitrary strati-
958 fication. Journal of Fluid Mechanics 912, A32.

959 Ernst, P. A., Subrahmanyam, B., Morel, Y., Trott, C. B., Chaigneau, A., 2023. Subsurface eddy detection optimized with potential vorticity from models in the arabian sea. Journal of Atmospheric and Oceanic Technology 40 (6), 677–700.

963 Ertel, H., 1942. On hydrodynamic eddy theorems. Physikalische Zeitschrift 43, 526-529.

965 Giordani, H., Lebeaupin Brossier, C., Léger, F., Caniaux, G., 2017. A PV-approach for dense water formation along fronts: Application to the NorthWestern Mediterranean. J. Geophys. Res.: Oceans 122 (2), 995–1015.

969 Gula, J., Blacic, T. M., Todd, R. E., 2019. Submesoscale coherent vortices in the Gulf Stream. Geophys. Res. Lett. 46.

971 Gula, J., Molemaker, M., McWilliams, J., 2015. Topographic vorticity gen-
eration, submesoscale instability and vortex street formation in the Gulf
Stream. Geophys. Res. Lett. 42, 4054–4062.

974 Gula, J., Molemaker, M., McWilliams, J., 2016. Topographic generation of
975 submesoscale centrifugal instability and energy dissipation. Nature Com-
munications 7, 12811.

977 Hallberg, R., Rhines, P., 1996. Buoyancy-driven circulation in an Ocean
978 basin with isopycnals intersecting the sloping boundary. J. Phys.
979 Oceanogr. 26, 913–940.
980 Hallberg, R., Rhines, P., 2000. Boundary sources of potential vorticity in geophysical circulations. Developments in Geophysical Turbulence, R. M. Kerr and Y. Kimura, Eds., Kluwer Academic, 51–65.

983 Haynes, P., McIntyre, M., 1987. On the evolution of vorticity and potential vorticity in the presence of diabatic heating and frictional or other forces. J. Atmos. Sci. 44 (5), 828–841.

986 Haynes, P., McIntyre, M., 1990. On the conservation and impermeability theorems for potential vorticity. J. Atmos. Sci. 47 (16), 2021–2031.

988 Held, I., Pierrehumbert, R., Garner, S., Swanson, K., 1995. Surface quasi-geostrophic dynamics. J. Fluid Mech. 282, 1–20.

990 Herbette, S., Morel, Y., Arhan, M., 2003. Erosion of a surface vortex by a seamount. J. Phys. Oceanogr. 33, 1664–1679.

992 Herbette, S., Morel, Y., Arhan, M., 2004. Subduction of a surface vortex under an outcropping front. J. Phys. Oceanogr. 34, 1610–1627.

994 Herbette, S., Morel, Y., Arhan, M., 2005. Erosion of a surface vortex by a seamount on the beta plane. J. Phys. Oceanogr. 35 (11), 2012–2030.

996 Holland, W. R., Keffer, T., Rhines, P., 1984. Dynamics of the Ocean in general circulation: The potential vorticity field. Nature 308, 698–705.

998 Hoskins, B. J., McIntyre, M. E., Robertson, A. W., 1985. On the use and significance of isentropic potential vorticity maps. Q. J. Roy. Met. Soc. 470, 877–946.

1001 Isern-Fontanet, J., Chapron, B., Lapeyre, G., Klein, P., 2006. Potential
1002 use of microwave sea surface temperatures for the estimation of ocean currents. Geophys. Res. Let. 33, L24608.
1003

1004 Kessouri, F., Renault, L., McWilliams, J. C., Damien, P., Bianchi, D., 2022. Enhancement of oceanic eddy activity by fine-scale orographic winds drives high productivity, low Oxygen, and low pH conditions in the Santa Barbara channel. J. Geophys. Res. : Oceans 127 (12), e2022JC018947.
1005
1006
1007

1008 Lapeyre, G., 2017. Surface quasi-geostrophy. Fluids 2 (1).
1009

1010 Lapeyre, G., Klein, P., Hua, B., 2006. Oceanic restratification forced by surface frontogenesis. J. Phys. Oceanogr. 36 (8), 1577–1590.

1011 Large, W., McWilliams, J., Doney, S., 1994. Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Rev. Geophys. 32, 363–403.
1012
1013

1014 Le Hénaff, M., Kourafalou, V., Morel, Y., Srinivasan, A., 2012. Simulating the dynamics and intensification of cyclonic Loop Current frontal eddies in the Gulf of Mexico. J. Geophys. Res. 117 (C02034).
1015
1016

1017 Legg, S., Jones, H., Visbeck, M., 1996. A Heton perspective of baroclinic eddy transfer in localized open ocean convection. J. Phys. Oceanogr. 26 (10), 2251–2266.
1018
1019

1020 Legg, S., Marshall, J., 1993. A Heton model of the spreading phase of open-ocean deep convection. J. Phys. Oceanogr. 23 (6), 1040–1056.
1021

1022 Lorenz, E. N., 1955. Available potential energy and the maintenance of the general circulation. Tellus 7 (2), 157–167.
1023
1024 Luyten, J. R., Pedlosky, J., Stommel, H., 1983. The ventilated thermocline. J. Phys. Oceanogr. 13 (2), 292-309.

1026 Marchesiello, P., McWilliams, J. C., Shchepetkin, A., 2003. Equilibrium structure and dynamics of the California current system. Journal of Physical Oceanography 33 (4), 753-783.

1029 Marshall, J., Nurser, G., 1992. Fluid dynamics of oceanic thermocline ventilation. J. Phys. Oceanogr. 22, 583-595.

1031 McWilliams, J., Flierl, G., 1979. Evolution of isolated, non-linear vortices. J. Phys. Oceanogr. 9 (6), 1155-1182.

1033 McWilliams, J. C., 1984. The emergence of isolated coherent vortices in turbulent flow. J. Fluid Mech. 146, 21-43.

1035 Meunier, T., Rossi, V., Morel, Y., Carton, X., 2010. Influence of bottom topography on an upwelling current: Generation of long trapped filaments. Ocean Modell. 41, 277-303.

1038 Molemaker, M. J., McWilliams, J. C., Dewar, W. K., 2015. Submesoscale instability and generation of mesoscale anticyclones near a separation of the California Undercurrent. J. Phys. Oceanogr. 45 (3), 613-629.

1041 Morel, Y., Darr, D., Tailandier, C., 2006. Possible sources driving the potential vorticity structure and long-wave instability of coastal upwelling and downwelling currents. J. Phys. Ocean. 36, 875-896.

1044 Morel, Y., Gula, J., Ponte, A., 2019. Potential vorticity diagnostics based on balances between volume integral and boundary conditions. Ocean Modelling 138, 23-35.
1047 Morel, Y., McWilliams, J., 1997. Evolution of isolated interior vortices in the Ocean. J. Phys. Ocean. 27 (5), 727-748.

1049 Morel, Y., McWilliams, J., 2001. Effects of isopycnal and diapycnal mixing on the stability of Ocean ic currents. J. Phys. Ocean. 31, 2280-2296.

1051 Morel, Y., Thomas, L., 2009. Ekman drift and vortical structures. Ocean Modell. 27, 185-197.

1053 Muller, P., 2006. The Equations of Oceanic Motions. Cambridge University Press.

1055 Nakamura, N., 1995. Modified Lagrangian-mean diagnostics of the stratospheric Polar vortices. part i. formulation and analysis of GFDL SKYHI GCM. Journal of Atmospheric Sciences 52 (11), 2096-2108.

1058 Napolitano, D. C., Alory, G., Dadou, I., Morel, Y., Jouanno, J., Morvan, G., 2022. Influence of the gulf of guinea islands on the atlantic equatorial undercurrent circulation. Journal of Geophysical Research: Oceans 127 (9), e2021JC017999.

1062 Pedlosky, J., 1987. Geophys. Fluid Dyn. Springer, New York, 710 pp.

1063 Renault, L., Hall, A., McWilliams, J. C., 2016a. Orographic shaping of U.S. West Coast Wind Profiles During the Upwelling Season. Climate Dynamics, 1-17.

1066 Renault, L., Masson, S., Arsouze, T., Madec, G., McWilliams, J. C., 2020. Recipes for how to force oceanic model dynamics. Journal of Advances in Modeling Earth Systems 12 (2), e2019MS001715.
1069 Renault, L., Molemaker, M. J., Mcwilliams, J. C., Shchepetkin, A. F., Lemarié, F., Chelton, D., Illig, S., Hall, A., 2016b. Modulation of wind work by oceanic current interaction with the atmosphere. Journal of Physical Oceanography 46 (6), 1685–1704.

1073 Rhines, P., 1986. Vorticity dynamics of the Ocean in general circulation. Annual review of Fluid Mechanics 18, 433–497.

1075 Rhines, P., Young, W., 1982a. Homogenization of potential vorticity in planetary gyres. J. Fluid Mech. 122, 347–367.

1077 Rhines, P., Young, W., 1982b. A theory of the wind-driven circulation. I, mid-Ocean Gyres. J. Mar. Res. 40, 559–596.

1079 Ripa, P., 1991. General stability conditions for a multi-layer model. J. Fluid Mech. 222, 119–137.

1081 Roed, L. P., Shi, X. B., 1999. A numerical study of the dynamics and energetics of cool filaments, jets, and eddies off the iberian peninsula. Journal of Geophysical Research: Oceans 104 (C12), 29817–29841.

1084 Rossi, V., Morel, Y., Garcon, V., 2010. Effect of the wind on the shelf dynamics: formation of a secondary upwelling along the continental margin. Ocean Modell. 31, 51–79.

1087 Saenz, J. A., Tailleux, R., Butler, E. D., Hughes, G. O., Oliver, K. I. C., 2015. Estimating Lorenz’s reference state in an ocean with a nonlinear equation of state for seawater. Journal of Physical Oceanography 45 (5), 1242–1257.

1091 Schneider, T., Held, I., Garner, S., 2003. Boundary effects in potential vorticity dynamics. J. Atmos. Sci. 60 (8), 1024–1040.
1093 Scotti, A., Passaggia, P.-Y., 2019. Diagnosing diabatic effects on the avail-
1094 able energy of stratified flows in inertial and non-inertial frames. Journal
1095 of Fluid Mechanics 861, 608-642.

1096 Shchepetkin, A., McWilliams, J., 2005. The Regional Oceanic Modeling Sys-
tem (roms): a split-explicit, free-surface, topography-following-coordinate
1097 Oceanic model. Ocean Modell. 9, 347-404.
1098

1099 Shchepetkin, A. F., McWilliams, J. C., 1998. Quasi-monotone advection
1100 schemes based on explicit locally adaptive dissipation. Monthly Weather
1101 Review 126 (6), 1541 – 1580.

1102 Song, Y., Haidvogel, D., 1994. A semi-implicit ocean circulation model using
1103 a generalized topography-following coordinate system. Journal of Compu-
tational Physics 115 (1), 228-244.
1104

1105 Straub, D. N., 1999. On thermobaric production of potential vorticity in the
1106 ocean. Tellus A: Dynamic Meteorology and Oceanography.

1107 Sutyrin, G., Flierl, G., 1994. Intense vortex motion on the beta-plane -
1108 Development of the beta-gyres. J. Atmos. Sci. 51 (5), 773-790.

1109 Tailleux, R., 2013a. Available potential energy and exergy in stratified fluids.
1110 Annual Review of Fluid Mechanics 45 (1), 35-58.

1111 Tailleux, R., 2013b. Available potential energy density for a multicomponent
1112 Boussinesq fluid with arbitrary nonlinear equation of state. Journal of
1113 Fluid Mechanics 735, 499-518.

1114 Tailleux, R., 2018. Local available energetics of multicomponent compress-
1115 ible stratified fluids. Journal of Fluid Mechanics 842, R1.
1116 Talley, L. D., 1988. Potential vorticity distribution in the North Pacific. J. Phys. Oceanogr. 18 (1), 89-106.

1118 Taylor, J., Ferrari, R., 2010. Buoyancy and wind-driven convection at mixed layer density fronts. J. Phys. Oceanogr. 40, 1222-1242.

1120 Thomas, L., Taylor, J., Ferrari, R., 2013. Symmetric instability in the Gulf Stream. Deep-Sea Res. 91, 96-110.

1122 Thomas, L. N., 2005. Destruction of potential vorticity by winds. J. Phys. Oceanogr. 35 (12), 2457-2466.

1124 Tseng, Y. H., Ferziger, J. H., 2001. Mixing and available potential energy in stratified flows. Phys. Fluids 13, 1281-1293.

1126 Vic, C., Roullet, G., Capet, X., Carton, X., Molemaker, M. J., Gula, J., 2015. Eddy-topography interactions and the fate of the Persian Gulf Outflow. J. Geophys. Res.: Oceans 120 (10), 6700-6717.

1129 Wagner, G. L., Young, W. R., 2015. Available potential vorticity and wave-averaged quasi-geostrophic flow. Journal of Fluid Mechanics 785, 401-424.

1131 Wenegrat, J. O., Thomas, L. N., Gula, J., McWilliams, J. C., 2018. Effects of the submesoscale on the potential vorticity budget of Ocean Mode Waters. J. Phys. Oceanogr. 48 (9), 2141-2165.

1134 Winters, K. B., D'Asaro, E. A., 1996. Diascalar flux and the rate of fluid mixing. Journal of Fluid Mechanics 317, 179-193.