# Eddy saturation in a reduced two-level model of the atmosphere

Melanie Kobrasa,b, Maarten H. P. Ambaumc and Valerio Lucarinia,b CONTACT
Melanie Kobras. Email<EMAIL_ADDRESS>aDepartment of Mathematics
and Statistics, University of Reading, Reading, UK; bCentre for the
Mathematics of Planet Earth, University of Reading, Reading, UK; cDepartment
of Meteorology, University of Reading, Reading, UK

###### Abstract

Eddy saturation describes the nonlinear mechanism in geophysical flows
whereby, when average conditions are considered, direct forcing of the zonal
flow increases the eddy kinetic energy, while the energy associated with the
zonal flow does not increase. Here we present a minimal baroclinic model that
exhibits complete eddy saturation. Starting from Phillips’ classical quasi-
geostrophic two-level model on the beta channel of the mid-latitudes, we
derive a reduced order model comprising of six ordinary differential equations
including parameterised eddies. This model features two physically realisable
steady state solutions, one a purely zonal flow and one where, additionally,
finite eddy motions are present. As the baroclinic forcing in the form of
diabatic heating is increased, the zonal solution loses stability and the eddy
solution becomes attracting. After this bifurcation, the zonal components of
the solution are independent of the baroclinic forcing, and the excess of heat
in the low latitudes is efficiently transported northwards by finite eddies,
in the spirit of baroclinic adjustment.

###### keywords:

eddy saturation; baroclinic instability; baroclinic adjustment; Phillips’
model; frictional adjustment

##  1 Introduction

The equilibrium volume transport of the Antarctic Circumpolar Current (ACC) is
known to be balanced by three contributions, namely the input of momentum at
the ocean surface by the wind, the downward transport of this momentum by
eddies and the bottom drag as the opposite force to the wind (Munk and Palmén
1951; Nadeau and Ferrari 2015). Straub (1993) was first to suggest that the
ACC volume transport is independent of the wind stress at the ocean surface
although the wind was assumed to be sufficiently strong. Since then more
studies, using resolved turbulent ocean eddies, have confirmed this finding
(e.g. Munday, Johnson, and Marshall (2013); Nadeau and Ferrari (2015)).
According to those studies, an increase in wind stress yields an increase of
the vertical shear in the channel and the flow becomes baroclinically
unstable, thus generating stronger eddies instead of a higher volume
transport, a process called eddy saturation.

Marshall et al. (2017), inspired by a previous model of variability in
atmospheric storm tracks (Ambaum and Novak 2014), explained the physical
principles of eddy saturation using a simple model with just three
ingredients, including a zonal momentum budget, a closure relation between the
eddy form stress and eddy energy, and an eddy energy budget. In their model,
both the vertical shear and the volume transport of the ACC are predicted to
be independent of the forcing by the surface wind and instead controlled by
the requirement of a sufficiently unstable vertical shear to overcome the
stabilizing role of the eddy energy dissipation. Moreover, the model explains
the increase of eddy energy with wind stress and they conclude an analogy to
the interaction between wave activity and baroclinicity in the original model
of atmospheric storm tracks (Ambaum and Novak 2014).

However, the process driving the variability of storm tracks in the mid-
latitude atmosphere differs from the transport in the ocean in several ways.
While mechanical stress dominates the oceanic processes described above, the
baroclinic instability in the atmosphere is a result of the horizontal
temperature gradient between the equator and poles, which is primarily due to
the presence of a stronger radiative forcing at low rather then high
latitudes. At all levels of the atmosphere, the meridional temperature
gradient is proportional to the vertical shear of the mean flow by thermal
wind balance and is the source of available potential energy for the eddies
which are associated with the storm tracks. The corresponding poleward eddy
heat fluxes in the storm tracks act to weaken the baroclinicity and therefore
the mean flow (e.g. Pedlosky (1979); Holton and Hakim (2013)). The latter
process is often referred to as baroclinic adjustment (Stone 1978). The Lorenz
energy cycle provides a comprehensive view of energetics of the climate system
that takes into account forcing, dissipation, and exchange of energy between
available potential and kinetic form, and between energy pertaining to the
mean flow and energy pertaining to eddy motions (Lorenz 1967; Peixoto and Oort
1992; Lucarini et al. 2014).

The heuristic model proposed by Ambaum and Novak (2014) describes this
baroclinic interaction between mean flow and eddy activity. In its steady
state, the model predicts a two-way equilibration of storm tracks to
extratropical thermal forcing and eddy friction: baroclinicity is independent
of the thermal forcing but proportional to the eddy dissipation, whereas storm
track activity is independent of eddy dissipation but proportional to thermal
forcing of large-scale baroclinicity and is therefore reminiscent of the eddy
saturation phenomenon in the ACC (Novak, Ambaum, and Harvey 2018).

Within the Lorenz energy cycle framework, one can see the eddy saturation
mechanism as that a forcing acting on the zonal fields results, on average, in
an increase in the reservoir of eddy energy (in both kinetic and available
potential form), whereas the reservoir of zonal energy stays largely
unaffected. Additionally, one can view the eddy saturation mechanism as an
extreme form of the baroclinic adjustment process; see also discussion in
Lucarini, Speranza, and Vitolo (2007).

The simplest model that can incorporate baroclinic processes together with
diabatic heating and surface friction is Phillips’ two-level quasi-geostrophic
model on the $\beta$-plane (Phillips 1956). The present paper uses this model
to derive a set of ordinary differential equations that are able to provide a
minimal yet meaningful model of the above described interaction between mean
flow and eddy activity. The novelty of this work is how the stability
properties of the derived models are determined. In contrast to the usual
approach of normal mode instability analysis (see e.g. Pedlosky (1979);
Hoskins and James (2014)), which focuses on studying the stability properties
of the zonal flow, we use here methods of dynamical systems theory, which
allow us to show the existence of a second attracting steady state instead of
growing normal mode baroclinic instabilities. Another novel aspect is that
this second steady state exhibits the above described eddy saturation
properties in a model with parameterised eddies.

In section 2 we briefly review the two-level model and its reduction which
follows work by Phillips (1956) and by Thompson (1987). In section 3 we
introduce the non-dimensionalisation used to define our reduced model and
determine the steady states. The stability and dependence on relevant
parameters thereof is analysed and physically explained in section 4 and
compared to the normal mode baroclinic instability analysis of Phillips’ two-
level model in section 5. In section 6 a quality factor is introduced to
describe the oscillatory behaviour of the model and section 7 summarizes the
results and discusses them in comparison to the current literature.

##  2 The model equations

We use the two-level quasi-geostrophic model in pressure coordinates of
Phillips (1956), consisting of two vorticity equations, coupled by a
thermodynamic energy equation and including surface friction and a
$\beta$-plane approximation, i.e. $f=f_{0}+\beta y.$ The equations are

$\displaystyle\left(\frac{\upartial}{\upartial
t}+\bm{V}_{1}\cdot\nabla\right)(\beta y+\nabla^{2}\psi_{1})-f_{0}$
$\displaystyle=A\nabla^{2}\nabla^{2}\psi_{1},$ (1a)
$\displaystyle\left(\frac{\upartial}{\upartial
t}+\bm{V}_{3}\cdot\nabla\right)(\beta
y+\nabla^{2}\psi_{3})+f_{0}\frac{\omega_{2}}{\delta p}$
$\displaystyle=A\nabla^{2}\nabla^{2}\psi_{3}-\kappa\nabla^{2}\psi_{4},$ (1b)
$\displaystyle\left(\frac{\upartial}{\upartial
t}+\bm{V}\cdot\nabla\right)\frac{\psi_{3}-\psi_{1}}{\delta
p}+\frac{R\overline{\Gamma}_{2}}{p_{2}f_{0}}\omega_{2}$
$\displaystyle=-\frac{R}{p_{2}f_{0}c_{\rm p}}J_{2},$ (1c)

where $\psi$ is the streamfunction, ${\bm V}=(u,v)$ the horizontal velocity
vector, $\omega={\rm d}p/{\rm d}t$ the vertical velocity and $\nabla$ is the
horizontal gradient. The pressure is denoted by $p$, where $\delta
p=500\,\text{hPa}$ is the pressure difference between the two vorticity levels
as well as the pressure-thickness of each layer. Subscripts 1, 2, 3 and 4
denote the pressure levels 250, 500, 750, and 1000 hPa respectively. The
fields are defined in the spatial domain $(x,y)\in[0,L]\times[0,W]$, where $L$
and $W$ are the length and width of the $\beta$-plane channel. We remark that,
as customary, the $x$ coordinate is aligned with the longitude and the $y$
coordinate is aligned with the latitude. The fields above are defined for non-
negative times $t\geq 0$ starting from smooth initial conditions. Clearly, we
have that $A(0,y,t)=A(L,y,t)$ $\forall y\in[0,W]$ and $\forall t\geq 0$ for
all fields $A$.

The parameters $A$ and $\kappa$ describe the eddy diffusion and the surface
friction diffusion. Additionally, $J$ is the diabatic heating and
$\overline{\Gamma}=RT/pc_{\rm p}-\upartial T/\upartial p$ the basic state
static stability, where $T$ is the temperature, $c_{\rm p}$ the isobaric
specific heat capacity of dry air and $R$ the specific gas constant for dry
air.

The diabatic heating $J$ is given by the sum of radiative and diffusive
contributions. The radiative contribution represents the net local radiative
energy gains and losses and is chosen to vary linearly in $y,$

$J_{\text{rad}}=2H\left(1-\frac{2y}{W}\right),$ (1b)

where $H$ is the mean rate of heating per unit mass for $y\in[0,W/2]$ (or
cooling for $y\in[W/2,W]$). The diffusive contribution represents the effect
of lateral eddy diffusion of temperature at level 2:

$J_{\text{diff}}=c_{\rm p}A\nabla^{2}T_{2}=\frac{2f_{0}c_{\rm p}p_{2}}{\delta
pR}A\nabla^{2}\psi_{\rm T}.$ (1c)

We remark that the parameters $A$ and $\kappa$ control the dissipation
processes that remove energy from the system. The parameter $H$ controls the
input of energy in the form of an increase in temperature difference between
low and high latitudes, hence fuelling the Lorenz energy cycle of the system
by production of zonal mean available potential energy. The competition
between dissipation and forcing determines the degree of instability and
turbulence of the system, which, after transients, reaches a statistically
steady state, the model climate.

Next, we define barotropic and baroclinic potential vorticities $q_{\rm
m}=\nabla^{2}\psi_{\rm m}+\beta y$ and $q_{\rm T}=(\nabla^{2}-\lambda^{2}_{\rm
R})\psi_{\rm T}$, where $\lambda^{2}_{\rm R}=2p_{2}f^{2}_{0}/(\delta
p^{2}\overline{\Gamma}_{2}R)$ is the inverse square of the Rossby radius of
deformation and the barotropic and baroclinic streamfunctions are defined as
$\psi_{\rm m}=(\psi_{1}+\psi_{3})/2$ and $\psi_{\rm T}=(\psi_{1}-\psi_{3})/2$.

