Thibault Jougla

David G Dritschel

On the energetics of a two-layer baroclinic flow

, a vertical shear arising from latitudinal heating variations is imposed on the flow and maintained by thermal damping. Idealised convection between the upper and lower layers is implemented by adding cyclonic/anti-cyclonic pairs, called hetons, to the flow, though the qualitative flow evolution is evidently not sensitive to this or other small-scale stochastic forcing. A very wide range of simulations have been conducted. A characteristic simulation which exhibits alternation between two different phases, quiescent and turbulent, is examined in detail. We study the energy transfers between different components and modes, and find the classical picture of barotropic/baroclinic energy transfers to be too simplistic. We also discuss the dependence on thermal damping and on the imposed vertical shear. Both have a strong influence on the flow evolution. Thermal damping is a major factor affecting the stability of the flow while vertical shear controls the number of jets in the domain, qualitatively through the Rhines scale L Rh = U/β.

Introduction

General features of planetary atmospheres such as alternating bands and currents, exemplified by the jet streams on Earth and more prominently on Jupiter and the giant gas planets, are known to be strongly affected by planetary rotation, specifically the latitudinal variation in planetary potential vorticity (PV). Nowhere else is this better seen than on Jupiter, a rapidly rotating planet with a particularly active atmosphere [START_REF] Rogers | The giant planet Jupiter[END_REF]. Observations of Jupiter in fact go back to ancient times. Babylonian astronomers first recorded the appearance of the planet [START_REF] Sachs | Babylonian observational astronomy[END_REF], and for them it was a bright star in the celestial sky. In 1630 with the improvement of optics, Jupiter's bands were first captured by Niccolo Zucchi. Then, Hooke and Cassini in 1664 and 1665 detected spots, now known as vortices [START_REF] Rogers | The giant planet Jupiter[END_REF]. Hence, the presence of bands and vortices was already known at this time. Much more recently, spacecraft observations have revealed complex turbulent dynamics occurring in Jupiter's atmosphere, composed of quasi-zonal jets and myriads of vortices like the Great Red Spot (see [START_REF] Dowling | Potential vorticity and layer thickness variations in the flow around jupiter's great red spot and white oval bc[END_REF] and [START_REF] Marcus | Jupiter's great red spot and other vortices[END_REF]).

The outer Jovian atmosphere is structured in zonal bands, called belts and zones.

They are separated by jets which are strong currents flowing eastwards around the planet. Contrary to Earth's jets, Jovian jets are remarkably straight, with few large scale meanders. Hubble offers us a visualisation of the Jovian atmosphere's external layer where jets and vortices are clearly present (see [START_REF] Simon | First results from the hubble opal program: Jupiter in 2015[END_REF] † and for previous observations see [START_REF] Limaye | Jupiter: New estimates of the mean zonal flow at the cloud level[END_REF] and [START_REF] Porco | Cassini imaging of jupiter's atmosphere, satellites, and rings[END_REF]). Observations of Jupiter have in turn inspired researchers to develop idealised models of the Jovian atmosphere (see e.g. [START_REF] Dowling | Dynamics of jovian atmospheres[END_REF] and [START_REF] Ingersoll | Dynamics of Jupiter's atmosphere[END_REF]). Such models generally differ as regards the vertical structure of the Jovian atmosphere. Actual data are sparse: we have only a single vertical profile collected by the Galileo probe (see [START_REF] Atkinson | The galileo probe doppler wind experiment: Measurement of the deep zonal winds on jupiter[END_REF]), which cannot be argued to be representative. To gain at least a qualitative understanding of the dynamics of Jupiter and of the gas giants generally, two main modelling approaches have been taken: a shallow model and a deep model. The shallow model (typically employing the rotating shallow water equations in a single layer) suggests that jets emerge essentially through an upscale cascade of geostrophic turbulence [START_REF] Williams | Planetary circulations: 1. Barotropic representation of Jovian and terrestrial turbulence[END_REF], [START_REF] Cho | The emergence of jets and vortices in freely-evolving shallow-water turbulence on a sphere[END_REF], [START_REF] Kaspi | Formation of jets by baroclinic instability on gas planet atmospheres[END_REF], [START_REF] Showman | Numerical simulations of forced shallow-water turbulence: effects of moist convection on the large-scale circulation of Jupiter and Saturn[END_REF], [START_REF] Scott | Non-robustness of the two-dimensional turbulent inverse cascade[END_REF]). By contrast, the deep model suggests that the flow is organised in deep concentric rotating cylinders, with jets their external representation, (see [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF], [START_REF] Busse | A simple model of convection in the jovian atmosphere[END_REF]). [START_REF] Kaspi | The deep wind structure of the giant planets: Results from an anelastic general circulation model[END_REF]) studied a 3D non-hydrostatic general circulation model, with a large density variation under the anelastic approximation. Unlike in models employing the Boussinesq approximation valid for weak density variations, they observe the development of baroclinic shear (a strong increase of horizontal wind speed with height). Likewise, [START_REF] Liu | Mechanisms of jet formation on the giant planets[END_REF] and [START_REF] Spiga | Waves and eddies simulated by highresolution global climate modeling of saturn's troposphere and stratosphere[END_REF] used general circulation models to study the atmospheres of the gas giants. [START_REF] Liu | Mechanisms of jet formation on the giant planets[END_REF] were able to reproduce main features of the large-scale circulation on both Jupiter and Saturn, and found that zonal jets extend down to a depth where drag is assumed to become important. [START_REF] Spiga | Waves and eddies simulated by highresolution global climate modeling of saturn's troposphere and stratosphere[END_REF] focused on Saturn and included specifically the roles of aerosol layers, ring shadowing and internal heat fluxes. They found that zonal jets develop in the troposphere and equatorial oscillations occur in the stratosphere. In summary, it has been argued that many models, having widely different formulations, capture qualitative features of planetary circulations on the gas giants.

On Earth, both the atmosphere and the oceans exhibit similar dynamical structures. The jets occurring in the oceans, however, are typically strongly meandering and unsteady (see [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF], [START_REF] Kamenkovich | Role of eddy forcing in the dynamics of multiple zonal jets in a model of the north atlantic[END_REF] and [START_REF] Berloff | On latency of multiple zonal jets in the oceans[END_REF]). The most well-known oceanic jets are the strong oceanic currents like the Gulf Stream and the Kuroshio. The jets occurring in the atmosphere also exhibit large-scale meanders but are generally much more well-defined than their oceanic counterparts. [START_REF] Arbic | Baroclinically unstable geostrophic turbulence in the limits of strong and weak bottom ekman friction: Application to midocean eddies[END_REF]) and [START_REF] Thompson | Two-layer baroclinic eddy heat fluxes: zonal flows and energy balance[END_REF] studied the impact of bottom friction on oceanic simulations. [START_REF] Arbic | Baroclinically unstable geostrophic turbulence in the limits of strong and weak bottom ekman friction: Application to midocean eddies[END_REF] examined the energetics behind the formation of mid-ocean eddies, while [START_REF] Thompson | Two-layer baroclinic eddy heat fluxes: zonal flows and energy balance[END_REF] focused on the role of the baroclinic (vertically varying) mode on the eddy heat flux. [START_REF] Venaille | Ribbon turbulence[END_REF], using a two-layer quasi-geostrophic model on the f -plane, found that for a small bottom friction the flow forms approximately barotropic (height independent) large-scale structures, which persist over a time-scale inversely proportional to the bottom friction damping rate. However, for high bottom friction, the upper layer flow dominates at leading order in the inverse damping rate. In this case, the upper flow behaves like a single-layer shallow-water flow (a 1 1 2 -layer model) and exhibits coherent jets separating regions of nearly homogeneous PV.