Since we want to obtain evolution equations for the longitudinally averaged or
mean flow (denoted by overbars) and longitudinally varying disturbances or
eddies (denoted by primes) separately, we introduce the zonal mean
$\overline{A}(y,t)=1/L\int_{0}^{L}\mathrm{d}xA(x,y,t)$ and define the
deviation from such mean as $A^{\prime}(x,y,t)=A(x,y,t)-\overline{A}(y,t)$ for
any field $A$.

Now the evolution equations for the mean zonal wind and mean thermal wind or
shear are obtained by taking half the sum and the difference of equations (1a)
and (1b) where we introduced the geostrophic assumption of horizontally non-
divergent flow, yielding the geostrophic relation $u=-\upartial\psi/\upartial
y$ and $v=\upartial\psi/\upartial x$ for the velocities in all terms but the
one containing the Coriolis parameter. Substituting equation (1c), taking the
derivative with respect to $y$ and applying the definitions and assumptions
above, we arrive at

$\displaystyle\frac{\upartial\overline{u}_{\rm m}}{\upartial
t}-\frac{1}{\lambda_{y}^{2}}\frac{\upartial^{3}}{\upartial
y^{3}}\overline{v^{\prime}_{\rm m}u^{\prime}_{\rm
m}}-\frac{1}{\lambda_{y}^{2}}\frac{\upartial^{3}}{\upartial
y^{3}}\overline{v^{\prime}_{\rm T}u^{\prime}_{\rm
T}}=-\lambda_{y}^{2}A\overline{u}_{\rm m}-\frac{\kappa}{2}(\overline{u}_{\rm
m}-2\overline{u}_{\rm T}),$ (1da)
$\displaystyle\frac{\upartial\overline{u}_{\rm T}}{\upartial
t}-\frac{1}{\lambda_{y}^{2}+\lambda_{\rm R}^{2}}\frac{\upartial^{3}}{\upartial
y^{3}}\overline{v^{\prime}_{\rm T}u^{\prime}_{\rm
m}}-\frac{1}{\lambda_{y}^{2}+\lambda_{\rm
R}^{2}}\frac{\upartial^{3}}{\upartial y^{3}}\overline{v^{\prime}_{\rm
m}u^{\prime}_{\rm T}}-\frac{\lambda^{2}_{\rm R}}{\lambda^{2}_{y}+\lambda_{\rm
R}^{2}}\frac{\upartial^{2}}{\upartial y^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}$ (1db)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad=-\lambda^{2}_{y}A\overline{u}_{\rm
T}+\frac{\lambda_{y}^{2}}{\lambda_{y}^{2}+\lambda_{\rm
R}^{2}}\frac{\kappa}{2}(\overline{u}_{\rm m}-2\overline{u}_{\rm
T})+\frac{\lambda^{2}_{\rm R}}{\lambda_{y}^{2}+\lambda_{\rm
R}^{2}}\frac{2\delta pRH}{f_{0}p_{2}c_{\rm p}W},$

where $\overline{u}_{\rm m}$ is the mean zonal wind and $\overline{u}_{\rm T}$
the mean thermal wind or mean zonal shear. Here we additionally assumed a
specific meridional shape of the mean wind, namely
$\upartial^{2}\overline{u}/\upartial y^{2}=-\lambda_{y}^{2}\overline{u}$ where
$\lambda_{y}^{2}$ is an unspecified wavenumber.

Following Thompson (1987), we next restrict the eddy component of the
barotropic and baroclinic streamfunction to be of the form

$\psi^{\prime}_{\rm m,T}(x,y,t)=A_{\rm m,T}(t)\sin kx\sin ly+B_{\rm
m,T}(t)\cos kx\sin ly,$ (1e)

where the wavenumber $k$ remains unspecified and $l=\pi/W$. We note that
$\nabla^{2}\psi^{\prime}_{\rm m,T}=-\lambda^{2}_{\nabla}\psi^{\prime}_{\rm
m,T}$ with $\lambda^{2}_{\nabla}=k^{2}+l^{2}=k^{2}+\pi^{2}/W^{2}$. Using this
assumption, all eddy product terms in equations (1da) and (1db) besides the
net poleward eddy temperature flux $v^{\prime}_{\rm m}\psi^{\prime}_{\rm T}$
cancel. Using the explicit expression (1e) for $v^{\prime}_{\rm m}$ and
$\psi^{\prime}_{\rm T}$ we obtain

$\displaystyle\overline{v^{\prime}_{\rm m}\psi^{\prime}_{\rm T}}$
$\displaystyle=\frac{k}{2}(A_{\rm m}B_{\rm T}-A_{\rm T}B_{\rm
m})\sin^{2}\frac{\pi y}{W},$ (1fa)
$\displaystyle\frac{\upartial^{2}}{\upartial y^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}$ $\displaystyle=\frac{\pi^{2}k}{W^{2}}(A_{\rm
m}B_{\rm T}-A_{\rm T}B_{\rm m})\left(\cos^{2}\frac{\pi y}{W}-\sin^{2}\frac{\pi
y}{W}\right),$ (1fb)

so at $y=W/2$,

$\frac{\upartial^{2}}{\upartial y^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}=-\frac{2\pi^{2}}{W^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}.$ (1g)

We see that the net poleward heat transport attains its maximum at $y=W/2$ and
vanishes at $y=0$ and $y=W$. Indeed, this qualitatively corresponds to what is
observed in more complex versions of the same model (Lucarini, Speranza, and
Vitolo 2007) and to what is observed in the actual climate (Peixoto and Oort
1992). Following again Thompson (1987), we assume that the mean thermal wind
$\overline{u}_{\rm T}$ has the same meridional structure, and inspecting
equation (1da) we see that the same must be true for the mean zonal wind
$\overline{u}_{\rm m}$. Additionally, from the assumption about the zonal mean
wind, it implies that $\lambda^{2}_{y}=2\pi^{2}/W^{2}$. Finally, equations
(1da) and (1db) can be substantially simplified to

$\displaystyle\frac{\upartial\overline{u}_{\rm m}}{\upartial t}$
$\displaystyle=-\lambda_{y}^{2}A\overline{u}_{\rm
m}-\frac{\kappa}{2}(\overline{u}_{\rm m}-2\overline{u}_{\rm T}),$ (1ha)
$\displaystyle\frac{\upartial\overline{u}_{\rm T}}{\upartial t}$
$\displaystyle=-\frac{\lambda_{y}^{2}\lambda_{\rm
R}^{2}}{\lambda_{y}^{2}+\lambda_{\rm R}^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}-\lambda^{2}_{y}A\overline{u}_{\rm
T}+\frac{\lambda_{y}^{2}}{\lambda_{y}^{2}+\lambda_{\rm
R}^{2}}\frac{\kappa}{2}(\overline{u}_{\rm m}-2\overline{u}_{\rm
T})+\frac{\lambda^{2}_{\rm R}}{\lambda_{y}^{2}+\lambda_{\rm
R}^{2}}\frac{2\delta pRH}{f_{0}p_{2}c_{\rm p}W},$ (1hb)

which is now a closed system by an evolution equation for the net poleward
temperature flux.

This equation and evolution equations for three other occurring eddy
correlation terms, namely the mean meridional kinetic energy
$\overline{v^{\prime 2}_{\rm m}}$, the temperature variance
$\overline{v^{\prime 2}_{\rm T}}$ and the cross-correlation between
temperature and geopotential $\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm
T}}$, are derived following to some extend the technique by Thompson (1987),
with the difference that our model includes a surface friction term. The
poleward temperature flux equations is

$\displaystyle\frac{\upartial}{\upartial t}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}}$ $\displaystyle=\frac{\lambda^{2}_{\rm
R}+\lambda^{2}_{y}-\lambda^{2}_{\nabla}}{\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R}}\overline{u}_{\rm T}\overline{v^{\prime 2}_{\rm m}}-\frac{\lambda^{2}_{\rm
R}}{\lambda^{2}_{\nabla}(\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}\left(\beta+\lambda^{2}_{y}\overline{u}_{\rm
m}\right)\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm
T}}+\frac{\lambda^{2}_{\nabla}-\lambda^{2}_{y}}{\lambda^{2}_{\nabla}}\overline{u}_{\rm
T}\overline{v^{\prime 2}_{\rm T}}$
$\displaystyle\qquad-\left(2A\lambda^{2}_{\nabla}+\frac{\kappa(3\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}{2(\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}\right)\overline{v^{\prime}_{\rm m}\psi^{\prime}_{\rm T}}.$ (1i)

The evolution equations for the three other eddy correlation terms can be
found in appendix A.

##  3 The non-dimensional model

Up to this point the meridional wavenumber $k$ of the eddy component of the
streamfunctions remained unspecified. In principle, one has that $k=2p\pi/L$,
where $p$ is a non-vanishing natural number. Now, we let $k$ be a multiple of
the longitudinal wavenumber $l$ and we define $k=jl=j\pi/W$ where $j$ is a
parameter describing the aspect ratio of the eddies, that is – for simplicity
– assumed to be a non-negative number. Indeed, $j$ can only assume discrete
values, but since $L\gg W$, they are closely spaced and so the discrete nature
of $j$ is neglected in what follows. Hence, $j$ larger than 1 means that the
eddies are elongated in the meridional direction. If $j$ is smaller than 1 but
positive, the eddies are elongated in the longitudinal direction.

To simplify the notation, we introduce the dimensionless variables

$M=\frac{\lambda^{2}_{\rm R}}{2\beta}\overline{u}_{\rm m},\
S=\frac{\lambda^{2}_{\rm R}}{2\beta}\overline{u}_{\rm T},\
T=\frac{\lambda^{5}_{\rm R}}{2^{3}\beta^{2}}\overline{v^{\prime}_{\rm
m}\psi^{\prime}_{\rm T}},\ K=\frac{\lambda^{4}_{\rm
R}}{2^{2}\beta^{2}}\overline{v^{\prime 2}_{\rm m}},\ V=\frac{\lambda^{4}_{\rm
R}}{2^{2}\beta^{2}}\overline{v^{\prime 2}_{\rm T}},\ X=\frac{\lambda^{4}_{\rm
R}}{2^{2}\beta^{2}}\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm T}}$

and define the time variable $\tau=t\beta/\lambda_{\rm R}$, so that equations
(1ha), (1hb), (1i) and (1ua)-(1uc) take the form

$\displaystyle\frac{{\rm d}M}{{\rm d}\tau}$ $\displaystyle=-aM+\alpha S,$
(1ja) $\displaystyle\frac{{\rm d}S}{{\rm d}\tau}$
$\displaystyle=\frac{\gamma\alpha}{2}M-bS-4\gamma T+\eta H,$ (1jb)
$\displaystyle\frac{{\rm d}T}{{\rm d}\tau}$ $\displaystyle=\mu
SK-\frac{\delta\lambda^{2}_{\rm R}}{2\lambda_{y}^{2}}X-\delta MX-
cT+\frac{j^{2}-1}{j^{2}+1}SV,$ (1jc) $\displaystyle\frac{{\rm d}K}{{\rm
d}\tau}$
$\displaystyle=\frac{2^{2}j^{2}(j^{2}-1)}{j^{2}+1}\frac{\lambda^{2}_{y}}{\lambda^{2}_{\rm
R}}ST-dK+2\alpha X,$ (1jd) $\displaystyle\frac{{\rm d}V}{{\rm d}\tau}$
$\displaystyle=\frac{2^{2}j^{2}\lambda^{2}_{y}\mu}{\lambda^{2}_{\rm R}}ST-
eV+\alpha\zeta X,$ (1je) $\displaystyle\frac{{\rm d}X}{{\rm d}\tau}$
$\displaystyle=j^{2}\delta\left(1+\frac{2\lambda^{2}_{y}}{\lambda^{2}_{\rm
R}}M\right)T-cX+\alpha V+\frac{\alpha\zeta}{2}K,$ (1jf)

where the definitions of the dimensionless constants $\alpha$, $\gamma$,
$\delta$, $\epsilon$, $\nu$ and $\eta$ can be found in appendix B.

For a non-zero heating rate $H$ this system exhibits three steady states. The
first one is a zonal steady state where all eddy components are zero:

$P^{0}=(M^{0},S^{0},0,0,0,0)=\frac{2}{2ab-\gamma\alpha^{2}}\eta
H(\alpha,a,0,0,0,0).$ (1k)

Here, $M^{0}$ and $S^{0}$ are always positive because the denominator of the
pre-factor and all parameter values including the heating rate are positive
within the considered parameter space (see appendix C). The second steady
state, denoted by $P^{*}=(M^{*},S^{*},T^{*},K^{*},V^{*},X^{*})$, is given by

$\displaystyle M^{*}$ $\displaystyle=\frac{\alpha}{a}S^{*},$ $\displaystyle
S^{*}$
$\displaystyle=\frac{j^{2}a}{A}\left(\left(B^{2}+4CD\right)^{1/2}+B\right),$
$\displaystyle T^{*}$
$\displaystyle=-\frac{2ab-\gamma\alpha^{2}}{8a\gamma}S^{*}+\frac{\eta}{4\gamma}H,$
(1l)

where $K^{*}$, $V^{*}$, $X^{*}$ and the new parameters $B$, $C$ and $D$ can be
found in appendix D. Note that $B$, $C$ and $D$ are independent of the heating
$H$ and therefore the steady state mean zonal wind and shear are independent
of the heating.

The third steady state of the system can be excluded as unphysical since the
values of $K$ and $V$ are negative in the given parameter space despite being
non-dimensional variables for $\overline{v_{\rm m}^{\prime 2}}$ and
$\overline{v_{\rm T}^{\prime 2}}$, respectively, which being squared real
quantities must always be non-negative. Therefore the system exhibits two
physical solutions, in the following referred to as the zonal ($P^{0}$) and
the eddy steady state ($P^{*}$).

##  4 Physical mechanisms of zonal and eddy saturated steady state

The linear stability of these two solutions is determined by numerically
calculating the eigenvalues of the Jacobian matrix of system (1ja)-(1jf) at
each steady state. Therefore, the parameters in table C were kept fixed,
whereas the surface friction diffusion $\kappa$, the eddy diffusion $A$, the
mean rate of heating $H$ and the parameter $j$ describing the eddy aspect
ratio were varied within the ranges given in table C.

In the following the notion of an attracting (stable) and repelling (unstable)
steady state is used to describe the linear stability of the two steady
states. Hence, if the eddy steady state is attracting, the model converges to
the steady state with finite eddy contributions. This state then breaks the
zonal symmetry of the zonal state: the phase of the eddies is not determined
but the eddy correlation statistics are fixed in time. On the other hand, an
attracting zonal steady state describes a state of the atmosphere where eddy
activity decays until the flow is purely zonal, and once in that state remains
zonal.

The considered parameter space is divided into two opposite stability regions
where either the zonal or the eddy steady state is attracting and the
respective other state is repelling. Keeping $\kappa$, $A$ and $j$ fixed but
increasing the heating rate $H,$ a transcritical bifurcation occurs (see for
example Guckenheimer and Holmes (1983)) and the system switches stability from
the zonal to the eddy steady state. From equation (1k) it is clear that the
mean zonal wind and shear of $P^{0}$ grow linearly with the heating $H$.
However, this proportionality is lost for the eddy steady state, where $M^{*}$
and $S^{*}$ are independent of the heating rate. This can be seen directly
from equation (1l) and from the left hand side of figure 1, where contours of
the mean zonal wind at the top level for the respective attracting steady
state are shown. Such an insensitivity of the mean zonal wind to the forcing
by the heating rate closely resembles the eddy saturation mechanism discussed
in the introduction.

In contrast to that, the (non-dimensional) net poleward heat transport $T^{*}$
grows linearly with $H$ (see again equation (1l) and figure 1 on the right)
and is non-negative in its attracting parameter region only. For values of $H$
where the zonal steady state is attracting, $T^{*}$ is negative and therefore
the net eddy heat transport is equatorward for the eddy state. This would
imply that the system transports heat from cold towards warm regions, which
suggests a condition that is thermodynamically not realisable. This further
clarifies that the eddy steady state is physically irrelevant in its repelling
region. Furthermore, it is clear from equation (1l) that the steady state
value of the net poleward heat transport only depends on the shear and
incoming heat but not on the three other eddy correlation variables $K^{*}$,
$V^{*}$ and $X^{*}$. Hence, the dynamics of the system are resembled by the
mean zonal shear (proportional to the mean zonal wind) and the net poleward
heat transport, and $K^{*}$, $V^{*}$ and $X^{*}$ do not need to be considered
separately.

Figure 1:  Contours of the mean zonal wind at the top level $\overline{u}_{1}$
(left) and the net poleward eddy temperature flux (right) for the respective
attracting steady state; the dashed region is the parameter space where the
zonal steady state $P^{0}$ is attracting, the non-dashed region where the eddy
steady state $P^{*}$ is attracting; note that the range of $H$ in the right
figure is four times the range in the left; $\kappa=4\times 10^{-6}\text{
s}^{-1}$, $j=2.5$.

The above described stability switch is also dependent on the eddy diffusion
parameter $A$. This is illustrated in figure 1, where a higher eddy diffusion
yields a larger heating rate at which the zonal steady state loses stability.
The relation between eddy diffusion and heating rate in $P^{0}$ is determined
by equation (1k). Replacing $S^{0}$ by the mean poleward temperature gradient
at the middle level $\overline{T}_{2}$ yields

$G(A,\kappa)\left(-\frac{\upartial\overline{T}_{2}}{\upartial
y}\right)=\frac{4\lambda^{2}_{\rm R}}{c_{\rm
p}W\left(\lambda^{2}_{y}+\lambda^{2}_{\rm R}\right)}H\quad\text{with}\quad
G(A,\kappa)=\lambda^{2}_{y}A\frac{2\lambda^{2}_{y}A+(2\gamma+1)\kappa}{2\lambda^{2}_{y}A+\kappa}.$
(1m)

Since in the considered parameter space $\kappa$ is much smaller than one, the
function $G$ is approximately linear in $A$. In physical terms, for a fixed
equator to pole mean temperature gradient, a linear increase of the heating
rate $H$ can be balanced by a linearly increased eddy diffusion. In this
sense, the stabilization of the zonal state against heating is always realized
by an eddy heat transport, be it parameterised, as an eddy diffusion, or
explicit, as in eddy heat transport.

It is also clear from (1m) that the surface friction $\kappa$ is not able to
balance the heating and the stability analysis shows that a change in $\kappa$
does not influence the critical value of heating for realistic values of
$\kappa$ (see figure 3 on the left).

These opposed roles of the zonal and eddy steady state describe two physically
related mechanisms of the atmosphere to compensate incoming heat. In the
attracting region of the zonal steady state the relatively low heating can be
compensated by the eddy diffusion. However, this balance is no longer
achievable for increasing $H$ and decreasing $A$ so that in the eddy steady
state the eddies have to reach a finite size to perform poleward transport of
heat and compensate the atmospheric radiative energy imbalance. In the actual
atmosphere – or in more complex models – this corresponds to the case where
the flow is baroclinically unstable, and baroclinic cyclones grow and decay as
manifestation of an active Lorenz energy cycle and, at the same time,
transport heat poleward. Although it is a similar mechanism, the compensation
of the heating does not depend on the friction diffusion $\kappa$. This is
explained later in this section.

Another parameter that influences the stability regions of the steady states
is the aspect ratio $j\in{\mathbb{R}}^{+}$ of the eddy shape. For the previous
analysis $j=2.5$ was chosen, which corresponds to eddies elongated by a factor
of $2.5$ in the meridional direction. For this aspect ratio, and even more
meridionally elongated eddies, a stability switch always occurs and the larger
$j$, the smaller is the heating rate $H$ at this switching point. In contrast
to that, for $j$ smaller than one, corresponding to longitudinally elongated
eddies, there is no exchange of stability and the zonal steady state is
attracting for the whole parameter space under consideration.

For an aspect ratio above one, the value of $H$ at which the zonal state loses
stability is additionally dependent on the eddy diffusion parameter. The left
hand side of figure 2 shows the values of $H$ at which the switch occurs as a
function of $j$ for several orders of magnitude of $A$. For a vanishing small
eddy diffusion (e.g. $A$ on the order of $10\text{ m}^{2}\text{ s}^{-1}$) the
eddy steady state is attracting for the whole range of $H$ and an aspect ratio
larger than one. Increasing $A$ yields an exchange of stability for higher
$H$, which was already seen in figure 1, and additionally an increasing value
of $j$ is required for the eddy steady state to become attracting.

Hence, for meridionally more elongated eddies the eddy steady state becomes
attracting at a lower heating rate. It is indeed well known that, in the full
quasi-geostrophic flow, baroclinic instability is facilitated for these
geometrical conditions (Pedlosky 1979). Only for a very large aspect ratio $j$
above five and for a large eddy diffusion the critical value of $H$ increases
again. As before, besides this behaviour at very high aspect ratios, a change
in the surface friction does not change this critical heating rate.

Figure 2:  The value of the heating rate $H$ at which the stability switch
occurs (left) and the mean zonal wind at the top level $\overline{u}_{1}$
(right) as a function of the eddy shape aspect ratio $j$, shown for several
values of the eddy diffusion $A$ in $\text{m}^{2}\text{ s}^{-1}$; $H=3.5\times
10^{-3}\text{KJ ton}^{-1}\text{ s}^{-1}$ (right); $\kappa=4\times
10^{-6}\text{ s}^{-1}$ (both).