To develop a conceptual understanding of Jovian jets and vortices, tyre (2015) recently proposed a 1 1 2 -layer quasi-geostrophic model based on the Dowling and Ingersoll model (see [START_REF] Dowling | Jupiter's great red spot as a shallow water system[END_REF]), with pre-defined jets (created by an imposed lower boundary shape interface) and forced by the injection of cyclone/anticyclone pairs. With this idealised Jovian weather-like model, they managed to reproduce observed features such as nearly steady and straight (zonal) jets together with long-lived vortices.

The goal of the present study is to better understand the vertical structure of jets and vortices in a turbulent planetary atmosphere. A two-layer quasi-geostrophic shallowwater β-channel model is used to examine the formation and evolution of jets and vortices, and in particular their sensitivity to vertical shear, forcing and damping. Following [START_REF] Panetta | Baroclinic eddy fluxes in a one-dimensional model of quasi-geostrophic turbulence[END_REF], a vertical shear is imposed by relaxing to a thermal equilibrium temperature gradient -however we do not include Ekman damping or any bottom friction (see [START_REF] Venaille | Ribbon turbulence[END_REF] for the impact of bottom friction on a two-layer quasigeostrophic model). Bottom friction is not relevant in shallow models of the gas giants [START_REF] Thomson | Jupiter's unearthly jets: a new turbulent model exhibiting statistical steadiness without large-scale dissipation[END_REF]. The vertical shear induces baroclinic instability, whose nonlinear equilibration attempts to reduce the shear. This competition between thermal forcing and baroclinic instability gives rise to the formation of baroclinic jets and, in extreme cases, to stepped PV profiles or 'staircases'. Such staircase formation has already been demonstrated in the single-layer barotropic context (see [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF] and [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]). Additionally, in the present model small-scale stochastic forcing is imposed, crudely mimicking convective processes. This is done by adding hetons or small cyclonic/anti-cyclonic pairs in the two layers. Unlike in the model of [START_REF] Thomson | Jupiter's unearthly jets: a new turbulent model exhibiting statistical steadiness without large-scale dissipation[END_REF], jets do not emerge from an imposed lower boundary shape but emerge naturally from the competition between baroclinic instability and thermal damping.

Section §2 presents the model, the equations, and the boundary conditions. Section §3 presents the results, starting with an analysis of a characteristic simulation, followed by a detailed examination of energy transfers in section §4. We then focus on a turbulent event in §5 and discuss the wider parameter dependence in §6. Finally, section §7 presents our conclusions and a few ideas for future research.

Model formulation

To study the formation and evolution of jets and vortices in turbulent planetary atmospheres, we make use of a two-layer quasi-geostrophic (QG) β-channel model (see [START_REF] Phillips | A simple three-dimensional model for the study of large-scale extra-tropical flow patterns[END_REF]), with a rigid bottom and a free upper surface (see figure 1). It is governed by the set of equations:

Dq 1 Dt = F 1 ; Dq 2 Dt = F 2 (2.1a
)

q 1 = βy + ∇ 2 ψ 1 - f 0 δ 1 H 1 ; q 2 = βy + ∇ 2 ψ 2 - f 0 (δ 2 -δ 1 ) H 2 (2.1b) δ 1 = f 0 (ψ 1 -αψ 2 ) g(1 -α) ; δ 2 = f 0 ψ 2 g (2.1c)
where D/Dt = ∂/∂t+u•∇ is the material derivative, q i is the quasi-geostrophic potential vorticity (QGPV) in layer i, F i includes all forcing and damping on the QGPV in layer i, β is the linear gradient of the Coriolis frequency f = f 0 + βy, ψ i is the streamfunction in layer i, δ i is the displacement of the upper interface of layer i, H i is the mean depth of layer i, ρ i the uniform density in layer i, α = ρ 2 /ρ 1 is the density ratio, and g is the acceleration due to gravity. We eliminate all reference to f 0 , g and the total mean depth H = H 1 + H 2 by specifying the mean baroclinic Rossby deformation wavenumber, kd , defined through

k2 d = 1 L 2 d = f 2 0 H g(1 -α)H 1 H 2 (2.2)
and making use of the fractional mean layer depths h i = H i /H whose sum is unity (h 1 + h 2 = 1). Then the only explicit parameters are h 1 , kd , α and β.

As explained e.g. in [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF], the QGPV expressions in equations 2.1 derive from a series expansion of the shallow-water Rossby-Ertel PV, q i = (ζ i + f )/(h i (1 + hi )), where ζ i is the relative vorticity in layer i, and hi is the fractional layer depth perturbation defined by hi = (δ i -δ i-1 )/H i , with δ 0 = 0. The leading-order non-constant terms, for ζ i /f 1 and hi 1, define the QGPV:

q i = βy + ζ i -f 0 hi . (2.3)
The interface displacements δ i are related to the layer pressure perturbations p i through hydrostatic balance (cf. [START_REF] Mohebalhojeh | Contour-advective semi-lagrangian algorithms for many-layer primitive-equation models[END_REF]):

p i = g n j=i (ρ j -ρ j+1 )δ j (2.4)
where here n = 2 and ρ n+1 = 0. From geostrophic balance, we have

p i = ψ i f 0 ρ i (2.5)
together with u i = -∂ψ i /∂y and v i = ∂ψ i /∂x. Using 2.4 and 2.5, we can re-express 2.3 for each layer as

q 1 = βy + ∇ 2 ψ 1 + h 2 k2 d (αψ 2 -ψ 1 ) ; q 2 = βy + ∇ 2 ψ 2 + h 1 k2 d (ψ 1 -ψ 2 ).
(2.6)

These may be regarded as 'inversion relations' providing the flow field u i in each layer (via ψ i ) from the instantaneous distribution of QGPV (q 1 and q 2 ).

Stratification

Atmospheric stratification is modelled by an exponential decrease of the density profile ρ = ρ 0 e -z/Hρ , where ρ 0 is the density at the bottom of the domain and H ρ is the scale height (such a profile corresponds to an isothermal basic state), see [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF], p.40. The impact of strong variations in density stratification has previously been studied by [START_REF] Fu | Nonlinear energy and enstrophy transfers in a realistically stratified ocean[END_REF] and [START_REF] Smith | The scales and equilibration of midocean eddies[END_REF] in a two-layer quasi-geostrophic ocean model. Such variations inhibit energy transfers to the barotropic mode, thereby favouring a relatively strong upper layer flow. The same enhanced baroclinic response is found in the present results below.

The scale height is defined through c ρ H ρ = H 1 + H 2 = H, with c ρ being the depth of the model in scale heights. We take equal layer depths H 1 = H 2 = c ρ H ρ /2, and define ρ i to be the mean density in each layer:

ρ 1 = H -1 1 H1 0 ρ dz = ρ 0 H ρ H 1 (1 -e - H 1
Hρ ) (2.7a)

ρ 2 = H -1 2 H H1 ρ dz = ρ 0 H ρ H 2 (e - H 1
Hρ -e

-H Hρ ) (2.7b)
Then the density ratio α = ρ 1 /ρ 2 is prescribed through α = e -cρ/2 . The densities ρ 1 and ρ 2 are not explicitly needed in the model, only their ratio α is.

Stochastic Forcing

Material changes of the QGPV in each layer occur through stochastic forcing and thermal damping, represented by the terms F 1 and F 2 in (2.1). The stochastic forcing models in a very simple way unresolved convective motions through the spatially-random addition of 'hetons', which helps to destabilise the flow. A heton is a pair of oppositesigned PV anomalies (cyclonic in the lower layer and anti-cyclonic in the upper layer, resulting from convergence and divergence respectively) which carry heat (see [START_REF] Carton | Hydrodynamical modeling of oceanic vortices[END_REF]). In order to avoid any net vorticity input in either layer, a compensating uniform vorticity is added (mimicking subsidence). Each heton has a fixed radius R = 0.05 and PV amplitude q heton on the highest baroclinic mode (there is no projection on the lowest mode). The frequency at which hetons are added is controlled by a prescribed (potential) enstrophy input rate η. The dependence of the flow evolution on η has been studied by using values ranging from 0.01 to 100. The dependence on the heton's radius R has also been examined but is found to have a minor effect. Notably, simulations initialised with small non-zonal perturbations and no heton forcing produce qualitatively similar results.

Thermal Damping and Vertical Shear

The thermal forcing of a planetary atmosphere generally varies from the equator to the poles. Equatorial regions typically receive more solar radiation than the poles (at lower levels). The resulting latitudinal temperature gradient implies a vertical gradient of the horizontal velocity, by thermal wind balance under the QG approximation. As in [START_REF] Panetta | Baroclinic eddy fluxes in a one-dimensional model of quasi-geostrophic turbulence[END_REF], this vertical shear is represented in a 2-layer QG model as a uniform westward flow in the lower layer and a uniform eastward flow in the upper layer, as illustrated in figure 2. Indeed, a negative northward temperature gradient at lower levels (in the northern hemisphere) implies a positive vertical gradient of velocity, here westward in the lower layer and eastward in the upper layer. The mean zonal velocity is kept fixed at its thermal equilibrium value at the y-boundaries. In the interior, the mean zonal velocity may vary as a result of baroclinic and barotropic instability. From equation 2.6, the QGPV is linearly related to the streamfunction, which in turn is linearly related to the velocity. Adding a uniform velocity to the flow corresponds to adding a linear streamfunction profile ψ ∝ y and hence a linear PV profile in 2.6. Thus, vertical shear corresponds to adding a linear profile to the pre-existing background PV, βy, so that the total PV is, say, i βy in each layer (i = 1, 2). Below, we take 1 to be a control parameter and deduce 2 from the requirements that (i) the vertical shear projects only on the vertical mode with the highest deformation wavenumber (hereafter, the 'baroclinic' mode), and that (ii) it has zero mass average (ρ 1 h 1 ū1 + ρ 2 h 2 ū2 = 0). Full details may be found in appendix §A.

To allow baroclinic instability, the potential vorticity gradients in the two layers must have opposite signs, implying 1 2 < 0. Baroclinic instability however acts to reduce the vertical shear, leading to a quasi-stable flow. To prevent this and to allow cycles alternating between turbulent and quiescent states, thermal damping is imposed to help maintain the large-scale vertical shear. The thermal damping here crudely corresponds to radiative forcing by the sun or internally. It acts to restore the layer interfaces δ i to their initial thermal equilibrium profiles δ eq,i at the rate r. This leads to nonconservative changes in the QGPV through the relation between the QGPV and interface displacement, namely (2.3) with hi = (δ i -δ i-1 )/H i and δ 0 = 0. To eliminate reference to f 0 , g etc., it is convenient to scale the interface displacements so they have units of vorticity:

δi ≡ f 0 δ i H .
(2.8)

Owing to the relations (2.1), together with the definition of k2 d in (2.2), we can express the scaled interface displacements as

δ1 = h 1 h 2 k2 d (ψ 1 -αψ 2 ) ; δ2 = h 1 h 2 k2 d (1 -α)ψ 2 .
(2.9)

Then, thermal damping contributes the following terms to F 1 and F 2 in (2.1):

F 1 ← r δ1 -δeq1 /h 1 ; F 2 ← r ( δ2 -δ1 ) -( δeq2 -δeq1 /h 2 .
(2.10)

Boundary Conditions

We employ a channel model, with periodic conditions in x (longitude) and free-slip rigid boundaries in y (latitude). On the y-boundaries, the zonally-averaged zonal velocity is held fixed to the value imposed by the vertical shear, as in [START_REF] Panetta | Baroclinic eddy fluxes in a one-dimensional model of quasi-geostrophic turbulence[END_REF]. The model has a rigid bottom surface and a free top surface, as illustrated schematically in figure 1.

Non-dimensionalisation

We choose to employ dimensionless quantities only. To this end, we take the width of the domain to be L x = 2π and its breadth to be L y = π, and, most importantly, ensure that these dimensions are much larger than the characteristic Rossby deformation length L d = 1/ kd = 0.05. † The total height of the model H is taken to be c ρ = 2 density scale heights H ρ , corresponding to a density ratio of ρ 2 /ρ 1 = α = e -cρ/2 0.37, though c ρ = 0.2 and 0 have also been studied (see below). For c ρ = 2 and kd = 20, and for equal layer depths h 1 = h 2 = 0.5, the deformation wavenumbers of each vertical mode are k d1 8.87 and k d2 17.93 (see appendix §A for details).

Parameter choices

We set the planetary vorticity gradient β to be 8π, corresponding to a Rossby wave period of unity for a barotropic disturbance with a wavenumber of 4. Recall that the initial PV gradients in each layer are i β, and the choice of 1 controls the strength of the vertical shear. Here, we consider two main values of 1 , namely -0.1 and -0.5, corresponding to 2 2.81 and 3.79 respectively. For 1 = -0.5, the initial mean velocities in the two layers are ū1 -0.094 and ū2 0.256 (see appendix §A for details). The enstrophy input rate η controlling the injection of hetons ranges between 0.01 and 100. We consider values of the heton PV q heton ranging from 1 × β/ kd to 1200 × β/ kd , whereas we fix the heton radius R at 0.05 unless otherwise stated. We vary the thermal damping rate r between 0 and 0.5. Section §6 examines the effect of variations in the parameters. A full list of the dimensionless physical parameters and the values used are provided in Table 1.

Over 100 simulations have been run using a variety of grid dimensions N x × N y , for N x = {128, 256, 512} and N y = N x /2 + 1, over 10 4 units of time (and up to 10 5 units for very small r). Note: an extra grid point is needed in the y direction to include the boundaries and ensure each grid box is square.

Numerical Model

All simulations have been carried out using the 'Combined Lagrangian Advection Method' (CLAM) developed by [START_REF] Dritschel | The combined Lagrangian advection method[END_REF]. This numerical method, using a pseudo-spectral method at large scales, and Lagrangian contour advection at small scales, is both highly accurate and efficient (see [START_REF] Dritschel | Two-dimensional magnetohydrodynamic turbulence in the small magnetic prandtl number limit[END_REF] for a recent demanding comparison). All numerical parameters follow the recommended choices outlined in [START_REF] Fontane | The hypercasl algorithm: a new approach to the numerical simulation of geophysical flows[END_REF] except that twice as many PV contour levels (here 80) are used to represent the initial PV variation in each layer. The PV contour interval remains fixed in time.

Results

We begin by presenting in detail a characteristic flow simulation. The flow illustrated alternates, irregularly, between a more disturbed 'turbulent' phase and a more zonal 'quiescent' phase. The initial baroclinically-unstable zonal state rapidly breaks down and never fully recovers. The instability equilibrates into a structured flow, a quiescent phase, where regions of nearly homogeneous PV are separated by eastward jets in the upper layer and westward jets in the lower layer. Next, and as a result of the accumulative effect of thermal damping, this structured state destabilises and breaks down, again through baroclinic instability but on a much different basic state than was present initially. The flow becomes more turbulent and much less organised. It then recovers to a structured quiescent phase only to break down again into a turbulent phase, and so on. We pay particular attention to the flow behaviour around these turbulent phases, and we study the changes in energy (kinetic/potential, zonal/eddy, barotropic/baroclinic) and zonalflow stability which occur.

A characteristic simulation

The chosen simulation has been run on a 512 × 256 'inversion' grid (the effective grid resolution is 16 times finer in each direction), with a thermal damping rate r = 0.01, a stratification parameter c ρ = 2, a lower-layer fractional PV gradient 1 = -0.5, a heton PV q heton = 1 × β/k d , and an enstrophy input rate η = 0.1.

At the initial time, the flow exhibits a linear slope in PV, opposite in each layer. This flow is baroclinically unstable (as explicitly shown in appendix §A), and within a few units of time leads to a strongly non-zonal, eddying, turbulent flow. This subsequently collapses and organises into bands of nearly uniform PV separated by jets. This quiescent phase, illustrated in figure 3 at t = 3000 typically lasts hundreds of time units, or several damping time scales r -1 .

The central panels in figures 3a and 3b show the PV field at this time in the upper and the lower layers, respectively. In the upper layer, four distinct bands of nearly homogeneous PV are present, separated by three jets flowing eastward. In the lower layer, by contrast, there are three bands of nearly homogeneous PV separated by two westward jets. As discussed in [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF] and [START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF] in the context of a single-layer model, the 'equivalent latitude' PV profile y e (q) is particularly useful for identifying jets, here seen by the nearly flat portions of the curves shown in the right panels of these figures. The formation of a staircase profile in y e (q) indicates the presence of homogeneous regions (dy e /dq → ∞) and jets (dy e /dq → 0). The function y e (q) is obtained by re-arranging the PV monotonically in each layer. In the upper layer, y e (q) = y min + A 2 (q)/L x where A 2 (q) is the area occupied by PV values having Left: normalised latitude 2y/π vs zonally-averaged zonal velocity ūi (y, t). Centre: PV field q i (x, y, t) over the entire domain. Right: equivalent latitude y e (q, t) vs normalised PV q = (q -q min )/(q max -q min ).

q > q min . In the lower layer, where the mean PV gradient is generally reversed, instead y e (q) = y min + A 1 (q)/L x where A 1 (q) is the area occupied by PV values having q < q max . The right panels in figures 3a and 3b show y e (q) versus normalised PV, defined by q ≡ (q -q min )/(q max -q min ), for each layer. Both layers exhibit a near staircase profile, but it is most distinct in the upper layer. The three main small gradient portions of the curve seen here correspond to the strong eastward jets seen in the zonally-averaged zonal velocity ū2 (upper left figure), while the four high gradient portions correspond to the nearly homogeneous PV regions. The zonal flow in these regions is weakly westwards. In the lower layer, the two (central) small gradient regions in y e (q) correspond to westward jets (see lower left panel). Between the westward jets in this layer one also sees eastward jets of comparable magnitude. These are induced by the strong PV gradients in the upper layer. PV inversion couples the layers together so that the flow in either layer depends on the PV in both layers. Here however the lower layer PV is relatively weak, so it has much less impact on the upper layer flow than has the upper layer PV on the lower layer flow.

Overall, the observed jets are not vertically coherent, implying that the baroclinicity of the flow is important. The two layers have different numbers of jets. Jets flow eastward in the upper layer and predominantly westward in the lower layer. Moreover, the upper layer jets are much faster than the lower layer jets. Additionally, signatures of the jets in one layer are visible in the PV field of the other layer. In the upper layer (central panel), the transitions between the two sky/pale blue bands, and between the pink-reddish bands (the two middle bands) coincide with the locations of the lower layer westward jets. These signatures appear in the equivalent latitude y e (q) (upper right panel) as small swerves in the middle of the two high gradient portions. In the lower layer, the signatures of the upper layer jets are less obvious, but are visible in the zonally averaged flow, an effect of PV inversion mentioned above. Jet signatures in PV are less evident in the lower layer in part because the flow is more disturbed at small scales there.

Figures 4a and4b illustrate the same flow at a later time t = 4000, during a turbulent phase. The flow is substantially less zonal, with large-amplitude excursions exhibited in layer flow, all at t = 9400. Left: normalised latitude vs zonally-averaged zonal velocity (Note that plot scale for u are the same as in figure 3). Centre: PV field over the entire domain. Right: equivalent latitude ye(q, t) vs normalised PV.

the PV field in both layers. Remarkably, nearly homogeneous regions persist in the upper layer though they are much more disturbed than in the quiescent phase. Moreover, the jets are considerably weaker and virtually absent from the lower layer flow. The 'jets' seen in the upper layer are not steady and show no obvious relation to the equivalentlatitude PV profile y e (q, t). A series of short-lived vortices can be seen, especially in the central regions of the flow. These formed from the buckling and breaking of the jets at earlier times and re-merge with the jets subsequently. Two movies of the PV evolution are available in the online supplementary material.

Various statistics summarising the time evolution of the flow are provided in figure 7. From top to bottom we show the equivalent PV q e (y) (the inverse of the function y e (q)) in the upper and lower layers as Hövmoller diagrams, the zonal and eddy kinetic energy in the first 'barotropic' and second 'baroclinic' vertical modes, the zonal and eddy potential energies, the maximum growth rate of the instantaneously zonally-averaged flow versus zonal wavenumber, and the meridionally-integrated spectrum of the available potential energy (APE) versus zonal wavenumber. In each plot, time is in the abscissa.

The equivalent PV Hövmoller diagrams allow one to identify major changes in the flow structure, such as a latitudinal shift of the jets and homogeneous regions. We can see here the entire evolution of the flow, starting with the first rapid re-organisation of the flow from its unstable initial state, to the alternation between quiescent phases and (relatively short) turbulent phases.

The initial linear PV variations i βy are short lived after the flow destabilises and evolve into nearly homogeneous regions separated by jets, or staircase profiles, by t = 500. This is the beginning of the first quiescent phase. Four nearly homogeneous regions are clearly identifiable in the upper layer and between them three jets are present. The lower layer at this stage has three nearly homogeneous regions with two westward jets separating them. Over the next 1600 units of time, the innermost homogeneous regions narrow slowly but progressively, due to the accumulative effects of thermal damping (this does not occur in the absence of damping). Near t = 2100, the quasi-zonal flow is disrupted and a short turbulent phase occurs, lasting no more than approximately 100 days, which is the thermal damping time-scale, r -1 . The flow then recovers, this time significantly more rapidly than it did following the initial instability, and enters a second quiescent phase, illustrated in figure 3 at t = 3000.

Over the next 1900 units of time, the same progressive thinning of the innermost homogeneous regions occurs, disrupted by a turbulent phase around t = 4000, illustrated in figure 4. This sequence of long quiescent and short turbulent phases continues with irregular frequency almost to the end of the simulation. At around t = 9000, the flow passes through a turbulent event and organises into an anomalous quiescent phase, similar to that seen after the initial instability near t = 0, see figure 5. Instead of having a symmetric flow with a centered jet, the flow exhibits a shift of the jets and homogeneous regions as shown in figure 6.

This anomalous quiescent phase is significantly more active and non-zonal than previous such phases. It is also found in other simulations and is discussed further in §B below.

The energy diagnostics closely correlate with the variations seen in the equivalent PV Hövmoller diagrams. Each turbulent phase is characterised by a strong peak in the nonzonal, eddy component of the energy and a dip in the zonal component. For the barotropic (mode 1) kinetic energy, the eddy energy (EKE1) is always comparable or greater than the zonal part (ZKE1). The baroclinic (mode 2) kinetic energy is smaller, with the eddy part (EKE2) surpassing the zonal part (ZKE2) only during turbulent phases. On the other hand, the potential energy is much larger: even the eddy potential energy (EPE) is greater than any component of the kinetic energy. The zonal potential energy (ZPE) is more than 10 times any component of kinetic energy. This dominance of ZPE is seen in all simulations conducted. Note: here we mean the 'available' potential energy (APE) relative to a resting basic state (APE measures the mean-square displacement of the layer interfaces).

The maximum growth rate of the instantaneous zonally-averaged flow (second panel from the bottom in figure 7) shows that the turbulent phases coincide with periods of relatively strong baroclinic instability, involving zonal wavenumbers k x mainly between 5 and 20. Note that the instability is not in fact purely baroclinic, except in the initial stages when the mean flow is a simple baroclinic shear. In general, there is a projection on the barotropic mode as well when the zonal mean flow is more complex (see e.g. [START_REF] James | Suppression of baroclinic instability in horizontally sheared flows[END_REF], who examined the influence of barotropic shear on baroclinic instability). It appears that the instability precedes the growth in the eddy energy components, but in fact it occurs after these components have already grown significantly (this is shown below in detail). Hence, it is too simplistic to think that the thermal damping has led to an unstable zonal flow state, which then becomes turbulent. The non-zonal flow must also play a role, but this cannot be assessed simply through a linear stability analysis.

Finally, in the bottom panel in figure 7 we show the meridionally-integrated APE spectrum as a function of zonal wavenumber k x (between 1 and 25) and time. At the very earliest stages, t 10, baroclinic instability excites wavenumbers centred on k x = 11 (as found in the linear stability analysis). Very soon after, a broad spectrum of disturbances emerges, with dominant power residing at much lower wavenumbers, around k x = 4 and to a lesser extent k x = 6. This is typical of quiescent phases, and indicates a significant and persistent amount of eddy activity even during these phases. The turbulent phases on the other hand are characterised by the excitation of a broad range of wavenumbers, both lower and higher. This spreading is characteristic of turbulent flow evolution. Each turbulent phase appears to first spread APE to higher wavenumbers then build up energy at low wavenumbers, predominantly k x = 2 and 3. This is short-lived however as the flow returns to a quiescent phase.

Energy Transfers

We next look more closely at the energy transfers taking place between different energy components (kinetic/potential, barotropic/baroclinic, zonal/eddy) to gain a better understanding of the flow changes occurring when entering and leaving turbulent phases. To be more quantitative, we define a turbulent phase as a period when the ZPE is smaller than the mean ZPE. Other definitions have been tried, but this simple definition adequately identifies turbulent phases, as seen in figure 7.

To study energy transfers, we examine two different kinds of instantaneous correlations between the six energy components ZKE1, ZKE2, ZPE, EKE1, EKE2 and EPE. Note that it is not possible to express the rate of change of any one component solely in terms of the other components, so one must seek other means to examine the energy transfers taking place. We have chosen correlations, separated into turbulent and quiescent phases, together with directions of energy change, as explained below.

The first correlation is the standard one, defined for two time series {x 1 , x 2 , ..., x n } and {y 1 , y 2 , ..., y n } by where x = n -1 n i=1 x i and ȳ = n -1 n i=1 y i . The second one, which we call the 'global correlation', is defined by the numerator of (4.1) only. This gives a better measure of the energy magnitudes involved and, in particular, helps to identify the key components involved in the energy transfers taking place. In the results presented, data have been sampled every 5 units of time, providing energy time series of n = 2001 values.

r xy = n i=1 (x i -x)(y i -ȳ) n i=1 (x i -x) 2 n i=1 (y i -ȳ) 2 (4.1)
Figure 8 shows three scatterplots of different energy components and the standard correlation matrix. First, the correlation matrix exhibits three distinct areas. The two red areas show that zonal energies are positively correlated with each other, and that eddy energies are also positively correlated with each other. The blue area shows that zonal energies are negatively correlated with eddy energies. This finding agrees with that found in figure 7 for the evolution of the energy components, namely that when there is a decrease of any zonal energy the other zonal energies also decrease and the eddy energies increase, and vice versa.

In the scatterplots, further information can be obtained. First of all, we see the strong, tight correlation between ZKE2 and ZPE (top left figure), consistent with the value of 0.94 seen in the correlation matrix. Moreover, the eddy energies are all tightly and positively correlated (only EKE2 vs EPE is presented here in the bottom left figure). The scatterplots exhibit a racquet shape, very narrow at the base, corresponding to the quiescent phases, and more spread at the top, corresponding to the turbulent phases. Moreover, there is a definite pattern of evolution: the flow leaves a quiescent phase with relatively low EPE, gains EPE during the peak of the turbulent phase, then returns to the quiescent phase with relatively high EPE.

Similar looping patterns are present in the zonal versus eddy and in the ZKE1 versus zonal scatterplots. These patterns are more widely spread than seen between the eddy energies themselves. If we look at ZPE versus EPE (top right figure), we see that ZPE diminishes as EPE grows when leaving a quiescent phase, following nearly a straight line. It then returns along a curved line to recover from a turbulent phase. During this recovery, EPE falls relatively quickly while ZPE does not show much growth until near the end of the turbulent phase. This slow recovery of ZPE is likely due to the weak thermal damping. The nearly straight lines, leaving quiescent phases and entering turbulent phases, are all closely parallel. The more eccentric one, with black arrows, corresponds to the initial stages where baroclinic instability acts on an idealised basic state. All the other straight lines are parallel to this one because the same phenomenon happens: ZPE is converted into EPE due to baroclinic instability. However the zonal energy at the beginning of subsequent turbulent phases is lower than it was in the initial phase, explaining the gap between the different lines.

A difficulty in the interpretation of these results arises because the scatterplots are twodimensional projections of a six-dimensional phase space. The loop patterns seen imply that other components not shown in a given 2D cross section are involved in the energy transfers taking place, and moreover that the transfers are not occurring simultaneously (this is highlighted in §5). An alternative interpretation, discussed next, focuses on the energy components exhibiting the greatest variations.

The standard correlation scatterplots in figure 8 show the correlation between different energy components, but not the energy transfer. To see the energy transfer, it is necessary to use the same plot limits in all scatterplots and to correlate the energy components using only the numerator of (4.1). These global correlations are exhibited in figure 9 in the same format as figure 8 for comparison. The main energy transfers involve just a few components, ZPE, EPE and to a weaker extent EKE1 and ZKE1. The other components have essentially no impact on the global energy transfer. In particular, some components which show a strong standard correlation instead show a very weak global correlation, such as ZKE2 versus ZPE. This means that ZKE2 and ZPE evolve similarly but there is scarcely any energy transfer between them.

The dominant energy transfer occurs between ZPE and EPE and is strongly anticorrelated, as expected. ZPE is also anti-correlated with EKE1, implying that both EPE and EKE1 are the dominant components uptaking any change in ZPE. This is consistent also with the positive correlation between EPE and EKE1. Finally, and perhaps surprisingly, only one correlation between zonal energies, ZPE versus ZKE1, is notable. For a relatively large gain in ZPE, there is a modest increase in ZKE1. Evidently, the thermal restoration of the sloping layer interfaces favours an increase in barotropic energy. Restoration of a linear slope would not have this effect, as a linear slope only generates a baroclinic flow, by construction. Instead, it must be that the eddying motions present even during quiescent phases continue to provide inhomogeneous PV mixing, maintaining sharp jets [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF][START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]). These jets project on both barotropic and baroclinic modes, unlike the initial vertical shear flow.

Focus on a turbulent phase

The detailed behaviour of the flow around a turbulent phase is examined in this section. The panels in figure 10 show the time evolution of various diagnostics over 1800 t 2800. The first three panels show energy components linearly rescaled to lie between 0 and 1 over the entire simulation duration (0 t 10000). In the third panel, the complement of the rescaled zonal energies is shown so that they can be better compared with the rescaled eddy energies (this panel combines the first two). The fourth panel shows the maximum growth rate for each zonal wavenumber k x (in the abscissa), and the bottom panel shows the meridionally-integrated APE spectrum, also versus k x .

In the top panel, the rescaled zonal energies ZKE2 and ZPE follow the same trend, especially on the approach to the turbulent phase. This is why the standard correlation between ZKE2 and ZPE is very high. However, ZKE1 decreases before the two other energies and recovers sooner. Similarly, in the second panel for the rescaled eddy energies, there is a small offset between the different components. EKE1 and EKE2 grow first (with EKE2 slightly ahead) and next EPE follows, with almost the same trend. This offset between the different energy components explains the looping or racquet patterns observed in the standard correlation scatterplots in Table 1.

The third panel compares the rescaled zonal and eddy energies. This highlights the differences in evolution of the various components. The eddy energies are seen to rise and fall rapidly, over a timescale comparable but shorter than the damping period r -1 = 100. The mode 1 zonal kinetic energy ZKE1 appears to lead all components into the turbulent phase, but due to the re-scaling, the energy changes involved are actually smaller than those taking place in ZPE and EPE (cf. figure 7). The onset of the turbulent phase is best associated with the steep rise in all the curves just before t = 2100, with the dominant energy transfer between ZPE and EPE. Other rising events happen later, e.g. between t = 2300 and 2400, without leading to any major disruption of the flow.

It therefore appears that the disturbance must reach a threshold amplitude before triggering a turbulent phase (this has been confirmed by examining many such events). In the PV field (not shown, but see the movie available to the online supplementary material), wave amplitudes along the jets are seen to slowly grow, eventually (sometimes after several false starts) leading to wave breaking and vortex formation. This strongly disturbed state lasts only a relatively short time before the vortices are either torn apart by the jet shear or recaptured by the jets, and the wave amplitudes gradually subside. The slow return to a quiescent state is most apparent in the behaviour of ZPE and ZKE2: the thermal damping slowly acts to rebuild the sloping layer interfaces, at the same time enhancing the baroclinic shear (in ZKE2). In this return to a quiescent phase, the principal energy changes (i.e. in ZPE) are mainly brought about by thermal damping. However, the rapid decay in all eddy energies, occurring over a period of just 20 to 30 units of time, cannot be explained simply by thermal damping. Instead, part of this decay must be due to inhomogeneous turbulent mixing, constrained by the mean PV gradient in each layer. The mixing acts to homogenise PV and sharpen jets, processes which are only enhanced by non-zonal variations [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF][START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF]. The mixing continues until most of the eddy energy (principally at higher wavenumbers where thermal damping is weak) is exhausted and converted into jets (qualitatively this is the argument behind the Rhines scale L Rh = U/β, where U characterises the eddy velocities and β is the background PV gradient, see [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF]). The low wavenumber part of the eddy energy gives rise to meanders, which unlike the high wavenumber part are affected by the thermal damping. Subsequently, thermal damping weakens the meanders and establishes quasi-zonal jets around t = 2200. Thereafter, thermal damping gradually modifies them, slowly bringing them closer together at the centre of the channel, as seen in the top two panels of figure 7. Eventually, this sets up the conditions for a further turbulent phase.

The maximum growth rate of the flow (shown in the fourth panel of 10) increases in the expected range of baroclinic instability wavenumbers, i.e. wavenumbers around k x = 11, at the onset of the turbulent phase. Thereafter, at the peak of the turbulent phase around t = 2150, a wide range of wavenumbers are excited. Notably, just before this peak there is a decrease in the growth rate for high wavenumbers (18 to 20), that will only recover during the next quiescent phase. While it is not very clear in this figure, this feature is common to all the turbulent phases.

The APE spectrum (bottom panel) before the turbulent phase shows significant and persistent excitation of wavenumbers 4 to 6. During the first part of the turbulent phase, there is an increase in power for small wavenumbers, from 2 to 4, and more generally for all wavenumbers. This is the effect of wave breaking, vortex generation and subsequent turbulent mixing. The increased excitation of all wavenumbers is matched by the growth of EPE shown in the second panel. Likewise, as the excitation subsides, so does the EPE, reaching a minimum around the end of the turbulent phase at t = 2350.

From this analysis, we can better understand the energy variations occurring during the onset and decay of a turbulent phase. The loop patterns exhibited in the correlation scatterplots (cf. Table 1) are associated with the asynchronous growth and decay of the various energy components. What we find generically, across many events and simulations examined, is that at the onset of a turbulent phase ZKE1 decreases marginally sooner than both ZKE2 and ZPE, which remain highly correlated until near the end of the turbulent phase. ZKE1 recovers much sooner than either ZKE2 or ZPE, and this appears to be an eddy effect. Looking back at figure 7, we can see that the dips in ZKE1 occurring during a turbulent phase only last while there is significant eddy energy. By contrast, ZKE2 and ZPE take longer to recover, well after the eddy energies have all but dissipated. Another generic feature worth noting is that, at the onset of a turbulent phase, EKE2 appears to increase slightly before EKE1, which is then followed by EPE. The recovery is also in the same order, with EPE recovering last. The recovery is also rapid, mainly a nonlinear effect of inhomogeneous mixing. Thermal damping acts most strongly at largest scales, and gradually weakens the low-wavenumber meanders on the jets. This helps to re-establish strong quasi-zonal jets, which then slowly shift due to thermal damping until conditions are established for the next turbulent phase.

Dependence on parameters

Parameter Sweep

To study the dependence of the flow evolution on the different parameters involved, a wide range of simulations have been conducted. The main parameters consist of the thermal damping coefficient r, the lower layer fractional PV gradient 1 (directly controlling the vertical shear), the stratification coefficient c ρ , the enstrophy input rate η, and the heton PV q heton . Other parameters have been varied, like the heton radius R and topographic forcing, but these have been found to have only a minor influence. Table 2 lists the parameters used for the main simulations conducted in this study. In the following subsections, we focus on the parameters having the greatest influence: the thermal damping r and the vertical shear parameter 1 . We then briefly summarise the influence of the remaining parameters.

Thermal damping

An important parameter controlling the flow evolution is thermal damping. A wide range of thermal damping rates r = {0, 0.0001, 0.001, 0.01, 0.1, 0.5} have been examined, mainly for different values of stratification coefficient c ρ and vertical shear parameter 1 . Thermal damping acts to restore the initial vertical shear. Simultaneously, baroclinic instability, induced by excessive vertical shear, leads to turbulent mixing which tends to reduce the vertical shear. There is thus a competition between thermal damping and baroclinic instability which nearly always results in the formation of jets through inhomogeneous PV mixing.

Thermal damping is found to have a major impact on the evolution of the flow. This is exhibited in figure 11 comparing the evolution of zonal potential energy (ZPE) for r = 0, 0.001, 0.01 and 0.1, keeping all other parameters at their default values. The green curve corresponds to the characteristic simulation examined already. At high damping rates, r 0.1, the flow appears to be very stable: the initial ZPE is reduced by only 12% for r = 0.1 and fluctuates close to a constant value thereafter. In fact, small- scale disturbances are ever present due to the perpetual creation of baroclinic instability, but no large-scale disruption like that seen in figure 4 for the characteristic simulation ever occurs. At high damping, the flow exhibits many weak jets whose number remains constant from the earliest times. By contrast, with no damping (r = 0), the imposed vertical shear is rapidly cancelled, suppressing any further baroclinic instability. The flow in this case is mainly dominated by large-scale oscillations with a few weak jets, see figure 12. For small but non-zero damping (r = 0.001 and 0.01), there is a competition between the damping trying to restore the vertical shear and baroclinic instability trying to break it down, leading to strong, widely-spaced jets. The damping has a direct impact on the intensity of the jets (their zonal velocity) and jet spacing, as seen in figure 12. As in the characteristic case presented in §3.1, the flow is more unstable and variable, and exhibits large-scale fluctuations. The damping rate strongly affects the recovery time from each turbulent phase, leading to a wider spacing of turbulent events.

Vertical shear

The initial vertical shear is directly controlled by the lower layer relative PV gradient, 1 (see appendix §A). When 1 < 0, PV gradients are reversed in the two layers, and the vertical shear is sufficient to induce baroclinic instability. The more negative 1 is, the greater is the vertical shear and thus the potential for instability.

Vertical shear has a direct influence on the number of homogeneous regions and jets that develop, see figure 13. Increased vertical shear implies stronger baroclinic instability and thus higher eddy velocities, u rms . This on its own would imply an increased jet spacing, qualitatively, through the Rhines scale L Rh = u rms / |q y |, except that the mean PV gradient |q y | also increases. Nevertheless, we find that higher vertical shear leads to both a wider jet spacing and more intense jets, see figure 14. Regarding the oscillations between quiescent and turbulent phases, higher vertical shear decreases their frequency. When the vertical shear is too weak to permit baroclinic instability ( 1 0), the flow evolution is radically different than what we have seen so far. Jets, if at all present, are very weak and erratic.

Other parameters

Many other simulations have been conducted, looking at the effect of stratification, the PV of the injected hetons, the enstrophy input rate, etc. Regarding stratification, the simulations presented in this paper are relevant to an atmospheric-like situation.

We have also run simulations for oceanic-like cases (c ρ → 0). The main difference is a greater propensity for meandering at small scales, such as seen in oceanic currents like the Kuroshio or the Gulf Stream. Large scale structures are still present, though they are weaker and significantly disturbed by small scale perturbations. Jets are also less prominent. As in the oceans, the flow is less structured into zonal bands but contains many more coherent vortices. These findings are broadly consistent with those found by [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF]; [START_REF] Kamenkovich | Role of eddy forcing in the dynamics of multiple zonal jets in a model of the north atlantic[END_REF]; [START_REF] Berloff | On latency of multiple zonal jets in the oceans[END_REF] who examined a forced idealised two-layer ocean with Ekman damping.

Regarding the heton PV and the enstrophy input rate, we have examined the effect of adding a few high intensity hetons versus adding numerous low intensity hetons. We thought that adding high intensity hetons would favour the formation of long-lived vortices, but this is not true. Inevitably any intense vortices which do form drift north or south and get torn apart by or incorporated into the intense jets. Overall, these parameters have only a small impact on the flow evolution. Simulations with different parameters exhibit the same number of homogeneous regions and jets, with comparable jet velocities and spacing.

Conclusions

We have investigated the emergence, nonlinear dynamics and energetics of jets and vortices in a two-layer β-plane quasi-geostrophic channel model. An extensive range of numerical simulations have been carried out using the Combined Lagrangian Advection Model [START_REF] Dritschel | The combined Lagrangian advection method[END_REF], at unprecedented resolution and for very long times (10 4 to 10 5 model days). We have primarily focused on the atmospheric case where there is a strong density difference between the two layers and where thermal damping acts to maintain vertical shear against baroclinic instability. As found in the original study of this type by [START_REF] Panetta | Baroclinic eddy fluxes in a one-dimensional model of quasi-geostrophic turbulence[END_REF], the competition between damping and baroclinic instability results in quasi-zonal jets. In the present study, additionally, we have included the effect of stochastic forcing, specifically adding small-scale baroclinic vortices called 'hetons' randomly at a prescribed enstrophy input rate. These hetons are meant to crudely model the effects of convection thought to occur in the atmospheres of the gas giant planets (see [START_REF] Thomson | Jupiter's unearthly jets: a new turbulent model exhibiting statistical steadiness without large-scale dissipation[END_REF] and references therein). The details of the stochastic forcing however do not seem to matter: qualitatively similar evolution is observed without such forcing when the initial flow is weakly perturbed.

The combination of thermal damping, heton forcing and atmospheric stratification leads to a 'baroclinic life cycle' (see [START_REF] Feldstein | Barotropic decay of baroclinic waves in a two-layer beta-plane model[END_REF][START_REF] Thorncroft | Two paradigms of baroclinicwave life-cycle behaviour[END_REF][START_REF] Esler | The turbulent equilibration of an unstable baroclinic jet[END_REF])), in which jets emerge through baroclinic instability of the vertically-sheared flow and subsequently break down and reform, often in an irregular manner. The jets form through inhomogeneous potential vorticity mixing [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF][START_REF] Scott | The structure of zonal jets in geostrophic turbulence[END_REF], in which potential vorticity is nearly homogenised in bands, but between which the potential vorticity abruptly jumps, forming a 'staircase' profile. Each jump is associated with a jet, eastward in the upper layer and westward in the lower layer, a direct effect of potential vorticity inversion [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF].

Both thermal damping and heton forcing continuously modify the jets that form through baroclinic instability. Thermal damping works to increase the vertical shear (favouring baroclinic instability), while heton forcing creates disturbances which eventually trigger the breakdown of the jets. In most cases, this breakdown is characterised by wave amplification and breaking on the quasi-zonal jets, leading to vortex detachment and turbulence. In this 'turbulent phase', the upper layer jets often persist but become strongly non-zonal and time dependent, while the lower layer jets are nearly obliterated before they reform. The turbulent phase lasts typically one thermal damping period, approximately, but only part of the recovery is due to thermal damping. From a detailed analysis of the energetics, it is found that the eddies (or non-zonal disturbances) created at the onset of a turbulent phase convert their energy into jets, arguably by inhomogeneous potential vorticity mixing (an inviscid mechanism). Thermal damping, by contrast, acts mainly at large scales to dampen long-wave disturbances and straighten the jets. Thereafter, thermal damping slowly shifts the jets and builds the mean vertical shear, setting up the conditions for subsequent instability and breakdown.

Each turbulent phase is followed by a 'quiescent' phase characterised by nearly zonal jets with weak small-scale disturbances. These quiescent phases often exhibit a gradual latitudinal migration of the jets. In one particular example studied in detail, two distinct quiescent phases were observed, one of which was significantly more robust (long lived) than the other (further analysis may be found in appendix B.) The quasi-zonal jets occurring during the quiescent phases are typically highly baroclinic. There are different numbers of jets in each layer located at different mean latitudes. Moreover, the upper layer jets are predominantly eastward, while the lower layer jets are predominantly westward, or a mix between the two due to the strong influence of the upper layer potential vorticity on both layers.

The wider dependence of the flow evolution on thermal damping, vertical shear, stratification, heton intensity and enstrophy input rate has also been studied. Increasing thermal damping weakens the turbulence created during baroclinic instability, leading to less temporal variation especially of the zonal energy components. While baroclinic instability still occurs (for sufficient vertical shear), it only manages to create weak jets for strong thermal damping (i.e. for a 10 day or shorter damping period). On the other hand, no damping results in the destruction of the initial vertical shear and no further chance for baroclinic instability. Weaker damping leads to baroclinic life cycles, here oscillations between turbulent and quiescent phases. Both phases extend in proportion to the damping period, though the recovery of the zonal part is significantly slower than the recovery of the eddy part.

Increasing vertical shear makes the flow more unstable, generally creating fewer, more intense jets. Changing the stratification to a form more appropriate to the oceans, we find that the jets exhibit significant meandering and are less well defined. Moreover, many more vortices are present which last for much longer times than found for atmospheric stratification. These results are consistent with the two-layer ocean study of [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF]; [START_REF] Kamenkovich | Role of eddy forcing in the dynamics of multiple zonal jets in a model of the north atlantic[END_REF]; [START_REF] Berloff | On latency of multiple zonal jets in the oceans[END_REF]. The remaining model parameters (heton intensity, enstrophy input rate) have been found to have only a weak impact on the flow evolution.

A conspicuous feature of the atmospheric-like simulations we have conducted is the absence of long-lived vortices, despite an extensive exploration of parameter space. Vortices do form but they inevitably collide with a jet and are either incorporated or destroyed. A major question then is: is there a parameter regime within this two-layer quasi-geostrophic model favouring long-lived vortices? Or, is it necessary to go beyond the quasi-geostrophic model and consider ageostrophic effects (e.g. as in a shallow-water context)? Or again, is spherical geometry fundamental? Does the Great Red Spot, for instance, depend on its proximity to the equatorial region, whose dynamics can be starkly different from quasi-geostrophic? Are two-layers sufficient? We are presently working to answer some of these questions by building and applying a highly-accurate multi-layer shallow water numerical model. perturbation amplitudes Q(y) and Ψ (y) together with the phase speed c = σ/k x :

(ū i -c)Q i + dq i dy Ψ i = 0, i = 1, 2 (A 20) d 2 Ψj
dy 2 -K 2 j Ψj = Qj , j = 1, 2 (A 21) (A 22) where i refers to layers, j and the hats refer to vertical modes, and K 2 j ≡ k 2 x + k2 dj . To combine these, we eliminate Q i in the first set of equations using Qj in the second set projected onto layers:

Q i = ci1 Q1 + ci2 Q2 = d 2 Ψ i dy 2 -ci1 K 2 1 Ψ1 -ci2 K 2 2 Ψ2 . (A 23)
Next, we eliminate Ψ1 and Ψ2 by projecting them onto modes:

Ψj = c j1 Ψ 1 + c j2 Ψ 2 . (A 24)
Using this in the expression for Q i above, we obtain

Q i = d 2 Ψ i dy 2 -B i1 Ψ 1 -B i2 Ψ 2 (A 25)
where B ij ≡ ci1 K 2 1 c 1j +c i2 K 2 2 c 2j . Finally, inserting this into (A 20) to eliminate Q i results in two coupled, second-order equations for the layer streamfunction amplitudes: This is an eigenvalue problem for the phase speed c. It is solved numerically by discretising the second derivative terms by centred finite differences, and using the boundary conditions Ψ i = 0 on both channel walls. This results in a block-tridiagonal generalised eigenvalue problem, which we have solved using the LAPACK routine DGGEV, The results of the analysis above have been verified by matching them with the exact results obtained for the initially uniform PV gradients dq i /dy = i β and uniform layer velocities ūi , i = 1, 2. In this special case, we can seek solutions for Ψ i ∝ sin(k y (y -y min ))

where k y is the meridional wavenumber (an integer for L y = π). This reduces (A 26) to a simple algebraic system, from which c is determined from a quadratic equation. It is straightforward to show that (baroclinic) instability always occurs when the PV gradients have opposite signs in the two layers, 1 2 < 0 (details omitted).

Figure 15 shows the growth rates as a function of k x and k y , obtained using the parameters of the characteristic simulation. The most unstable disturbances have small meridional wavenumbers k y , and zonal wavenumbers k x between 11 and 13. For a given k x , the maximum growth rate is positive until k x > 13. Hence, there are always disturbances with low k x which are unstable (note: k x = 0 is excluded). These disturbances simply have higher meridional wavenumbers k y . Baroclinic instability occurs in a band of total wavenumbers k 2

x + k 2 y between approximately 11 and 13 here.

Figure 1 :

 1 Figure 1: Vertical cross section of the model. Here, H 1 and H 2 are the mean layer thicknesses, ρ 1 and ρ 2 the layer densities, and δ 1 and δ 2 the displacements of the middle and upper interfaces.

Figure 2 :

 2 Figure 2: Illustration of the initial undisturbed state of a two-layer quasi-geostrophic verticallysheared flow, with uniform westward flow in the lower layer and uniform eastward flow in the upper layer.

Figure 3 :

 3 Figure3: The top figures show the upper layer flow, and the bottom figures show the lower layer flow, all at t = 3000. Left: normalised latitude 2y/π vs zonally-averaged zonal velocity ūi (y, t). Centre: PV field q i (x, y, t) over the entire domain. Right: equivalent latitude y e (q, t) vs normalised PV q = (q -q min )/(q max -q min ).

Figure 4 :

 4 Figure 4: The top figures show the upper layer flow, and the bottom figures show the lowerlayer flow, all at t = 4000. Left: normalised latitude vs zonally-averaged zonal velocity (Note that plot scale for u are the same as in figure3). Centre: PV field over the entire domain. Right: equivalent latitude ye(q, t) vs normalised PV.

Figure 5 :

 5 Figure 5: The top figures show the upper layer flow, and the bottom figures show the lower

Figure 7 :

 7 Figure 7: (a,b) Hövmoller diagrams of the equivalent PV q e (t, y), with latitude in the ordinate. (c,d,e) Energy components (with the eddy part in blue, the zonal part in black and the red line the time averaged APE). (f) log 10 of the maximum growth rate of the zonally-averaged flow versus zonal wavenumber k x . (g) log 10 of the meridionallyintegrated spectrum of available potential energy versus zonal wavenumber k x . The coloured contours correspond to specific values of the spectrum as indicated on the colourbar. In all of these plots time is in the abscissa.

Figure 8 :

 8 Figure 8: (a, b, c) Scatterplots of different energy components in the characteristic simulation, for (a) ZPE vs ZKE2, (b) EPE vs ZPE, and (c) EPE vs EKE2. Note: the mean values are subtracted and the scales are adjusted to the range of energy values observed. The arrows represent a third of the distance between adjacent points. Black arrows: initial phase. Cyan arrows: quiescent phases. Red arrows: turbulent phases. The lines show the best fit of the data (cyan for quiescent phases, red for turbulent phases, and blue overall), after optimal rotation of the data and a least-squares fit. The lengths of the lines are proportional to the spread in the data. (d) The standard correlation matrix is shown, computed using (4.1) for the different energy components.

Figure 9 :

 9 Figure 9: (a, b, c) Scatterplots of different energy components, (a) for ZPE vs ZKE2, (b) forEPE vs ZPE, (c) for EPE vs EKE2, in the characteristic simulation as in Table1but using fixed, identical scales in each plot to emphasise which energy components exhibit the greatest variations. (d) The global correlation matrix is shown, computed using only the numerator of (4.1) for the different energy components.

Figure 10 :

 10 Figure 10: (a,b,c) Rescaled energy components (with mode 1 kinetic energy in black, mode 2 kinetic energy in blue and available potential energy in red) versus time for 1800 t 2800. In (a) only zonal components are shown, in (b) only eddy components are shown, while in (c) all components are shown with the zonal ones inverted for comparison. (d) log 10 of the maximum growth rate of the zonally-averaged flow versus zonal wavenumber k x . (e) log 10 of the meridionally-integrated spectrum of the available potential energy versus k x . The coloured countours correspond to specific values of the spectrum as indicated on the colourbar. In each panel, time is in the abscissa.

Figure 11 :

 11 Figure 11: The evolution of zonal (upper plots) and eddy (lower plots) potential energy (ZPE and EPE) for different thermal damping rates r: no damping, black; 0.001, blue; 0.01 (reference), red; and 0.1, green.

Figure 12 :Figure 13 :

 1213 Figure12: The left figure shows the time-averaged between t = 6500 and t = 7300 (quiescent phase) zonally-averaged zonal velocity ūi(y, t), for different thermal damping rates r: no damping, black; 0.001, blue; 0.01 (reference), red; and 0.1, green.

Figure 14 :

 14 Figure 14: Upper layer (left figure) and lower layer (right figure) zonally-averaged zonal velocity ūi(y, t) averaged between t = 3800 and t = 4500 (quiescent phase), for different vertical shears 1: -0.2 purple; -0.3 red; -0.4 green; and -0.5 orange.

Figure 15 :Figure 16 Figure 17 :

 151617 Figure15: Growth rate as a function of x and y wavenumbers, respectively k x and k y , for the initial flow state used for the characteristic simulation.

Table 1 :

 1 Physical parameters used in the model simulations.

	Parameter	Symbol Value or range considered
	Longitudinal domain length	Lx	2π
	Latitudinal domain length	Ly	π
	Rossby deformation length	L d	0.05
	1st deformation wavenumber	k d1	8.87
	2nd deformation wavenumber	k d2	17.93
	Upper to lower layer density ratio	α	0.37 to 1
	Depth in density scale heights	cρ	0 to 2
	Fractional layer depths	h1, h2	0.5
	Planetary vorticity gradient	β	8π
	Lower layer PV gradient /β	1	-0.5 to -0.1
	Upper layer PV gradient /β	2	2.81 to 3.79
	Thermal damping rate	r	0 to 0.5
	Enstrophy input rate	η	0.01 to 100
	Heton PV amplitude	q heton	1.257 to 1508
	Heton radius	R	0.05

  Upper layer (left figure) and lower layer (right figure) zonally-averaged zonal velocity ūi(y, t) at two different times t = 3000 blue and t = 9400 orange.

				t = 9400											t = 9400
	0.4			t = 3000						0.4						t = 3000
	0.2										0.2
	0.0										0.0
	0.2										0.2
	0.4										0.4
	0.4	0.2	0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4	0.4	0.2	0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4
	Figure 6:

Table 2 :

 2 Numerical and physical parameters for the set of experiments conducted in this study.

	Nx	cρ	1	η	q1 × k d /β	r	other
	512	2	-0.5	0.1	1	0, 10 -4 , 10 -3
						10 -2 , 0.1, 0.5
	512	2	-0.1	0.1	1	0, 10 -3 , 10 -2
						0.1, 0.5
	512	0.2	-0.1	0.1	1	0, 10 -3 , 10 -2
						0.1, 0.5
	512	0.2	-0.1	0.01, 1	1	0.01
	512	2	-0.1	0.01, 1	1	0.01
		0, 2	-0.5	10, 100 10, 30, 100, 300	0.01
	256	2	-0.1	100	100, 300, 1200		Heton radius
	256	2	-0.1	100	300	0.001	R = 0.025 rather
	256	0	-0.1	100	100	0.01	than 0.05
	256	0	-0.5	100	100	0.01	h1 = 0.9
	256	0	-0.5	100	100	0.01	h1 = 0.95
	256	5	-0.5	100	100	0.01
	256	2	2.0	1	1	0.001	opposite
	256	0	1.5	1	1	0.001	vertical shear
	256	2 -0.5, -0.4, -0.3	0.1	1	0.01
			-0.2, -0.1
	128, 512 2	-0.5	0.1	1	0.01
	256	0	1	100	100	0.01	with and without
							topography

  11 Ψ 1 -ū1 B 12 Ψ 2 = c d 2 Ψ 1 dy 2 -B 11 Ψ 1 -B 12 Ψ 2 21 Ψ 1 -ū2 B 22 Ψ 2 = c d 2 Ψ 2 dy 2 -B 21 Ψ 1 -B 22 Ψ 2 .

	ū1 -ū1 B ū2 d 2 Ψ 1 dy 2 + dq 1 dy d 2 Ψ 2 dy 2 + dq 2 dy -ū2 B (A 26)

† Email address for correspondence: tj30@st-andrews.ac.uk

† A more representative measure of L d is 2π/ kd = π/10, yet this is still much smaller than Ly.

Acknowledgements

The authors wish to acknowledge the LabEx OSU G@2020 and the University of St Andrews for the funding of TJ.

Appendix A. Layer model decomposition and related analyses

A.1. Vertical modes

The inversion of (2.6) for the layer streamfunctions ψ 1 and ψ 2 is accomplished through a projection onto vertical modes, i.e.

where a hat indicates a vertical mode, and the c ij are the projection coefficients, determined as follows. Using the same projection for the PV anomaly q i -βy in (2.6), we obtain

for each vertical mode (j = 1 or 2). The objective is to make the terms involving c j1 and c j2 equal to -γ j k2 d ψj for each j, since then we have simple Helmholtz-type equations to invert for ψj , namely

where k 2 dj ≡ γ j k2 d is the squared deformation wavenumber for mode j. Equating then the terms above with -γ j k2 d ψj and re-arranging, we obtain

As this must be true for all ψ 1 (x, y, t) and ψ 2 (x, y, t), the constant coefficients must vanish:

The only nontrivial solution results by taking the determinant of this linear system to be zero, leading to

which gives the prefactor γ j in k 2 dj ≡ γ j k2 d :

We have taken the minus sign in front of the square root for the lower mode, since it is purely 'barotropic' with γ 1 = 0 in the Boussinesq limit α → 1. The other mode is called 'baroclinic' and always has the higher deformation wavenumber. Note:

Regarding the projection coefficients c j1 and c j2 , by convention we take c 11 + c 12 = 1 for the 'barotropic' mode j = 1. This, along with (A 5) gives

For the 'baroclinic' mode j = 2, we take c 22 = 1. Again using (A 5), we find

The inverse of these coefficients are needed to express modes in terms of layers, i.e. as

where ∆ = c 11 c 22 -c 12 c 21 .

A.2. Vertical shear

A vertical shear (a uniform flow ūi in each layer i) is imposed by taking the lower layer mean PV to be q1 = -1 βy, specifying 1 , and requiring the mass-integrated momentum to be zero, or h 1 ū1 +αh 2 ū2 = 0. The latter requirement determines the constant 2 in the associated upper layer PV q1 = 2 βy, as well as the mean flows ū1 and ū2 , as explained next.

We start with (2.6) for the mean flow, for which ψi = -ū i (y -y c ), where y c is the centreline of the channel. For this flow, the relative vorticity ∇ 2 ψi is zero. After cancelling the common y factor and some re-arrangement, (2.6) implies

We can solve these for the mean flows ū1 and ū2 in terms of 1 and 2 , but we also want to obtain 2 from the condition of zero momentum, h 1 ū1 + αh 2 ū2 = 0. Omitting the details, this gives

The mean velocities ûj associated with each mode j are

from which we obtain the mean layer velocities by projection:

Initially, the flow is in thermal equilibrium, with uniform velocity ūi in each layer i. Associated with this flow, by thermal wind balance, the layer interfaces adopt a linear profile δ eq,i ∝ y -y c . In terms of the scaled displacements δeq,i = f 0 δ eq,i /H, upon using their definitions in terms of streamfunctions (2.9), we find

These are the profiles to which the interfaces are relaxed to subsequently, at a specified rate r.

A.4. Linear stability

We next provide brief details of the linear stability analysis used in the main body of the paper to calculate the maximum growth rate. We start with the instantaneous zonallyaveraged PV qi and zonal velocity ūi in each layer i = 1, 2. In this subsection, these profiles are general functions of y for each sampled time t. In the numerical code, they are provided at each y grid point (typically 129 or 257 grid points, including the edges). The PV gradient dq i /dy, needed in the analysis below, is computed using centred differences in y. Notably, only the interior grid point values are needed, since the perturbation satisfies homogeneous boundary conditions (required to ensure zero meridional velocity at each boundary).

In the linear stability analysis, the base flow is considered to be both x and tindependent. Hence, we can consider infinitesimal perturbations of the form

where k x is a prescribed zonal wavenumber (necessarily integer for a domain of length L x = 2π), σ is the frequency (the imaginary part of which, if non-zero, is the growth rate), while Q(y) and Ψ (y) are the perturbation amplitudes (eigenfunctions) to be determined along with σ (the eigenvalue).

Plugging the above forms into the conservative form of the governing equations (2.1) -with F i = 0 -and in the PV inversion relations (2.6), then linearising, and finally cancelling the common e i(kxx-σt) } factor, we obtain the following eigen-system for the

Appendix B. Evolution variations

To illustrate the variability occurring between different simulations of the same flow, two other nearly identical simulations have been conducted. All parameters were kept the same as in the characteristic simulation previously illustrated, only the inversion grid resolution was halved to 256 × 129. Additionally, a different random number seed was used in each simulation.

The two simulations exhibit striking differences. In the first case, see figure 16, we see oscillations between turbulent and quiescent phases, sometimes with a shift of the jets and homogeneous regions, while retaining the jet spacing. The flow alternates between the two different quiescent phases encountered in the characteristic simulation. The 'centred' state exhibits jets slowly coming together toward the center of the domain. By contrast, the 'shifted' state does not exhibit any evident jet drift and moreover appears to be much shorter lived. Arguably, the jet configuration of this state is less robust than that in the centred state, though both the stability analysis of the zonal flow and the APE spectrum show no major differences between the states.

In the second simulation, see figure 17, the shifted quiescent state develops after the initial baroclinic instability and dominates the first third of the simulation. Turbulent events throughout this period collapse back to this state except for the last one. This prominent turbulent event instead collapses to a new configuration resembling the centred state seen over much of the characteristic simulation. This state then dominates the remainder of the evolution.


Adrien Chan-Hon-Tong

Nicolas Audebert

Pointing is sufficient for remote sensing images