The right hand side of figure 2 shows the mean zonal wind at the top level of
the eddy steady state as a function of the aspect ratio. Eddies close to a
circular shape, i.e. $j$ close to one, yield unrealistically high wind speeds
just above $215\text{ m }\text{s}^{-1}$. For more meridionally elongated
eddies the wind speed decreases for an aspect ratio of up to 5 and thereafter
increases again for very high eddy diffusion values. For $j\geq 2.5$ one gets
reasonably realistic values for the wind speed. However, before this upward
slope, the wind speed neither depends on the eddy diffusion nor on the surface
friction diffusion (the latter not shown).

For the net poleward heat transport the opposite holds: the longer the eddies
in the meridional direction, the higher the heat transport until it drops
again for very large $j$. Besides for these very high aspect ratios, the
surface friction diffusion $\kappa$ has again no impact on the values, whereas
a very high eddy diffusion $A$ yields a decreased heat transport for the whole
range of $j$.

As briefly mentioned earlier, the surface friction $\kappa$ could perhaps be
thought to act similarly to the eddy diffusion $A$, taking out energy from the
system and thus compensating for the heating $H$. However, the heating forces
a vertical zonal shear by thermal wind balance, but would not force a zonal
wind directly at the surface. The model can develop very low surface winds but
increase its shear by increasing upper level zonal winds. In this sense, the
surface friction does not directly compensate for the incoming heating $H$.

However, insofar that the heating induces surface winds, the surface friction
would represent a sink for those surface winds. In that case, one would expect
that a higher surface friction is able to dissipate the incoming energy and
therefore moves the stability switch to higher values of the heating. This
frictional adjustment can not be seen for realistic values of $\kappa$ (figure
3 on the left) but is observable for extremely high surface friction values
(figure 3 on the right).

Additionally, and in contrast to the eddy diffusion $A$, an extremely high
surface friction increases the mean zonal wind to unrealistically high values.
This behaviour of the mean flow was also observed by Novak, Ambaum, and Harvey
(2018) and Marshall et al. (2017) for the Antarctic Circumpolar Current.

Figure 3:  Contours of the mean zonal wind at the top level $\overline{u}_{1}$
for realistic (left) and extremely high (right) values of the surface friction
parameter $\kappa$; the dashed region is the parameter space where the zonal
steady state $P^{0}$ is attracting, the non-dashed region where the eddy
steady state $P^{*}$ is attracting; note that in the left figure the scale for
the surface friction is logarithmic; $A=10^{5}\text{ m}^{2}\text{ s}^{-1}$,
$j=2.5$.

##  5 Short and long wavelength cut-off

For sufficiently large eddy diffusion $A$, the left hand side of figure 2
shows a short (large $j)$ and long wavelength (small $j$) cut-off of the
instability of the zonal solution. In other words, our model only exhibits
finite-size eddies within a certain range of wavelengths and the zonal flow
stabilizes for very short and long waves.

A similar result follows from the classical stability analysis of the Phillips
model presented in, for example, Holton and Hakim (2013). This is done for a
simpler version of the two-level model, without surface friction, eddy
diffusion or external forcing due to heating, and is obtained by performing a
linear perturbation analysis. Assuming wave-like solutions for the barotropic
and baroclinic development of the streamfunction, they obtain a dispersion
relation for the phase speed of the waves and from this derive conditions for
instability of the zonal flow and the growth of perturbations. The zonal flow
is stable for very short waves because they are inefficient in converting
available into potential energy, and for very long waves, because they are
stabilized by the differential rotation of the plane along latitudes
($\beta$-effect). This instability is similar to the repelling regime of the
zonal steady state in our model and the results of Holton and Hakim (2013) are
compared to ours in the following.

The first condition in Holton and Hakim (2013) for instability to occur is

$\displaystyle\overline{u}_{\rm T}>\frac{\beta}{\lambda_{\rm R}^{2}}.$ (1n)

This is also Phillips’ criterion for baroclinic instability of a two-level
model, see Phillips (1954). Using the values in table C we obtain a minimal
shear of approximately $3.65\text{ m}\text{ s}^{-1}$, meaning that for a shear
below this value the zonal steady state must always be attracting. In fact, it
can be shown that for values of $\overline{u}_{\rm T}$ below this level our
model is always in the attracting region of the zonal steady state and as the
shear is increased this state becomes repelling.

According to Holton and Hakim (2013), the minimum value of $\overline{u}_{\rm
T}$ for which $P^{0}$ becomes repelling occurs when
$k^{2}=\sqrt{2}\lambda_{\rm R}^{2}/2$, where $k$ is the meridional wavenumber.
Using the relation $k=j\pi/W$, this condition becomes

$\displaystyle j=\frac{2^{-1/4}W\lambda_{\rm R}}{\pi}\approx 5.6,$ (1o)

where we used again the values in table C. Analysing our model yields exactly
that value of $j$ at which the mean zonal shear is minimal. The corresponding
minimum is $\overline{u}_{\rm T}\approx 4.1\text{ m}\text{ s}^{-1}$ and is
similar to the value given by Holton and Hakim (2013). This value of $j$ also
yields the wavenumber of maximal instability of the purely zonal flow.

The short wavelength cut-off observed for large eddy diffusion $A$ in figure 2
on the left is given in Holton and Hakim (2013) by the critical wavelength

$\displaystyle\rm L_{\rm c}=\frac{2\pi}{\lambda^{2}_{\rm R}}\approx 3\times
10^{3}\text{ km}.$ (1p)

For waves shorter than this value, instability can not occur and the zonal
steady state remains attracting. In our model this value corresponds to
$j>6.6$ and matches roughly the value for which the zonal steady state becomes
attracting again if the eddy diffusion is very high and the heating rate is
approximately the standard value given in table C.

##  6 Q factor

The quality factor Q is a dimensionless parameter that is used in mechanics
and electronics to describe the oscillatory behaviour of a damped system. In
general, the Q factor is defined as the ratio between the natural undamped
frequency of a dynamical system and its damping coefficient. Therefore, a high
Q factor indicates a larger number of oscillations before the system reaches
the attracting steady state, in contrast to a low Q factor indicating that the
damping is more dominant and the attracting steady state is reached after
fewer oscillations. A Q factor below $1/2$ corresponds to an overdamped system
that does not oscillate at all when displaced from its steady state, but
returns to it by exponential decay.

This concept can be applied to the theory of damped oscillators, which satisfy
the linear equation

$\displaystyle\ddot{x}+\lambda\dot{x}+\omega^{2}_{0}x=0,$ (1q)

where $\omega_{0}^{2}$ is the frequency of free oscillations and $\lambda$ the
damping coefficient. The general solution of equation (1q) is

$\displaystyle x=A\exp\left(-\frac{\lambda}{2}t\right)\exp\left(\pm
i\left(\omega^{2}_{0}-\frac{\lambda^{2}}{4}\right)^{1/2}t\right)=A\exp\left(pt\right),$
(1r)

where A is a fixed complex amplitude and $p=-\lambda/2\pm
i\left(\omega^{2}_{0}-\lambda^{2}/4\right)^{1/2}$. The generally accepted
definition of the Q factor for the damped oscillator is

$\displaystyle Q=\frac{\omega_{0}}{\lambda}.$ (1s)

We now consider a linearised system obeying the differential equation
$\text{d}x/\text{d}t=\alpha x$. This is associated with a linear damped
oscillator if we take $p=\alpha=\alpha_{\rm r}+i\alpha_{\rm i}$, where
$\alpha_{\rm r}<0$ is the real and $\alpha_{\rm i}$ the imaginary part of the
eigenvalue. With this identification, we find $\alpha_{\rm r}=-\lambda/2$ and
$\alpha_{\rm i}=\pm\left(\omega_{0}^{2}-\lambda^{2}/4\right)^{1/2}$ and hence
the Q factor of a linearised system can be expressed as

$\displaystyle Q=\frac{\left(\alpha_{\rm i}^{2}+\alpha_{\rm
r}^{2}\right)^{1/2}}{-2\alpha_{\rm r}}=\frac{|\alpha|}{-2\alpha_{\rm r}}.$
(1t)

By construction, $Q\geq 1/2$ for eigenvalues with a non-zero imaginary part,
which is associated with the oscillatory behaviour of the system. For real
eigenvalues, the Q factor is exactly $1/2$ and the system is said to be
critically damped. Similar to an over-damped system, the steady state is
approached without oscillations.

The left hand side of figure 4 shows contours of the Q factor in the heating
rate $H$ versus eddy diffusion $A$ plane. As expected, Q is always greater
than $1/2$ so the system oscillates everywhere in the considered parameter
space. Furthermore, it can be seen that an increase of the eddy diffusion
parameter $A$ leads to a decrease of the Q factor. This corresponds to the
fact that the diffusion acts as a damping and stabilizing factor of the system
and leads to a reduction of the available energy. Therefore, the higher the
diffusion, the less the initial transient relies on the eddy heat flux to
compensate the heat imbalance. In contrast to that, a larger heating rate $H$
increases the Q factor since more heat needs to be transported away from the
equator and therefore more oscillatory cycles of the eddies are needed during
the initial transient. The contours of Q in the heating rate versus surface
friction $\kappa$ plane are similar (not shown) and can be interpreted in the
same way.

On the right hand side of figure 4 the Q factor is shown as a function of the
eddy aspect ratio $j$ for different values of $A$. The highest Q factor occurs
for the smallest $j$, i.e. for the most symmetric eddy shape, and the larger
$j$, the smaller is the Q factor. Hence, for eddies that are more elongated in
the meridional direction, the system undergoes fewer oscillatory cycles before
it reaches the attracting steady state. In physical terms, meridionally
elongated eddies are more efficient in balancing the heat in the atmosphere
than eddies close to symmetric. This also coincides with the findings in the
earlier chapters. Apart from that, the dependence on the diffusion parameter
$A$ is the same as in the contour plot on the left side.

Figure 4:  LEFT: Contours of the Q factor in the heating rate $H$ versus eddy
diffusion $A$ plane, the dashed region is the parameter space where the zonal
steady state $P^{0}$ is attracting, the non-dashed region where the eddy
steady state $P^{*}$ is attracting, $\kappa=4\times 10^{-6}\text{ s}^{-1}$,
$j=2.5$; RIGHT: Q factor as a function of $j$ for various values of the eddy
diffusion $A$ in $\text{m}^{2}\text{ s}^{-1}$, $\kappa=4\times 10^{-6}\text{
s}^{-1}$, $H=3.5\times 10^{-3}\text{KJ ton}^{-1}\text{ s}^{-1}$.

##  7 Discussion and conclusion

We used Phillips’ classical two-level quasi-geostrophic model on the
$\beta$-plane (Phillips 1956) and performed a model reduction to develop a
system of six ordinary differential equations describing the interaction
between the mean flow and shear, and the eddy activity in the mid-latitude
atmosphere. This has been achieved by imposing a specific shape on the eddy
streamfunction and thus looking at the evolution of the amplitude of the
atmospheric modes considered in our spectral projection.

Instead of taking the usual approach of performing a normal mode instability
analysis, we used classical dynamical systems theory and determined two
physically possible steady states, namely a zonal solution, where the flow is
purely zonal and has no eddy contribution, and an eddy solution, where zonal
flow and eddies coexist. The considered parameter space is divided into two
separate regions where either the zonal or the eddy solution is attracting,
while the alternative solution is unstable and repels nearby initial
conditions.

The two competing states are in close correspondence with two possible
physical mechanisms that allow the atmosphere to reach a steady state in such
a way that the imbalance in the diabatic heating between equator to the poles
is compensated. Relatively weak imbalances in the diabatic heating can be
balanced by diffusion and no eddies need to develop to flatten the meridional
temperature gradient. Indeed, the flow is, in this case, baroclinically
stable. In this state the mean flow increases for stronger forcing by diabatic
heating, because the latter imposes a stronger temperature difference between
low and high latitudes. This monotonicity is lost as soon as the heating
reaches a certain threshold, where the eddies start growing to compensate the
incoming heat because the damping effect of the eddy diffusion is no longer
sufficient. The zonal flow is then baroclinically unstable so that eddies grow
and reach a finite steady-state amplitude, in such a way that the resulting
poleward eddy heat flux weakens the baroclinicity and therefore the mean flow.
This process is referred to as baroclinic adjustment (Stone 1978).

In this baroclinically unstable eddy steady state the zonal flow is
independent of the external forcing by the heating. In contrast, the net
poleward heat transport by the eddies and the eddy energy increase with
incoming heat, so incoming energy in form of diabatic forcing of the zonal
mean flow goes directly into the eddies in form of eddy kinetic energy, where
it is dissipated. Therefore, the stability switch from zonal to eddy steady
state describes a transfer of the location of dissipation of incoming energy.
In terms of the Lorenz energy cycle, for sufficiently large heating, the zonal
flow can no longer act as the reservoir for the incoming energy, so a new
reservoir needs to develop, namely the eddies. Thus, the reduced model is
exhibiting complete eddy saturation.

In contrast to Munday, Johnson, and Marshall (2013), stating that eddy-
resolving models of the ACC lead to a much lower sensitivity of the volume
transport to increased forcing, Marshall et al. (2017) suggested the
possibility to capture the physics of eddy saturation in models with
parameterised eddies. Thereafter, Mak et al. (2017) used a highly idealised
model configuration of the ACC with parameterised eddies to obtain a near
total eddy saturation. For the atmosphere, Novak, Ambaum, and Harvey (2018)
found a model where the mean baroclinicity is nearly independent of the
thermal forcing, whereas the eddy intensity responds more strongly to the
forcing and is independent of the eddy dissipation. In this context, the model
presented in this work is the first to exhibit complete eddy saturation in a
model of the mid-latitude atmosphere with parameterised eddies.

Another aspect of the eddy saturation mechanism described in Marshall et al.
(2017) is the increase of the volume transport with the bottom drag, and one
would expect a similar effect of the surface friction in our model,
dissipating energy from the system and therefore compensating for the incoming
heat. This would lead to a stability switch at higher values of the heating
and to an increase of the zonal wind speed. However, this effect can only be
seen for unrealistically high values of the surface friction, perhaps
indicating a limitation of our model. Due to the two-level setup the model is
able to develop a high wind shear induced by the incoming heating by
increasing the upper level zonal wind but maintaining a very low surface wind.
Hence, the surface friction can not directly compensate for the incoming
heating in the model, although in the atmosphere the heating actually induces
stronger surface winds.

In addition to the zonal wind being independent of the incoming heating and
the surface friction (for realistic values) in the baroclinically unstable
eddy steady state, the wind speed is independent of the eddy diffusion as
well. Instead, an increased diffusion leads to a decrease of eddy energy,
indicating the damping effect of eddy dissipation which is known as
dissipative control (e.g. Novak, Ambaum, and Harvey (2018)), and is in
accordance with Marshall et al. (2017), stating that the equilibrium volume
transport is controlled by the ACC requiring sufficiently strong vertical
shear induced by external forcing to overcome the stabilising role of the eddy
dissipation.

Despite the simplicity of the model and the limiting assumptions made on the
shape of the eddies, our model exhibits complete eddy saturation, a regime of
the mid-latitude atmosphere that appears to be supported by numerical models.
The eddy saturation mechanism is explained as a switch of stability from a
zonal mean state to a mean state with finite eddy amplitude. This work is of
relevance to climate modeling, where the response of the mid-latitude storm
tracks to changes in the diabatic heating is still not fully understood.

## Acknowledgement

Melanie Kobras is supported by the U.K. Engineering and Physical Sciences
Research Council (Grant EP/R513301/1). Valerio Lucarini acknowledges the
support received from the EPSRC project EP/T018178/1 and from the EU Horizon
2020 project TiPES (grant no. 820970).

## References

  * Ambaum and Novak (2014) Ambaum, Maarten H. P., and Lenka Novak. 2014. “A nonlinear oscillator describing storm track variability.” _Quarterly Journal of the Royal Meteorological Society_ 140 (685): 2680–2684.
  * Guckenheimer and Holmes (1983) Guckenheimer, John, and Philip Holmes. 1983. _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_. 1st ed., Applied Mathematical Sciences. Springer, New York, NY.
  * Holton and Hakim (2013) Holton, James R., and Gregory J. Hakim. 2013. _An Introduction to Dynamic Meteorology_. 5th ed. Academic Press.
  * Hoskins and James (2014) Hoskins, Brian J., and Ian M. James. 2014. _Fluid Dynamics of the Midlatitude Atmosphere_. 1st ed. John Wiley & Sons, Ltd.
  * Lorenz (1967) Lorenz, Edward N. 1967. _The nature and theory of the general circulation of the atmosphere_. Vol. 218. World Meteorological Organization Geneva.
  * Lucarini et al. (2014) Lucarini, Valerio, Richard Blender, Corentin Herbert, Francesco Ragone, Salvatore Pascale, and Jeroen Wouters. 2014. “Mathematical and physical ideas for climate science.” _Reviews of Geophysics_ 52 (4): 809–859.
  * Lucarini, Speranza, and Vitolo (2007) Lucarini, Valerio, Antonio Speranza, and Renato Vitolo. 2007. “Parametric smoothness and self-scaling of the statistical properties of a minimal climate model: What beyond the mean field theories?” _Physica D: Nonlinear Phenomena_ 234 (2): 105–123.
  * Mak et al. (2017) Mak, J., D.P. Marshall, J.R. Maddison, and S.D. Bachman. 2017. “Emergent eddy saturation from an energy constrained eddy parameterisation.” _Ocean Modelling_ 112: 125–138.
  * Marshall et al. (2017) Marshall, David P., Maarten H. P. Ambaum, James R. Maddison, David R. Munday, and Lenka Novak. 2017. “Eddy saturation and frictional control of the Antarctic Circumpolar Current.” _Geophysical Research Letters_ 44 (1): 286–292.
  * Munday, Johnson, and Marshall (2013) Munday, David R., Helen L. Johnson, and David P. Marshall. 2013. “Eddy Saturation of Equilibrated Circumpolar Currents.” _Journal of Physical Oceanography_ 43 (3): 507–532.
  * Munk and Palmén (1951) Munk, W. H., and E. Palmén. 1951. “Note on the Dynamics of the Antarctic Circumpolar Current.” _Tellus_ 3 (1): 53–55.
  * Nadeau and Ferrari (2015) Nadeau, Louis-Philippe, and Raffaele Ferrari. 2015. “The Role of Closed Gyres in Setting the Zonal Transport of the Antarctic Circumpolar Current.” _Journal of Physical Oceanography_ 45 (6): 1491–1509.
  * Novak, Ambaum, and Harvey (2018) Novak, Lenka, Maarten H. P. Ambaum, and Ben J. Harvey. 2018. “Baroclinic Adjustment and Dissipative Control of Storm Tracks.” _Journal of the Atmospheric Sciences_ 75 (9): 2955–2970.
  * Pedlosky (1979) Pedlosky, Joseph. 1979. _Geophysical Fluid Dynamics_. Springer Study Edition. Springer-Verlag Berlin Heidelberg.
  * Peixoto and Oort (1992) Peixoto, Jose P., and Abraham H. Oort. 1992. _Physics of Climate_. American Institute of Physics.
  * Phillips (1954) Phillips, Norman A. 1954. “Energy Transformations and Meridional Circulations associated with simple Baroclinic Waves in a two-level, Quasi-geostrophic Model.” _Tellus_ 6 (3): 274–286.
  * Phillips (1956) Phillips, Norman A. 1956. “The general circulation of the atmosphere: A numerical experiment.” _Quarterly Journal of the Royal Meteorological Society_ 82 (352): 123–164.
  * Stone (1978) Stone, Peter H. 1978. “Baroclinic Adjustment.” _Journal of Atmospheric Sciences_ 35 (4): 561–571.
  * Straub (1993) Straub, David N. 1993. “On the Transport and Angular Momentum Balance of Channel Models of the Antarctic Circumpolar Current.” _Journal of Physical Oceanography_ 23 (4): 776–782.
  * Thompson (1987) Thompson, Philip D. 1987. “Large-Scale Dynamical Response to Differential Heating: Statistical Equilibrium States and Amplitude Vacillation.” _Journal of the Atmospheric Sciences_ 44 (8): 1237–1248.

##  Appendix A Additional evolution equations

The evolution equations for the mean meridional kinetic energy
$\overline{v^{\prime 2}_{\rm m}}$, the temperature variance
$\overline{v^{\prime 2}_{\rm T}}$ and the cross-correlation between
temperature and geopotential $\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm T}}$
are