Object detection in remote sensing images with center only

Introduction 1.Segment before detect

Since the publication of [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], deep learning is more and more becoming a common tool. High quality deep learning engines like PYTORCH and TENSORFLOW offer simple and efficient way to run a network on GPU and/or on the cloud. In this context, given a specific computer vision problem, designing an end to end deep networks allows to simplify the code and to get an automatic GPU acceleration.

Thus, end to end networks are useful even with slightly lower accuracy/quality that the state of the art for the given problem. As an example, [START_REF] Chen | Fast image processing with fully-convolutional networks[END_REF] learn end to end network for low level computer vision operator accelerating these low level operation thank to automatic GPU implementation.

Now, some problems are more likely to accept an end to end network. For example, object detection does not trivially accept one. Object detection is one of the oldest computer vision problem, and consists in producing a list of bounding boxes of all class instances from an image (see figure 1). So, the output of a detector is highly unstructured and inhomogeneous: unstructured because expected to be a list whose size is unknown, inhomogeneous because numerical representation of bounding boxes lead either to mix coordinate and scale or either to accept numerical pair of coordinates with large variances. Indeed, state of the art deep networks designed for object detections introduce a lot of non standard operations: pre computation of candidates [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF], ROI pooling trick [START_REF] Shaoqing Ren | Faster R-CNN: Towards real-time object detection with region proposal networks[END_REF] or anchor trick [START_REF] Redmon | Yolo9000: better, faster[END_REF].

Figure 1: Illustration of object detection: putting bounding boxes on object in the image

Inversely, semantic segmentation is as suitable as classification for end to end deep networks. Semantic segmentation consists in predicting a class for each pixel (or pixel cells of an image). So, output of semantic segmentation pipeline is a known size map with one channel per class (corresponding value can be seen as likelihood of belonging to each class). This is highly structured and homogeneous, allowing very straightforward end to end deep network like [START_REF] Badrinarayanan | Segnet a deep convolutional encoder decoder architecture for robust semantic pixelwise labelling[END_REF][START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] More interesting, in case of remote sensing images, an accurate segmentation mask can be converted into a detection output by considering the connected components: in a remote sensing images, there is only few perspective effects, so an object will rarely occulted an other -this way, all instances will be in different connected components. This one scale property is also the case for lot of medical images including cytology, but in medical images biological objects can merge into complex biological tissues. Inversely, in remote sensing image objects do not overlap each other (e.g. cars do not overlap even on a parking).

Thus, in remote sensing images, segment before detect is a rising paradigm whose an example can be found in [START_REF] Audebert | Segmentbefore-detect: Vehicle detection and classification through semantic segmentation of aerial images[END_REF]. We can also stress that, even on natural image, connexion between detection and segmentation leads to state of the art methods [START_REF] He | mask r-cnn[END_REF].