$\displaystyle\frac{\upartial\overline{v^{\prime 2}_{\rm m}}}{\upartial t}$
$\displaystyle=\frac{2k^{2}(\lambda^{2}_{\nabla}-\lambda^{2}_{y})}{\lambda^{2}_{\nabla}}\overline{u}_{\rm
T}\overline{v^{\prime}_{\rm m}\psi^{\prime}_{\rm
T}}-(2A\lambda^{2}_{\nabla}+\kappa)\overline{v^{\prime 2}_{\rm
m}}+2\kappa\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm T}},$ (1ua)
$\displaystyle\frac{\upartial\overline{v^{\prime 2}_{\rm T}}}{\upartial t}$
$\displaystyle=\frac{2k^{2}(\lambda^{2}_{y}+\lambda^{2}_{\rm
R}-\lambda^{2}_{\nabla})}{\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R}}\overline{u}_{\rm T}\overline{v^{\prime}_{\rm m}\psi^{\prime}_{\rm
T}}-\left(2A\lambda^{2}_{\nabla}+\frac{2\kappa\lambda^{2}_{\nabla}}{\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R}}\right)\overline{v^{\prime 2}_{\rm
T}}+\frac{\kappa\lambda^{2}_{\nabla}}{\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R}}\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm T}},$ (1ub)
$\displaystyle\frac{\upartial}{\upartial t}\overline{v^{\prime}_{\rm
m}v^{\prime}_{\rm T}}$ $\displaystyle=\frac{\lambda^{2}_{\rm
R}k^{2}}{\lambda^{2}_{\nabla}(\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}\left(\beta+\lambda^{2}_{y}\overline{u}_{\rm
m}\right)\overline{v^{\prime}_{\rm m}\psi^{\prime}_{\rm
T}}-\left(2A\lambda^{2}_{\nabla}+\frac{\kappa(3\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}{2(\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}\right)\overline{v^{\prime}_{\rm m}v^{\prime}_{\rm
T}}+\kappa\overline{v^{\prime 2}_{\rm T}}$
$\displaystyle\qquad+\frac{\kappa\lambda^{2}_{\nabla}}{2(\lambda^{2}_{\nabla}+\lambda^{2}_{\rm
R})}\overline{v^{\prime 2}_{\rm m}}.$ (1uc)

##  Appendix B Dimensionless constants of model (1ja)-(1jf)

$\displaystyle\alpha=\frac{\lambda_{\rm R}}{\beta}\kappa$
$\displaystyle,\quad\gamma=\frac{\lambda^{2}_{y}}{\lambda^{2}_{y}+\lambda^{2}_{\rm
R}},\quad\delta=\frac{\lambda^{2}_{\rm
R}}{\hat{j}(\hat{j}\lambda^{2}_{y}+\lambda^{2}_{\rm
R})},\quad\epsilon=\frac{3\hat{j}\lambda^{2}_{y}+\lambda^{2}_{\rm
R}}{\hat{j}\lambda^{2}_{y}+\lambda^{2}_{\rm R}},\quad\nu=\frac{\lambda_{\rm
R}\lambda^{2}_{y}A}{\beta},$ $\displaystyle\zeta$
$\displaystyle=\frac{\hat{j}\lambda^{2}_{y}}{\hat{j}\lambda^{2}_{y}+\lambda^{2}_{\rm
R}},\quad\mu=\frac{\lambda^{2}_{\rm
R}+(1-\hat{j})\lambda^{2}_{y}}{\hat{j}\lambda^{2}_{y}+\lambda^{2}_{\rm
R}}\quad\eta=\frac{\lambda^{5}_{\rm R}R}{\beta^{2}f_{0}c_{\rm
p}W(\lambda^{2}_{y}+\lambda^{2}_{\rm R})},$ (1v)

where

$\displaystyle\hat{j}=\frac{j^{2}+1}{2},\ a=\nu+\frac{\alpha}{2},\
b=\nu+\gamma\alpha,\ c=2\hat{j}\nu+\frac{\alpha\epsilon}{2},\
d=2\hat{j}\nu+\alpha,\ e=2\hat{j}\nu+2\alpha\zeta.$ (1w)

##  Appendix C Parameter space

Fixed parameter values for stability analysis.  Parameter $\lambda^{2}_{\rm
R}$ $\lambda^{2}_{y}$ $\beta$ $R$ $f_{0}$ $c_{\rm p}$ $W$ Value $4.39\times
10^{-12}$ $1.97\times 10^{-13}$ $1.6\times 10^{-11}$ $287$ $10^{-4}$ $1004$
$10^{7}$ Unit $\text{m}^{-2}$ $\text{m}^{-2}$ $\text{m}^{-1}\text{ s}^{-1}$
$\text{J K}^{-1}\text{ kg}^{-1}$ $\text{s}^{-1}$ $\text{J K}^{-1}\text{
kg}^{-1}$ m

Varying parameter values for stability analysis.  Parameter $\kappa$ $A$ $H$
$j$ Minimum $10^{-7}$ $1$ $0$ $0.1$ Standard $4\times 10^{-6}$ $10^{5}$
$3.5\times 10^{-3}$ - Maximum $8\times 10^{-6}$ $10^{5}$ $4\times 10^{-3}$
$6.8$ Unit $\text{s}^{-1}$ $\text{m}^{2}\text{ s}^{-1}$ $\text{KJ
ton}^{-1}\text{ s}^{-1}$ -

##  Appendix D Other variables of eddy steady state

Non-dimensional eddy steady state variables of the mean meridional kinetic
energy $K^{*}$, the temperature variance $V^{*}$ and the cross-correlation
between temperature and geopotential $X^{*}$:

$\displaystyle K^{*}$
$\displaystyle=\frac{2^{2}j^{2}(j^{2}-1)\lambda^{2}_{y}}{(j^{2}+1)d\lambda^{2}_{\rm
R}}S^{*}T^{*}+\frac{2\alpha}{d}X^{*},$ $\displaystyle
V^{*}=\frac{2^{2}j^{2}\mu\lambda^{2}_{y}}{e\lambda^{2}_{\rm
R}}S^{*}T^{*}+\frac{\alpha\zeta}{e}X^{*},$ $\displaystyle X^{*}$
$\displaystyle=\left(\frac{2j^{2}\alpha\lambda^{2}_{y}g}{a\lambda^{2}_{\rm
R}}S^{*}+j^{2}\delta\right)\left(c-\alpha\zeta\frac{d+e}{de}\right)^{-1}T^{*},$
(1af)

with constants

$\displaystyle B$
$\displaystyle=2\alpha\delta(g-f)\left(c-\alpha^{2}\zeta\frac{d+e}{de}\right)^{-1},\quad
C=\frac{\delta^{2}\lambda^{2}_{\rm
R}}{2\lambda^{2}_{y}}\left(c-\alpha^{2}\zeta\frac{d+e}{de}\right)^{-1}+c,$
$\displaystyle D$ $\displaystyle=8\frac{\lambda^{2}_{y}}{\lambda^{2}_{\rm
R}}\left(2a^{2}\mu\left(1-\frac{1}{\hat{j}}\right)\frac{d+e}{de}+\alpha^{2}fg\left(c-\alpha^{2}\zeta\frac{d+e}{de}\right)^{-1}\right),$
$\displaystyle f$
$\displaystyle=\frac{2a\mu}{d}+(1-\frac{1}{\hat{j}})\frac{a\zeta}{e}-\delta,\quad
g=\delta+\frac{(j^{2}-1)a\zeta}{(j^{2}+1)d}+\frac{2a\mu}{e}.$

Here, $B$, $C$ and $g$ are positive for the whole parameter space given by
tables C and C, $D$ and $f$ are non-negative for $j$ larger than 1 and for the
positivity of $D$ additionally $\kappa\geq 3\times 10^{-6}\text{ s}^{-1}$ is
required.


11institutetext: Centre for Digital Music, Queen Mary University of London, UK
11email<EMAIL_ADDRESS>
# Zero-shot Singing Technique Conversion

Brendan O’Connor 11 Simon Dixon 11 George Fazekas This research is funded by
the EPSRC and AHRC Centre for Doctoral Training in Media and Arts Technology
(EP/L01632X/1).11

###### Abstract

In this paper we propose modifications to the neural network framework, AutoVC
[17] for the task of singing technique conversion. This includes utilising a
pretrained singing technique encoder which extracts technique information,
upon which a decoder is conditioned during training. By swapping out a source
singer’s technique information for that of the target’s during conversion, the
input spectrogram is reconstructed with the target’s technique. We document
the beneficial effects of omitting the latent loss, the importance of
sequential training, and our process for fine-tuning the bottleneck. We also
conducted a listening study where participants rate the specificity of
technique-converted voices as well as their naturalness. From this we are able
to conclude how effective the technique conversions are and how different
conditions affect them, while assessing the model’s ability to reconstruct its
input data.

###### Keywords:

Voice synthesis, singing synthesis, style transfer, neural network, singing
technique, timbre conversion, conditional autoencoder, sequential training,
latent loss

##  1 Introduction

Voice conversion (VC) is the task of converting the timbre of the voice so
that the linguistic content is perceived to be spoken by a different person.
It has been explored in relation to both singing and speech, which both
possess different attributes consideration. Singing voice analysis is
considerably more focused on sustained notes, harmonic/rhythmic structure, and
relative pitch. In speech, these musical values are non-existent. Instead
there is greater emphasis on aperiodic aspects, such as consonant utterances
and rapidly shifting spectral envelopes. Tasks like VC and text-to-speech are
in far more demand in the industry than singing-related tasks, and have
therefore monopolised the spotlight in voice analysis and synthesis research.
The latest approaches towards VC achieving state-of-the-art conversions
utilise probabilistic machine learning techniques. Public domain speech
datasets also vastly overshadow singing datasets in size and availability
[11], and so there is still much to be explored in relation to singing
analysis and synthesis.

In this paper we tackle the task of singing technique conversion (STC) - the
task of converting a singing technique without affecting the perceived
identity of the singer, musical structure or linguistic content. We define
singing technique as the method of voice production to achieve different
timbres by adjusting the airflow, vocal folds, vocal tract shape, and
sympathetic vibrations in the body [5]. We regard STC as a variation of voice
conversion (VC), where the possibilities of voice transformation are
restricted to be within a realistic variance of timbre for any given singer.
We chose the term singing techniques as opposed to singing style, due to the
latter term’s inconsistent use in literature, often referring to a range of
very different audio and musical attributes due to its lack of reference to a
concrete audio or singing concept.

To achieve STC, we apply a neural network model in the form of the conditioned
autoencoder, AutoVC [17]. We discuss certain adaptations made to the
architecture and investigate the effects of training it on different
permutations of several datasets. To evaluate the model’s ability to perform
STC, we had participants rate the naturalness of the voice and guess what the
target singing technique was supposed to be. Examples of audio used in this
listening test can be found
online.111https://github.com/Trebolium/singing_technique_conversion

Real-time pitch correction algorithms have become commonplace in the music
industry and influence the characteristics of modern pop singers today. We
believe that the refined task of STC could have a similar influence on music
production as it opens up the possibility of artistically manipulating a
singer’s performance, rather than just quantising their pitch. Over the last 5
years, many machine-learning approaches have been proposed to tackle voice
transformation for speech (as discussed in the next section), but much less
attention has been given to transforming the expression of the singing voice.

##  2 Related Work

Recent research in VC has been based on neural networks, which have influenced
the frameworks proposed in this paper. [15] conditioned an autoencoder
(trained on linguistic data) on speaker embeddings generated from a separately
trained classifier network. During inference, these embeddings could be
replaced to achieve VC. AutoVC [17] adapted this method to work with
spectrograms, which will be described in detail in Section 4. This was
improved upon by conditioning the network on pitch contours to enforce prosody
during conversion [18], and further disentanglement was achieved for timbre,
pitch contours, rhythm and utterances simultaneously by utilising 3 separate
bottlenecks with different restrictions [19]. [26] achieve VC by using vector
quantisation to separate speaker and content information, and later utilised
U-nets [21] to compensate for information lost during vector quantisation.

The application of the variational autoencoder (VAE) is well suited for ‘many-
to-many’ conversions (where all examples used for inference are seen during
training). [16] use fully convolutional VAEs, conditioned on acoustic
features, to perform VC. They combine spectral features of both converted and
unconverted reconstructed audio in order to avoid over-smoothing - a known
issue with VAEs. While VAEs present an elegant framework, they produce
‘blurry’ results. Generative Adversarial Networks (GANs) have been known to
reproduce better quality reconstructions of images than VAEs. However as they
come without an autoencoder they are harder to train and suffer from ‘mode
collapse’, and there has not yet been an elegant proposal for combining VAEs
with GANs [22]. The use of VAEs has the added benefit of utilising
unsupervised learning, which bypasses the issue of low resources regarding
labelled singing datasets. [6] used Gaussian-Mixture VAEs (GMVAEs) for
controllable speech synthesis, modelling the different attributes of speech as
separate prior distributions before combining them in a VAE. For singing voice
conversion, [13] adapted AutoVC by conditioning the network on pitch contours
transposed to a suitable register for the converted singing, achievable
through the implementation of a vocoder. [8] utilised a Wasserstein-GAN
framework, using a decoder for pitch contours and another for generating
‘formant masks’. The product of these two decoders is the estimated mel-
spectrogram for singing. They later explored the capabilities of this
framework to achieve timbre and singing style disentanglement [9], where a
singing query is converted into a singer identity embedding and used to
condition both the pitch skeleton and formant-mask encoders on pitch
modulation style and singer timbre, respectively. [10] present the only other
research we know of that addresses STC. They use GMVAEs to model singer and
technique information to perform many-to-many conversions using a VAE
architecture that utilises a convolutional recurrent neural network (CRNN)
architecture.

The issue remains however, of what can be done with singing datasets which are
small and few. [13] notes that the generalisation of the AutoVC framework
allows it to be utilised as a Universal Background Model. [1] synthesise
monophonic singing datasets by superimposing pitch contours on existing speech
datasets. [3] use several autoencoder instances, trained separately on vocoder
spectral data and music mixtures, while being conditioned on shared content
embeddings and 1-hot speaker embeddings to produce a final network that is
singer-independent and generates monophonic singing from musical mixtures.
[12] generate novel speaker embeddings by combining embeddings from existing
singers as a method of data augmentation.

##  3 Architecture

We use the AutoVC framework [17] for singing technique transformation, due to
its elegant method of applying disentanglement. It is also capable of
converting between source and target examples that have not been seen in the
training datasets (zero-shot conversion). In AutoVC, a standard autoencoder
architecture is conditioned on speaker embeddings that uniquely describe the
timbre of a speaker to perform VC on spectrograms. These embeddings are
generated by a pretrained speaker verification network [24]. The spectrograms
are concatenated with these speaker embeddings, and fed through an encoder
$E_{c}$, after which the encoded information is again concatenated with
speaker embeddings before being fed to the decoder $D_{c}$. This conditioning,
combined with careful calibration of an appropriate bottleneck size, allows
the autoencoder to disentangle speaker timbre from utterance information.
AutoVC also contains a ‘postnet’ convolutional layer which is appended to the
decoder to further develop a refined spectrogram from the decoder’s output.
After training, the speaker embeddings concatenated at the bottleneck can be
swapped out to achieve VC. The loss function for AutoVC is a weighted
combination of the self-reconstruction loss for both the decoder
($L_{\mathit{decoder}}$) and the postnet ($L_{\mathit{postnet}}$) output
spectrograms, and the latent loss ($L_{\mathit{latent}}$). The latent loss
represents the difference between the bottleneck’s embedding $E_{c}(x)$ for
the input x and its reconstructed form $E_{c}(\hat{x})$. This is summarised in
Equation 1, where $\mu$ and $\lambda$ are empirically determined weights.
Further details of AutoVC’s architecture are given by [17], which we follow in
our implementation except for several adjustments discussed in this section.

$L_{\mathit{total}}=L_{\mathit{decoder}}+\mu L_{\mathit{postnet}}+\lambda
L_{\mathit{latent}}.$ (1)

We will herein refer to our implementation of AutoVC as AutoSTC to reflect its
purpose of STC. To facilitate this, we developed our own singing technique
encoder (STE) to replace the external speaker encoder that was used in the
original implementation. The STE is initially trained as a classifier. It
takes a mel spectrogram as input, which is split into chunks of 0.5 seconds.
These are fed in parallel through a neural network consisting of four
2D-convolutional layers (each of which is followed by batch normalisation,
ReLU activation and max-pooling), two dense layers, two BLSTMs, a simplified
attention mechanism [20], two more dense layers and finally a classification
layer. This architecture was adapted from the VAE used by [10] and influenced
by [4]. This network is able to achieve 86% accuracy when classifying singing
techniques within a test set of VocalSet (detailed in Section 4, while our
implementation of a 1D convolutional network on the waveform data as described
by [25] only scored 57%. During conversion, the STE’s embedding preceding the
classification layer are used for concatenation and conditioning with AutoVC
as described above in place of the external speaker encoder embeddings.

##  4 Training and Inference

The Vocalset dataset [25] used to train the STE consists of recordings of 20
singers performing several musical exercises with different singing
techniques. We chose a subset containing the techniques belt, straight,
vibrato, lip trill, vocal fry and breathy, trimming off excess files that
appear in one class but not the other, to yield a balanced class subset of
1182 examples (roughly 8K seconds). As the dataset is so small we only
partition it into training and test sets by 8:2.

As [13] showed that the sequence of training of different datasets is
important, AutoSTC was trained using subsets taken from VocalSet, VCTK [23]
and the raw singer recordings from MedleyDB [2] in various permutations. All
data was sampled at 16kHz and transformed into 80-bin mel spectrograms. While
being trained on one dataset, AutoSTC was simultaneously tested on test sets
from all three datasets in between training iterations (the VocalSet test set
was the same set omitted when training the STE). We recorded the number
iterations and loss values for each dataset where the loss showed no further
improvement into Table 1, and transferred the saved neural network parameters
of a nearby checkpoint to the proceeding dataset training session in the
sequence. We trained AutoSTC once for every permutation of the datasets. Table
1 shows that the order in which datasets are fed to the network does have a
considerable impact on its loss. The paths Vc->Vs->Md (spanning 750k training
steps) and Vc->Md->Vs (500k steps) led to the lowest loss values for MedleyDb
and VocalSet reconstruction respectively, and were used to train models that
generated the examples used in our listening test (see Section 5).

Table 1: Shows losses and training iterations (in parentheses) for VocalSet
(left) and MedleyDB (right) along alternative paths. The optimum training path
is highlighted in bold from left to right. Training that leads to an increase
in loss is indicated with a circumflex, at which point that path is abandoned.
For space, the dataset names are shortened as follows: VCTK:Vc, VocalSet:Vs,
MedleyDB:Md.

Loss-Iteration for Vs
---
|  | Vs | 0.0274(100k) | Md | ^
Vc | 0.0653(300k) | Md | 0.0386(150k) | Vs | 0.0268(50k)
|  | Vc | ^ | Md | -
Vs | 0.0347(150k) | Md | ^ | Vs | -
|  | Vs | 0.0290(50k) | Vc | ^
Md | 0.0500(200k) | Vc | ^ | Vc | -
Loss-Iteration for Md
---
|  | Vs | 0.0474(150k) | Md | 0.0265(100k)
Vc | 0.0479(500k) | Md | 0.0295(150k) | Vs | ^
|  | Vc | 0.0474(100k) | Md | 0.0301(50k)
Vs | 0.0562(150k) | Md | 0.0370(100k) | Vs | ^
|  | Vs | ^ | Vc | -
Md | 0.0367(150k) | Vc | ^ | Vc | -

We found the L1 loss between decoder/postnet self-reconstructions and the
input encouraged better convergence over L2 loss. We also tested the impact of
excluding the latent loss for 100K steps for both VC and STC tasks. Results
showed that training without latent loss performs significantly better for
both tasks. The loss for STC with latent loss was 0.0237 (and 0.0185 without
loss), while spectrograms were blurry and the audio lacked microtonal
variation or vibrato, leaving a ‘bubbliness’ artefact in its absence. Vowels
were also not consistently reproduced. These shortcomings however are less
noticeable for speech than singing. This result is worth highlighting, as
latent loss has been used consistently in frameworks of a similar nature [13,
15, 17].

Further preliminary trials allowed us to fine-tune AutoSTC’s bottleneck. We
analysed the resynthesised audio and noted that the net size of the feature
space was more indicative of audio quality than focusing on time/frequency
axis fine-tuning separately. We estimated the threshold to be a downsampling
factor of 16 for the time-axis, with each timestep containing 32 features.
Lower dimensionality representations resulted in deterioration in the
reconstructed audio in a very similar manner to that of the network with
latent loss included.

##  5 Experiment Design

To evaluate our proposed network’s ability to perform STC, we conducted a
listening study, where 19 participants evaluated the converted audio for
specificity and naturalness under the different conditions of models, gender,
and source/target techniques used. The first of the models (Vs1) was trained
on VocalSet alone and converted Vocalset data. The second (Vs2) was trained
using the optimum path for Vocalset presented in Table 1 to also convert
VocalSet data, while the third (M1) used the optimum path for MedleyDB to
convert MedleyDB data. Converted spectrograms were resynthesised using a
pretrained wavenet model provided by the author of the AutoVC
paper222https://github.com/auspicious3000/autovc. Each model produced 8 random
examples per participant, while adhering to a balanced representation of both
gender and subset (train/test) conditions. To evaluate naturalness, we asked
participants to consider how synthetic/realistic the voice itself sounded, and
rate them on a scale of 1 (very unnatural) to 5 (very natural). In a separate
task to evaluate specificity, participants were given a reference recording
featuring a converted singing technique along with 6 unlabelled candidate
recordings from the same singer to choose from. These candidate recordings
were randomly selected from the relevant dataset partition, so that they each
featured 1 of the 6 potential target techniques assigned to the reference
recording. Participants were asked to select one recording they thought
featured a singing technique closest to that of the reference recording, or
more if the answer seemed ambiguous. In the case where reference recordings
were converted MedleyDB examples, no ground truth labels existed, and so a
singer of the same gender was randomly chosen from Vocalset to represent the 6
candidate singing techniques instead. Each of these tasks was presented 24
times. 6 resynthesised recordings of unconverted audio were also evaluated for
naturalness. The interface was built using the Web Audio Evaluation Tool [7].

##  6 Results

The Mean Opinion Score (MOS) for unconverted data was 3.75 ± 0.34, and is
important to consider when analysing the results of the study. This highlights
the fact that a considerable amount of perceived naturalness has already been
lost during the wavenet resynthesis process, and that the MOS values for
technique conversion should be considered with this in mind. It is the
comparison between conditions that we are interested in.

To calculate the similarity score $S$ for each condition, we used the formula
in Equation 2, where $P_{n}$ is a binary vector reflecting a participant’s
true/false predictions (identifying whether each candidate technique was the
same as what was presented in the reference audio) for the $n$th task, $C_{n}$
is a 1-hot vector reflecting the correct technique for the task, and $N$ is
the total number of tasks in the given condition. The similarity score is an
average count of correct predictions weighted by the reciprocal of the number
of predictions made for the corresponding task.

$S=\frac{1}{N}\sum_{n=1}^{N}\frac{P_{n}.C_{n}}{||P_{n}||_{1}}.$ (2)

Figure 1 displays the results obtained from the listening study. The top graph
displays MOS values for naturalness, with whiskers indicating the confidence
intervals. The lower graph displays similarity scores. The combination of
these two graphs give us insight into how each of our models perform, and what
conditions influence the naturalness and specificity of the converted singing.

Figure 1: Top: Bar graph showing naturalness (MOS values and confidence
intervals) for all conditions. The colours group together the conditions for
(left to right): models, subsets, genders, source technique and target
technique. Bottom: Bar graph showing similarity scores determined by Equation
2, relative to correct answers.

We detected from a Spearman’s rank analysis that MOS and similarity scores
were not significantly correlated. Similarity scores across all conditions
measure higher than the chance level (0.16), which suggests that our models
have some success in converting to the target techniques. The condition of
source-technique groups does not significantly influence recognisability of
the converted singing technique. However, the data would suggest that the
features of the target techniques trill and breathy are significantly more
distinguishable than the rest. Vibrato scored the lowest for similarity,
suggesting that this was a particularly difficult technique to synthesise
convincingly. The reason for this is most likely due to the fact that
VocalSet, upon which the STE network was trained, contains numerous examples
labelled as belt and straight, while still featuring a considerable amount of
frequency modulation (a unique feature of vibrato), making it difficult for
AutoSTC to disentangle vibrato from other techniques effectively. It may also
be the case that AutoSTC has difficulty disentangling vibrato from pitch
contours. Alternatively it is possible that our models instead focused on
altering the phonation modes associated with vibrato, which would be
considerably less obvious to listeners than identifying whether frequency
modulation is occurring.

The inclusion of all datasets in training Vs2 seemingly diminished its ability
to accurately convert techniques (although the difference was not
statistically significant). The M1 model scored significantly worse than the
other models, which tells us that the features learned to generate technique
embeddings from the STE network were not generalisable to data outside the
dataset the VTE was trained on. There was also no statistically significant
difference between gender and subset similarity scores.

In regards to MOS results, the target technique trill scored lowest,
suggesting that conversions to a trill technique may sound unnatural. Vs2
samples were significantly higher than Vs1 and M1, which suggests that
providing the network with multiple datasets does improve its ability to
synthesise natural sounding data. The target technique condition vibrato
scored the highest, but as mentioned above, this may be because the network is
making changes more subtle than the frequency modulation which lessens the
amount of transformation required, causing less synthetic artefacts. It is
also perfectly possible that participants simply perceive the singing voice to
be more natural when vibrato is present.

##  7 Conclusion

In this paper we have presented a network for vocal technique classification,
and the first network to perform zero-shot conversion on singing techniques,
achieving above chance level for all tested conditions. We have demonstrated
that omitting latent loss and choosing the order in which AutoSTC was fed
different datasets significantly diminished its reconstruction loss, improving
its ability to reconstruct mel spectrograms. However we can conclude from the
results of the listening study that this does not have any significant effect
on AutoSTC’s ability to perform technique conversion and may even diminish it.
We therefore conclude that the features generated by supervised learning on
the labelled VocalSet dataset are not sufficient to generalise to recordings
of other singers. We also consider that the appearance of frequency modulation
in other techniques in VocalSet may have forced the network to give less
importance to this vibrato feature (we have however witnessed conversions
where frequency modulation was synthesised, but in very limited cases, so we
can not rule out the possibility that the AutoSTC framework is incapable of
converting singing technique features beyond their spectral filter
properties). The findings of our listening study are in agreement the vocal
timbre maps generated in our previous research [14].

Augmentation techniques such as those discussed in Section 2 may improve the
generalisation of the VTE to unseen data. We would also like to apply the
Generalised End-to-End Loss techniques from [24] to the VTE and fine-tune its
output embedding size. Due to shortcomings in labelled datasets, we will
explore unsupervised/semi-supervised networks such as VAEs. It may also be
worth investigating how AutoSTC performs when we condition it on further
attributes such as speaker identity, pitch contours and vowel sounds. As we
consider STC to be a restricted variation of VC and the fact that there are
considerably larger datasets for speech, it may also be worth exploring the
effects of pre-training an AutoVC framework for VC before switching its
speaker encoder for the singing technique encoder and training it for STC. In
future work we will also consider alternative options to the speech-trained
wavenet vocoder as this has introduced artefacts to the audio that likely
lowered MOS ratings for all audio. We have also observed that AutoSTC was
unintentionally able to remove vibrato from singing when underfitting, which
may be a capability worth fine-tuning in future work.

## References

  * [1] Basak, S., Agarwal, S., Ganapathy, S., Takahashi, N.: End-to-end Lyrics Recognition with Voice to Singing Style Transfer. In: 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 266–270 (2021) https://doi.org/10.1109/ICASSP39728.2021.9415096
  * [2] Bittner, R., Salamon, J., Tierney, M., Mauch, M., Cannam, C., Bello, J.: MedleyDB: A Multitrack Dataset for Annotation-intensive MIR Research. In: 15th International Society for Music Information Retrieval Conference, pp. 155–160. (2014)
  * [3] Chandna, P., Blaauw, M., Bonada, J., Gomez, E: Content Based Singing Voice Extraction from a Musical Mixture. In: 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 781-785 (2020)
  * [4] Choi, K., Fazekas, G., Sandler, M., Cho, K.: Convolutional Recurrentneural Networks for Music Classification. In: 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 2392–2396. (2017) https://doi.org/10.1109/ICASSP.2017.7952585
  * [5] Heidemann, K.: A System for Describing Vocal Timbre in Popular Pong. In: Music Theory Online, vol. 22, (2016) https://doi.org/10.30535/mto.22.1.2
  * [6] Hsu, W.-N., Zhang, Y., Weiss, R. J., Zen, H., Wu, Y., Wang, Y., Pang, R.: Hierarchical Generative Modeling for Controllable Speech Synthesis. In: The International Conference on Learning Representations, p. 27. New Orleans, LA, USA (2019)
  * [7] Jillings, N., Moffat, D., De Man, B., Reiss, J. D.: Web Audio Evaluation Tool: A browse-based listening environment. In: 12th Sound and Music Computing Conference. Maynooth, Ireland. (2015)
  * [8] Lee, J., Choi, H.-S., Jeon, C.-B., Koo, J., Lee, K.: Adversarially Trained End-to-end Korean Singing Voice Synthesis System. arXic preprint arXiv:1908.01919 (2019)
  * [9] Lee, J., Choi, H.-S., Koo, J., Lee, K.: Disentangling Timbre and Singing Style with Multi-Singer Singing Synthesis System. In: 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 7224–7228 (2020) https://doi.org/10.1109/ICASSP40776.2020.9054636
  * [10] Luo, Y.-J., Hsu, C.-C., Agres, K., Herremans, D.: (2019). Singing Voice Conversion with Disentangled Representations of Singer and Vocal Technique Using Variational Autoencoders. arXiv preprint arXiv:1912.02613 (2019)
  * [11] Meseguer-Brocal, G., Cohen-Hadria, A., Peeters, G.: Creating DALI,a Large Dataset of Synchronized Audio, Lyrics, and Notes. In: Transactions of the International Society for Music Information Retrieval, vol. 3 pp. 55-67 (2020) https://doi.org/10.5334/tismir.30
  * [12] Nachmani, E., Wolf, L.: Unsupervised Singing Voice Conversion. arXiv preprint arXiv:1904.06590 (2019)
  * [13] Nercessian, S.: Zero-shot Singing Voice Conversion. In: Proceedings of the International Society for Music Information Retrieval Conference (2020)
  * [14] O’Connor, B., Dixon, S., Fazekas, G.: An Exploratory Study on Perceptual Spaces of the Singing Voice. In: The 2020 Joint AI Conference on Music Creativity, vol. 1, Stockholm, Sweden (2020)
  * [15] Jia, Y., Zhang, Y., Weiss, R. J., Wang, Q., Shen, J., Ren, F., Wu, Y.: Transfer Learning from Speaker Verification to Multispeaker Text-To-Speech Synthesis. arXiv preprint arXiv:1806.04558 (2019)
  * [16] Kameoka, H., Kaneko, T., Tanaka, K., Hojo, N.: ACVAE-VC: Non-parallel many-to-many voice conversion with auxiliary classifier variational autoencoder. arXiv preprint arXiv:1808.05092 (2020)
  * [17] Qian, K., Zhang, Y., Chang, S., Yang, X., Hasegawa-Johnson, M.: AutoVC: Zero-Shot Voice Style Transfer with Only Autoencoder Loss. In: 36th International Conference of Machine Learning, vol. 47, pp. 5210-5219 (2019)
  * [18] Qian, K., Jin, Z., Hasegawa-Johnson, M., Mysore, G. J.: F0-Consistent Many-To-Many Non-Parallel Voice Conversion Via Conditional Autoencoder. In: 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 6284–6288 (2020) https://doi.org/10.1109/ICASSP40776.2020.9054734
  * [19] Qian, K., Zhang, Y., Chang, S., Cox, D., Hasegawa-Johnson, M.: Unsupervised speech decomposition via triple information bottleneck. In: 37th International Conference on Machine Learning, vol. 119, pp. 7792–7802 (2020)
  * [20] Raffel, C., Ellis, D. P. W.: Feed-Forward Networks with Attention Can Solve Some Long-Term Memory Problems. arXiv preprint arXiv:1512.08756 (2016)
  * [21] Ronneberger, O., Fischer, P., Brox, T: U-Net: Convolutional Networks for Biomedical Image Segmentation. In: N. Navab, J. Hornegger, W. M. Wells, A. F. Frangi (eds.) Medical Image Computing and Computer-Assisted Intervention, vol. 9351, pp. 234–241. (2015) https://doi.org/10.1007/978-3-319-24574-428
  * [22] Tolstikhin, I., Bousquet, O., Gelly, S., Schoelkopf, B.: Wasserstein Auto-Encoders. arXiv preprint arXiv:1711.01558 (2019)
  * [23] Veaux, C., Yamagishi, J., MacDonald, K.: CSTR VCTK Corpus: English multi-speaker Corpus for CSTR Voice Cloning Toolkit. (2017) https://doi.org/10.7488/ds/2645
  * [24] Wan, L., Wang, Q., Papir, A., Moreno, I. L.: Generalized End-to-End Loss for Speaker Verification. In: 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 4879–4883 (2018) https://doi.org/10.1109/ICASSP.2018.8462665
  * [25] Wilkins, J., Seetharaman, P., Wahl, A., Pardo, B. (2018).: VocalSet: A Singing Voice Dataset. 19th International Society for Music Information Retrieval Conference, Paris ,France pp. 468–474 (2018)
  * [26] Wu, D.Y., Chen, Y.H., Lee, H.Y.: One-Shot Voice Conversion by Vector Quantization. In: 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (2020)