Ground truth issue

Now, producing a semantic segmentation ground truth is so time consuming that large semantic segmentation datasets may only appear for market Segmenting is much more time consuming than just pointing objects i.e. giving the object centres (see figure 3).

In fact, even providing bounding boxes is time consuming: validating a box is twice more fast than drawing it in [START_REF] Russakovsky | Best of both worlds: human-machine collaboration for object annotation[END_REF]. More, [START_REF] Russakovsky | Best of both worlds: human-machine collaboration for object annotation[END_REF] does not take into account that drawing boxes leads to much more complex human computer interface than clicking points. Base on private little experiments with disconnected softwares, we currently think that drawing boxes is close to 20 times longer than just pointing objects.

With this order of magnitude in mind, we argue that it is relevant to aim to learn deep segmentation pipeline from pointing ground truth (centres of objects) instead of detection ground truth (boxes of objects) or, off course, segmentation ground truth (masks of objects).

Alternatively, there are lot of works aiming to reduce the need of human annotations: 0 shot learning [START_REF] Bucher | Improving semantic embedding consistency by metric learning for zero-shot classiffication[END_REF], self supervised training [START_REF] Doersch | Multi-task self-supervised visual learning[END_REF] or interactive verification instead of annotation [START_REF] Russakovsky | Best of both worlds: human-machine collaboration for object annotation[END_REF][START_REF] Dim | We don't need no bounding-boxes: Training object class detectors using only human verification[END_REF]. However, current non supervised trainings have today lower accuracy than supervised one. Also, interactive verification needs a very large software infrastructure to become possible and [START_REF] Dim | We don't need no bounding-boxes: Training object class detectors using only human verification[END_REF] currently does not really demonstrate that it decreases the human time needed to annotate. It just shows that it decreases the complexity of what is asked to human.

Inversely, our contribution is to offer a (quite simple) trick to learn deep segmentation pipeline with pointing dataset only, in remote sensing context.

This way, we build on the top of the [START_REF] Audebert | Segmentbefore-detect: Vehicle detection and classification through semantic segmentation of aerial images[END_REF] but removing the need for an expensive semantic segmentation ground truth.

These tricks are presented in the next section. Then, in section 3, we present datasets and results, before conclusion of section 4. Baseline: Natural segmentation training for convolutional neural networks consists to compare the output mask to the ground truth mask, like if each pixel was a classification data, averaging the loss over all the pixels [START_REF] Badrinarayanan | Segnet a deep convolutional encoder decoder architecture for robust semantic pixelwise labelling[END_REF][START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF].

So, typically, in segmentation, the output tensor p has a size of B × C × W × H (one likelihood per class and per pixel) and ground truth has a size B × W × H (with values from {0, ..., C -1}). The loss can be simply computed by seeing both those tensors like classification tensors with batch size being

B × W × H instead of B: loss = B,W,H b=1,w=1,h=1 - p[y[b]][h][w] + log( C c=1 exp(p[c][h][w])).
Now, such natural way is completely not able to handle pointing ground truth converted in mask by setting one foreground pixel per cen-tres. Indeed, such centre mask is too much unbalanced.

This natural way to learn a segmentation deep pipeline will be called multi classes or binary or raw training in the following of this paper, depending on the ground truth used to train (original mask for multi classes, binary mask for binary and centre mask for raw ).

Extending centres: The first trick to be able to learn from pointing ground truth is simply to fill a square of predefined size centred on each object centre. Then, training is done as usual.

The main problem with this trick is that pipeline will have to learn from noisy ground truth: background pixel close to a car will be tagged as foreground. But, balance will be globally the same than in the binary mask while needing only object centres.

Removing border: The second trick is highly inspired from [START_REF] Jiu | Human body part estimation from depth images via spatiallyconstrained deep learning[END_REF]. [START_REF] Jiu | Human body part estimation from depth images via spatiallyconstrained deep learning[END_REF] shows the relevancy of removing pixels close to a boundary of the ground truth segmentation. Precisely, [START_REF] Jiu | Human body part estimation from depth images via spatiallyconstrained deep learning[END_REF] offer to remove pixel from boundaries during training, but, show that it increases performance even during testing while boundaries are kept at this stage. This result is not trivial as breaking the symmetry between training and testing (remove border from train but not test) is usually considered as a bad practice. Authors from [START_REF] Jiu | Human body part estimation from depth images via spatiallyconstrained deep learning[END_REF] argues that by removing ambiguous pixel, segmentation pipeline is able to learn better and thus removing large errors localized far from border by not focusing on small border errors.

Unsupervised restoration: graph cut [START_REF] Xu | Object segmentation using graph cuts based active contours[END_REF] like algorithms can be applied to try to restore the mask from the known points (centres for foreground and pixel far from centres for background) by taking advantage of the spatial consistency of the segmentation.

Merging all of them: we presented 3 tricks (removing border, extending the centres and using a graph cut to try to restore the mask), but we can combine all of them. This strategy will be called coarse training in the following. Formally, the loss is: loss = b,w,h∈allowed

-p[ỹ[b]][h][w]+log( C c=1 exp(p[c][h][w]
)) with ỹ being the ground truth extended (with a constant kernel or a graph cut) and allowed the set of pixels far from border in ỹ. In practice, removing some pixels is done by adding a dumb class with class weight equal 0.

Keep the pooling: The idea of making the learning easier to increase performances from [START_REF] Jiu | Human body part estimation from depth images via spatiallyconstrained deep learning[END_REF] leads to an other trick: keeping a pooled segmentation. This trick consists in decomposing the image in cells corresponding to spatial pixel of the last layer of the network after all pooling layers. Then, a cell is considered as foreground as soon as it contain an object centre.

Currently, since the apparition of VGG [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF], large amount of works have been published to try to restore spatial dimension (dilated convolution, unpooling with index [START_REF] Badrinarayanan | Segnet a deep convolutional encoder decoder architecture for robust semantic pixelwise labelling[END_REF], ...). This spatial restoration is natural for pure segmentation. Now, in pointing context, restoring spatial resolution is not the goal. This is why we offer this simple trick which will be shown to be quite efficient. Indeed, we argue that it provides 2 advantages: it reduces the unbalance by dividing background instance by the cell size, and, it allows to have a more network friendly training as expected output is very close to the network structure.

This last trick will be called pooled training and consists formally to use the following loss:

loss = B,W ,H b=1,w=1,h=1 -p[y [b]][h][w] + log( C c=1 exp(p[c][h][w]))
with W , H being the size of the image after pooling (typically H = H 16 and y is the pooled ground truth.

Experiments

Datasets

We use four datasets in our benchmark: the data fusion contest 2015 [START_REF] Lagrange | Benchmarking classification of earthobservation data: from learning explicit features to convolutional networks[END_REF] (DFC2015), the ISPRS POSTDAM dataset [START_REF] Rottensteiner | The isprs benchmark on urban object classification and 3d building reconstruction[END_REF] (POSTDAM), the VEDAI dataset [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery : A small target detection benchmark[END_REF], and we also rely on SACLAY1 which is a private dataset (but planned to be released).

DFC2015 dataset2 provides 6 10000x10000 ortho images at very high resolution (5cm). [START_REF] Lagrange | Benchmarking classification of earthobservation data: from learning explicit features to convolutional networks[END_REF] provides a ground truth for pixelwise semantic segmentation in 7 classes including vegetation, building, road, car, boat. Here, we downscale image to 5120x5120 (thus resolution is around 10cm). We split train and test according to [START_REF] Lagrange | Benchmarking classification of earthobservation data: from learning explicit features to convolutional networks[END_REF].

POSTDAM dataset3 is very similar to DFC2015. It provides 38 ortho images at high resolution (between 5 and 10 cm). As for DFC2015, we also downscale each 6000x6000 image to 4096x4096.

VEDAI dataset is a detection dataset. Ground truth is composed of bounding boxes of vehicles. The dataset is composed of around 100 1024x1024 orthoimages at 12.5 cm of resolution. We do not use the label associated to the bounding boxes.

Finally, SACLAY is a private dataset formed with IGN ortho images, annotated for car pointing. The dataset contains around 20 5000x5000 pixel image at 20cm of resolution (we crop the 4608x4608 corresponding image).

Dataset conversion

All datasets are converted into vehicle pointing datasets (ground truth is a set of centres of vehicle).

For VEDAI, centres are either extracted from ground truth boxes. All boxes are considered as vehicles (car, truck, farming vehicle, plane: all are considered as vehicle).

For POSTDAM, ISPRS, we extracted connected component of vehicles and compute the centres. Thus, we lost the multi class aspect: we only keep background (all except vehicle vs foreground vehicle).

This label loss is important as it decreases performance. We compare detection results achieved by multi classes segmentation [START_REF] Audebert | Segmentbefore-detect: Vehicle detection and classification through semantic segmentation of aerial images[END_REF] versus binary segmentation ( [START_REF] Audebert | Segmentbefore-detect: Vehicle detection and classification through semantic segmentation of aerial images[END_REF] but from binary mask). Performance drops from 75% of Gscore to 44% on POSTDAM.

We will show in the later that training with pointing may lead to the same performance than training with binary segmentation, but both are lower than training with multi class segmentation. The question about why having multi classes helps the classifier is out of the scope of this paper, but, such kind of results are not new (see [START_REF] Ristin | From categories to subcategories: large-scale image classification with partial class label refinement[END_REF]).

For the coarse trick, we dilate centres with a kernel ranging from 15 to 31. We remove border with a kernel ranging from 9 to 17. Graph cut is done on each centre and not globally on the hole image without specific tuning. Currently, results are globally close when moving these pre processing parameters. Thus, for more reproducible results, we offer to use the real mask eroded by a kernel of size 11 as representative of coarse training.

At test time, each pipelines produces a segmentation mask which is post processed to produce a set of points. Post processing is erosion and dilatation (with kernel ranging from 7 to 11) and extraction of the centres of connected components.

Predicted points and ground truth points are matched like boxes in detection (criterion for matching is a euclidean distance less than 25 pixels). This leads to correct matches, false alarms and miss detections, and then to precision recall. Then, we form the Gscore (product of precision and recall).

Result

All pipelines are trained for segmentation in a UNET and/or FCN fashion [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] except for the pooled training which is only a VGG.

The main result of this experiment is that pooled training performances are only sightly lower than binary ones while pooled training needs only centres and binary needs masks. Also, pooled training outperforms other coarse training. And raw training just leads to nothing (raw is training from a mask while only 1 pixel tagged as foreground per object). Table 1 presents Gscore of all trainings.

An other results is that when matching pooled training mask against the original mask, we reach only slightly lower pixelwise performance than the original binary training. More precisely, we obtain much more coarse mask from spatial point of view, but, we catch more foreground connected component and have less false foreground connected component. Results are in table 2.

Conclusion

We offer a simple trick to learn segmentation pipelines from pointing datasets (ground truth is just the list of object centres). This trick is to virtually forward the pointing mask through the pooling layers of the network in order to reduce unbalance and to offer the network a expected output closer from this internal structure. Performance are in Gscore in % (here we count the number of pixels not the number of objects).

trick reaches interesting results: it even reaches only slightly lower results than pipeline trained from mask (while being learn from centres). These results invite to use pointing instead of boxes or masks when running out of time/fund to make a dataset ground truth.

Figure 2 :

 2 Figure 2: Illustration of semantic segmentation: give a class to each pixel

Figure 3 :

 3 Figure 3: Illustration of pointing dataset: ground truth is the set of object centres

Table 1 :

 1 The pipeline trained with our Gscore of training strategies under pointing metric for fcn architecture.Results are Gscore in %. Datasets are converted into pointing dataset: we keep the centres of each bounding boxes or connected components. raw consist to learn from mask containing one foreground pixel per objects. coarse consist to learn from mask where centres are extended and border as unknown. pooled consist to learn grid quantified mask.binary consist to learn from original binarized mask (available only for POSTDAM and DFC2015).

	dataset	raw coarse pooled binary
	DFC2015	0	20	26	31
	POSTDAM	0	30	41	44
	VEDAI	0	4	19
	SACLAY	0	3	26

Table 2 :

 2 Evaluation of segmentation mask learnt from segmentation or centres.

We thank Adèle Koeniguer for providing the ground truth for SACLAY

www.grss-ieee.org/.../data-fusion/2015-ieee-grss-data-fusion-contest/

http://www2.isprs.org/commissions/comm3/wg4/2d-sem-label-potsdam.html

Acknoledgment

We thank Adèle Koeniguer for providing the ground truth of SACLAY.