égayaient souvent mes journées à GTT, ainsi que pour ton soutien en fin de thèse. Jegan, merci pour tout le temps que tu m'as consacré durant mon passage à GTT, mais surtout pour ta bienveillance, ton estime, la confiance que tu m'as accordée, et tous les précieux enseignements que j'ai retirés de nos conversations ! Mohamed, merci pour ton honnêteté, ta sympathie, ta bonne humeur et ta curiosité. Enfin, Zahra, merci pour ta gentillesse, ta finesse, ton ouverture d'esprit et toutes tes attentions à mon égard. J'espère m'en être montré digne et en particulier t'avoir aidée et soutenue comme tu le méritais au cours de ces derniers mois.

Sur un tout autre plan, je remercie chaleureusement mes amis de la section « Musculation » du club sportif de l'Ecole polytechnique : Christophe, Marie, Marwane, Paul C., et Yannick, pour avoir contribué, chacun à votre façon, à me distraire aussi bien durant nos entraînements qu'en dehors. Un grand merci à Christophe et Marie pour les (longs !) entraînements du week-end ! Christophe, je garde un souvenir particulier de ton état d'esprit authentique, et de cet inébranlable enthousiasme que l'on partage à toujours tenter de repousser nos limites sans crainte des éventuels échecs. Marie, je salue ta gentillesse et ta patience ! Marwane, merci pour tes conseils avisés lorsqu'il s'agissait de parfaire ma technique sur certains exercices ! Un merci spécial à Paul C. pour tous les entraînements en extérieur, aussi bien à Gif-sur-Yvette qu'à Massy, par grand beau temps comme sous une pluie battante, mais avec le sourire quelles qu'aient pu être les circonstances ! Je pense que les poteaux des terrains de rugby de Gif se souviennent encore de nos passages improvisés aux anneaux ! J'espère que l'on parviendra un jour tous les deux à maîtriser ce fameux Graal que constitue la « full planche », voire, soyons ambitieux (il faut l'être), la « maltese » ! Enfin, Yannick, un grand merci, tu sais pourquoi. Je mentionne tout de même que je suis désolé pour les « quelques » grognements que j'ai parfois laissé échapper au cours de mes entraînements, ainsi que pour les quelques séismes occasionnés par mes séances de soulevé de terre...

Merci !

Je remercie tout d'abord mon directeur de thèse, Paul Billant, pour de très nombreuses raisons qu'il serait trop long d'énumérer ici. Je citerai néanmoins les principales : Paul, merci pour ta disponibilité sans failles et ta bienveillance en toutes circonstances, y compris dans les moments les plus délicats de cette thèse qui m'aura donné bien du fil à retordre ! Je retiens également la liberté que tu m'as laissée tout au long de la thèse, une liberté telle qu'elle m'a permis de continuer à m'épanouir dans mes passions et en particulier d'assurer la dizaine d'heures d'entraînement hebdomadaires dont j'ai tant besoin pour garder mon discernement et mon équilibre au quotidien. Enfin, comment ne pas mentionner le travail colossal abattu lors de cette étape très éprouvante qu'est la rédaction ? Si toutefois je ne devais retenir qu'une chose parmi tout ce que cette thèse a exacerbé en moi, ce serait sans aucun doute la ténacité, car c'est là ce qui a, en fin de compte, le plus d'importance, puisque cela continuera de me servir quotidiennement à l'avenir, aussi bien dans le cadre professionnel que sur un plan plus personnel. Je te remercie donc tout simplement pour m'avoir guidé patiemment et avec subtilité tout au long de cette longue mais inestimable progression. J'en profite pour saluer le petit comité des numériciens, à savoir Ernesto, Eunok, Gaétan, Guillaume, Jérémy, Léopold, Lutz, Mathieu B., Océane, et Tristan. Je dois cependant des remerciements spéciaux à Léopold, Océane et Tristan, pour toutes les choses que nous avons partagées ces dernières années, mais aussi pour votre gentillesse, votre intelligence, et votre discernement. Merci également à Adrien Saint-S., Camille, Caroline C., Caroline F., Emmanuel de L., Jean-Philippe et Vincent pour votre sympathie. Merci par ailleurs aux secrétaires, Delphine, Magali, Sandrine et Thérèse, pour faciliter toutes nos démarches administratives avec le sourire, ainsi qu'à Daniel et Toai pour vous assurer que les problèmes d'ordre numérique trouvent des solutions rapides et efficaces. Je me dois bien évidemment de remercier les membres de mon jury de thèse, et notamment mes deux rapporteurs, pour avoir accepté d'évaluer mes travaux, en particulier au vu des efforts et de l'attention que les modélisations théoriques du quatrième chapitre requièrent, mais aussi pour la curiosité et l'intérêt que vous avez tous témoigné à l'égard de mes travaux.

Ma thèse, bien qu'académique, ne s'est toutefois pas bornée au laboratoire puisque j'ai également eu la chance d'effectuer une activité complémentaire en entreprise en parallèle de celle-ci. Parmi toutes les personnes que j'ai pu côtoyer dans ce cadre-là, j'adresse des remerciements tout particuliers à Antoine, Jegan, Mohamed et Zahra. Une multitude de raisons me viennent à l'esprit pour vous exprimer ma profonde gratitude. Antoine, merci pour ton franc-parler, ton humanité et pour toutes tes histoires qui Chapter 1

Introduction

The context of the present dissertation is that of Geophysical Fluid Dynamics (GFD). GFD is a branch of fluid mechanics that is concerned with the dynamics of atmospheric and oceanic flows such as cyclones, jets, or streams for example.

Geophysical flows 1.Structure and classification

The atmosphere and the oceans both consist of shallow layers of fluid: their thickness is small relative to their horizontal extent. For example, hurricanes generally exhibit a diameter of the order of 100 km. In the vertical direction however, they remain confined to the bottom layer of the atmosphere, i.e. the troposphere whose thickness is about 10 km. Therefore, their aspect ratio is of the order of 0.1 and is thus well below unity. As for oceans, their maximum depth is reached in the Mariana Trench, and is also 10 km approximately. Since they extend over more than hundreds of kilometers in the horizontal, it is clear that their aspect ratio is small as well.

A direct consequence is that geophysical flows are classified in terms of their horizontal length scale. Hurricanes are planetary scale flows, whereas tornadoes for instance are much more local eddies. This variety of scales is illustrated in figure 1.1, which is a satellite image of the Oyashio Current that displays multiple eddies of various diameters.

The notion of coherent structures is also quite frequently mentioned in GFD. In most cases, these coherent structures refer to vortices, because they are ubiquitous components of geophysical flows as indirectly suggested previously.

Central physical effects

As underlined in many books about GFD [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF][START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF][START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF][START_REF] Cushman-Roisin | Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects[END_REF], the dynamics of atmospheric and oceanic flows is constrained by two physical effects, namely the ambient stable stratification and the rotation of the Earth.

The stratification is the result of density variations along the vertical. These changes in the density arise from fluctuating temperature and pressure in the atmosphere, and from varying temperature and salinity in the oceans. In general, the density decreases upwards, resulting in a configuration where lighter layers of fluid are located above 9 Figure 1.1: Photograph taken from space in March, 1992, of eddies in the Oyashio Current. Sea-ice allows visualizing the flow. It is reproduced from [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF]. the heavier fluid. Thus, the stratification is usually stable. If we consider a stratified flow at rest, then the isopycnals -i.e. the lines of constant density -are rigorously horizontal. Under the effect of any disturbance, the parcels of fluid will depart from their equilibrium position by oscillating vertically around this mean position, at a typical frequency called Brunt-Väisälä frequency. These oscillations are due to the fact that the stratification acts like a spring, whose stiffness would be the equivalent of the intensity of the stratification. A strong stable stratification therefore tends to inhibit vertical motions of large amplitude.

The dynamics of geophysical flows is also affected by the Earth's rotation through the existence of a Coriolis force. How influential stratification and rotation are actually depends on the length scale of the flow. Indeed, they are both significant at the largest geophysical scales. In contrast, there exist length scales where stratification prevails over the Earth's rotation, namely the atmospheric mesoscales, comprised between 1 km and 100 km globally, and the oceanic sub-mesoscales that range from 10 m to 10 km approximately. The present dissertation deals with the stratified regime only.

Notion of turbulence

Geophysical flows exhibit a large length scale even in the purely stratified regime. Thus, their Reynolds number is very large as well. In other words, inertia effects largely dominate the viscous dissipation. Such conditions are conducive to the generation of turbulence [START_REF] Lesieur | Turbulence in Fluids[END_REF], which is why we briefly discuss this concept in the following.

Global features

Turbulent motions share a variety of typical features. First, they are disordered: the motion of the fluid does not exhibit any clear trend because the fluid parcels, since they move erratically, cannot be clearly distinguished from each other. This is illustrated in figure 1.2, which shows a turbulent motion generated by a grid. In the upstream region at the left of the grid, the incident flow is uniform and laminar: the streamlines can be clearly distinguished. In contrast, in the downstream region, the streamlines are tightened. There is a first sub-region where they still can be identified, but the fact that they are closer to each other implies the existence of shear layers. Their destabilisation ultimately leads to a transition towards turbulence, after which the flow no longer exhibits any clear structure.

Then, turbulence is a chaotic phenomenon. For example, it is very sensitive to the conditions imposed at the onset of the motion. More generally, the local evolution of any turbulent flow is not deterministic and is therefore not predictable from the classical Navier-Stokes equations. For this reason, the dynamics is preferentially characterized on the basis of statistical averages.

Finally, turbulent flows are generally seen as complex flows made of a flurry of vortices of various sizes that interact with each other. This is at the origin of energy and enstrophy -the vorticity squared -transfers. The dynamics of these exchanges will be discussed in the following, because it actually depends on the nature of the turbulent motions.

Three-and two-dimensional turbulence Three-dimensional isotropic homogeneous turbulence

Isotropic homogeneous turbulence is the model case of reference. It has been considered at the origin by [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluids for very large Reynolds numbers[END_REF], who derived a number of theoretical results regarding its spectral properties. In particular, Kolmogorov has demonstrated the famous law which states that the kinetic energy spectrum E(k) varies like k -5/3 , where k ∼ 1/l is the wavenumber associated with the length scale l of the turbulent motion.

It is worth stressing that this spectrum is in fact very robust in the sense that it accounts for the dynamics of a wide variety of turbulent flows. Notably, stratified turbulence, which will be extensively discussed further, is characterized by this spectrum, although it is strongly anisotropic in sharp contrast with isotropic homogeneous turbulence.

In isotropic homogeneous turbulence, the energy is transferred from the largest towards the finest scales through a so-called direct energy cascade. The acting process is referred to as vortex stretching. Basically, the largest structures provide some part of their energy to smaller structures through non-linear inertial mechanisms and so on until the so-called Kolmogorov scale, where the energy is ultimately dissipated by viscosity. Indeed, this scale is defined as the scale at which the Reynolds number is of order unity. To put it another way, inertial effects are balanced by viscous effects below this threshold. Statistically, the turbulent motion is evenly distributed along the three dimensions.

Two-dimensional turbulence

Besides three-dimensional isotropic turbulence, another essential concept in turbulence is that of two-dimensional turbulence [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF][START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF][START_REF] Lilly | Numerical simulation of two-dimensional turbulence[END_REF], first for its strong theoretical interest, and then because geophysical turbulence resembles bidimensional turbulence in a certain regime [START_REF] Charney | Geostrophic turbulence[END_REF].

From a theoretical point of view, this kind of turbulence clearly departs from its three-dimensional counterpart in that the stretching of the vorticity is no longer possible geometrically. This is at the origin of the conservation of the vorticity and enstrophy in the inviscid limit. As a result, the energy transfers are tremendously affected, since a forward cascade (towards small scales as in three-dimensional turbulence) can no longer settle. A simple demonstration is provided by [START_REF] Augier | Turbulence en milieu stratifié, étude des mécanismes de la cascade[END_REF], which is essentially based on the fact that the enstrophy is proportional to the energy in spectral space. It therefore turns out that two-dimensional turbulence is characterized by an inverse energy cascade, namely from the smallest towards the largest scales.

Geophysical turbulence

Global issues

The early work by [START_REF] Lumley | The spectrum of nearly inertial turbulence in a stably stratified fluid[END_REF] dealt with the nature of stratified turbulence. This investigation as well as others [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF][START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF][START_REF] Lilly | Numerical simulation of two-dimensional turbulence[END_REF] participated in motivating a flurry of subsequent studies regarding the dynamics of turbulence in the atmosphere and the oceans.

In particular, [START_REF] Charney | Geostrophic turbulence[END_REF] focused on the dynamics of quasi-geostrophic turbulence, which takes place in atmospheric and oceanic flows whose horizontal length scale is large enough to ensure that the stratification as well as the Earth's rotation are strongly influential in determining the motion of the flow. In his paper, Charney underlined that quasi-geostrophic turbulence, although not rigorously bidimensional, resembles two-dimensional turbulence, hence the importance of this model case. This finding has naturally suggested that the energy cascades towards large scales in a certain range of length scales in the atmosphere and the oceans, thus asking the question of how energy can be dissipated at smaller scales. A global concern has therefore been to identify the nature and mechanisms of the energy cascades in geophysical media, i.e. to understand the physical processes that characterize purely stratified turbulence on one hand and stratified-rotating turbulence on the other hand. [START_REF] Gage | Evidence for a k -5/3 law inertial range in mesoscale two-dimensional turbulence[END_REF] argued that the atmospheric mesoscale kinetic energy spectrum results from an inverse energy cascade through two-dimensional turbulence as described by [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF]. [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] supported this analysis. Atmospheric mesoscales correspond to length scales that are essentially influenced by a strong stratification. It has been conjectured that stably stratified flows exhibit a two-dimensional dynamics when the stratification is strong, notably because vertical motions are then inhibited [START_REF] Riley | Direct numerical simulations of homogeneous turbulence in density-stratified fluids[END_REF][START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF][START_REF] Lilly | Stratified turbulence in the atmospheric mesoscales[END_REF][START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]. Thus, stratified turbulence is strongly anisotropic, and mainly concentrated in the horizontal direction. The latter result effectively tends to indicate that the dynamics of the energy and enstrophy transfers at atmospheric mesoscales can be properly explained on the basis of two-dimensional turbulence.

Stratified turbulence

Nevertheless, considering the two-dimensional limit is valid only if vertical gradients are not too intense. For this reason, several authors have questioned the hypothesis that strongly stratified flows behave as nearly two-dimensional flows. Indeed, Herring and Métais (1989); [START_REF] Métais | Numerical simulations of freely evolving turbulence in stably stratified fluids[END_REF]; [START_REF] Park | Turbulent mixing in stratified fluids: layer formation and energetics[END_REF]; [START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably stratified turbulence. Part 2. Large-scale and small-scale anisotropy[END_REF]; [START_REF] Waite | Stratified turbulence dominated by vortical motion[END_REF]; [START_REF] Lindborg | Stratified turbulence forced in rotational and divergent modes[END_REF] report that strongly stratified turbulence is made of thin horizontal layers, hence implying the existence of strong vertical variations. This layering is clearly visible in figure 1.3, which shows vertical cross-sections of typical turbulent velocity fields in presence of a strong stratification. [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] clearly confirmed the idea of strong vertical variations. This result has therefore motivated additional studies in order to clarify the nature of the energy cascade that is characteristic of stratified turbulence. Notably, [START_REF] Lindborg | Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence?[END_REF] and [START_REF] Cho | Horizontal velocity structure functions in the upper troposphere and lower stratosphere: 1. observations[END_REF] invalidated Gage's (Gage, 1979) and Lilly's proposal [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]) that the mesoscale range of the atmosphere is the seat of an inverse energy cascade due to bidimensional turbulence. Indeed, by studying structure functions, [START_REF] Lindborg | Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence?[END_REF] and [START_REF] Cho | Horizontal velocity structure functions in the upper troposphere and lower stratosphere: 1. observations[END_REF] have demonstrated that the energy is actually transferred downscale. Evidence supporting this theory have been brought in more recent papers such as [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF], [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], and [START_REF] Riley | Stratified turbulence: A possible interpretation of some geophysical turbulence measurements[END_REF]. Besides, this result has been extended in the presence of a background rotation since [START_REF] Lindborg | The effect of rotation on the mesoscale energy cascade in the free atmosphere[END_REF] has shown that this energy cascade can actually prevail as long as the Rossby number, which is the ratio of inertial to Coriolis effects, is larger than ∼ 0.1. Figure 1.3: Vertical cross-sections of horizontal velocity fields taken from Direct Numerical Simulations of stratified turbulence. The snapshots are reproduced from [START_REF] Lindborg | Stratified turbulence forced in rotational and divergent modes[END_REF].

Stratified-rotating turbulence

Stratified-rotating turbulence is influenced by rotation besides stratification. It therefore takes place beyond the mesoscale and sub-mesoscale ranges in the atmosphere and the oceans respectively. Thus, quasi-geostrophic scales are the seat of a stratifiedrotating turbulence. As already mentioned, [START_REF] Charney | Geostrophic turbulence[END_REF] discussed the similarities between quasi-geostrophic turbulence and bidimensional turbulence, hence suggesting the existence of an inverse cascade of energy in quasi-geostrophic turbulence as in two-dimensional turbulence.

More generally, stratified-rotating turbulence is characterized by an inverse cascade of energy [START_REF] Métais | Inverse cascade in stably stratified rotating turbulence[END_REF], from smaller towards larger scales. This is due to the existence of a competition between stratification and rotation, in the sense that they have antagonistic effects. Although they both tend to inhibit vertical motions, the vertical scale of the flow depends on their respective strengths: the stratification tends to diminish vertical scales whereas the rotation tends to stretch them. For example, applying a background rotation to three-dimensional turbulence generates vertical coherence [START_REF] Bartello | Coherent structures in rotating three-dimensional turbulence[END_REF]. More precisely, the dynamics of the energy transfers accross the different scales is complex in that downscale and upscale energy cascades coexist [START_REF] Pouquet | Geophysical turbulence and the duality of the energy flow across scales[END_REF]. Numerical studies of stably stratified and rotating turbulence have led to the conclusion that the amount of energy transferred downscale is dependent upon the magnitude of the rotation: a larger stratification enhances downscale transfers while a larger rate of rotation inhibits downscale transfers of energy [START_REF] Vallgren | Possible explanation of the atmospheric kinetic and potential energy spectra[END_REF][START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. [START_REF] Deusebio | Dimensional transition in rotating turbulence[END_REF] supported these results by showing that the rotation strengthens upscale transfers at the expense of the direct energy cascade, and that this is achieved through the suppression of the enstrophy production, i.e through the reduction of the dissipation. Interestingly, their results indicated that even in the presence of a rapid rotation, the inverse cascade of energy may eventually vanish when the vertical scale of the flow is large enough relative to the forcing scale at which the energy is injected, hence underlining the complexity of systems influenced by a strong stratification as well as a rapid rotation.

Routes to dissipation A global paradox

To briefly sum things up, the previous discussion has shown that:

1. Stratified rapidly rotating turbulence is the seat of an inverse cascade of energy, i.e. from smaller to larger scales. In the atmosphere and the oceans, this behaviour typically occurs at the largest scales.

2. Purely stratified turbulence is in contrast characterized by a forward cascade of energy as in three-dimensional turbulence: the energy is transferred from the largest to the smallest scales. This cascade takes place in an intermediate range of length scales in geophysical media, namely the mesoscale range in the atmosphere and the sub-mesoscale range in the oceans.

Thus, the dynamics is complex in that the energy is not directly transferred downscale throughout the whole range of geophysical scales as is the case in three-dimensional isotropic turbulence. Identifying the mechanisms responsible for the dynamical transition from a 2D-like turbulence at large scales and a more isotropic 3D-like turbulence at small enough scales actually remains a central question in the study of GFD [START_REF] Molemaker | Baroclinic instability and loss of balance[END_REF][START_REF] Müller | Routes to dissipation in the ocean: The 2D/3D turbulence conundrum[END_REF][START_REF] Ferrari | Ocean circulation kinetic energy: Reservoirs, sources, and sinks[END_REF]. Since quasi-geostrophic turbulence behaves like bidimensional turbulence as indicated by [START_REF] Charney | Geostrophic turbulence[END_REF], a transfer of energy towards small scales is possible only through so-called ageostrophic, i.e. unbalanced motions. Precisely, recent numerical simulations of stratified rotating turbulence [START_REF] Bartello | Quasigeostrophic and stratified turbulence in the atmosphere[END_REF][START_REF] Vallgren | Possible explanation of the atmospheric kinetic and potential energy spectra[END_REF][START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF][START_REF] Pouquet | Geophysical turbulence and the duality of the energy flow across scales[END_REF] have proved that an ageostrophic cascade of energy towards small scales exists. Although this mechanism seems promising to account for the dynamical transition between the large and small scale turbulent regimes, the nature of the non-linearities that are involved in this process remains unclear.

The zigzag and Kelvin-Helmholtz instabilities

Stratified fluids are the seat of three-dimensional instabilities, which are good candidates for the generation of small scales since they typically intervene in the transition towards turbulence.

A central instability that develops in the presence of an ambient stratification is the so-called zigzag instability. To briefly explain its principle, this instability acts on a pair of vertical vortices. It bends the pair as a whole in a zigzag as illustrated in figure 1.4. It is antisymmetric with respect to the plane that separates the two vortices. An important feature is that this instability does not saturate and results in the slicing of the pair in independent pancake dipoles [START_REF] Beckers | Dynamics of pancake-like vortices in a stratified fluid: experiments, model and numerical simulations[END_REF] in the ultimate stages of its evolution.

It has been discovered experimentally by Billant and Chomaz (2000a) and then studied theoretically and numerically (Billant and Chomaz, 2000b,c) in the case of counter-rotating vortices. These preliminary investigations have then been extended to the case of co-rotating vortices [START_REF] Otheguy | Elliptic and zigzag instabilities on co-rotating vertical vortices in a stratified fluid[END_REF][START_REF] Otheguy | Theoretical analysis of the zigzag instability of a vertical co-rotating vortex pair in a strongly stratified fluid[END_REF]. These analyses were finally reconciled and extended to the stratified-rotating case by [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 1. General stability equations[END_REF]; [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 2. Analytical and numerical analyses[END_REF].

Interestingly, [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF] have shown, by means of a numerical analysis, that the zigzag instability transfers the energy from large to small dissipative scales, and as such represents a contributor of the direct energy cascade that characterizes strongly stratified turbulence. Along with the fact that it slices vortices into smaller pancake-like vortices, and as such may explain the layering which is typical of stratified turbulence, this suggests that the zigzag instability is actually an essential mechanism in stratified flows.

Another major characteristic of this instability is that it is at the origin of an intense growth of the vertical shear. In fact, this intensification can even be sufficient to trigger the Kelvin-Helmholtz (KH) instability [START_REF] Thorpe | A method of producing a shear flow in a stratified fluid[END_REF][START_REF] Drazin | Kelvin-helmholtz instability of finite amplitude[END_REF] as a secondary instability, as reported by [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]; [START_REF] Augier | Onset of secondary instabilities on the zigzag instability in stratified fluids[END_REF]; [START_REF] Augier | Spectral analysis of the transition to turbulence from a dipole in stratified fluid[END_REF].

The KH instability is typically generated at the interface between two flows moving with constant and opposite velocities. Depending on the conditions, the interface can be destabilized and be the seat of roll-up events referred to as KH billows. A model experiment where this instability is studied is described in [START_REF] Thorpe | A method of producing a shear flow in a stratified fluid[END_REF]. An illustration of the billows resulting from the KH instability is provided in figure 1.5. In practice, the KH instability is likely to be triggered as soon as the vertical shear of the horizontal velocity of the flow is large enough according to Miles & Howard criterion [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Note on a paper of John W. Miles[END_REF]. As such, it can occur in any vertically sheared flow a priori.

Vertical decorrelation of a structure by an environmental shear flow

A mechanism of generation of small scales has been proposed by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] for strongly stratified flows. It consists in a process of vertical decorrelation of a structure by an environmental, vertically sheared flow. Under the assumptions that the horizontal velocity of the structure is purely advected by the shear flow and that no vertical velocity arises, Lilly demonstrated that the vertical shear of the horizontal velocity grows algebraically as time evolves. Hence, on the grounds of Miles & Howard criterion (Miles, [START_REF] Thorpe | A method of producing a shear flow in a stratified fluid[END_REF]). 1961;[START_REF] Howard | Note on a paper of John W. Miles[END_REF], one would expect the shear instability to be triggered unconditionally after a certain time that notably depends on the intensity of the advection due to the shear flow. This analysis suggests that stratified turbulence should be rapidly turned into three-dimensional turbulence via the shear instability.

Several studies have revealed that the KH instability effectively exists in strongly stratified turbulence, as predicted by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]. More precisely, [START_REF] Smyth | Length scales of turbulence in stably stratified mixing layers[END_REF]; [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF]; [START_REF] Godoy-Diana | Vertical length scale selection for pancake vortices in strongly stratified viscous fluids[END_REF]; [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]; [START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF] indicate that the so-called buoyancy Reynolds number is crucial for characterizing the dynamics of stratified turbulence, and especially for discriminating whether conditions are conducive to the onset of the shear instability or not. Indeed, [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] report the existence of KH billows in their simulations when this quantity is larger than the unit. On the contrary, they observe that the shear instability leaves place to wave-like localized disturbances as the buoyancy Reynolds number is lessened below the unit. These findings are consistent with [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF], who argued that the local Richardson number scales as the inverse of the buoyancy Reynolds number. Thus, if the buoyancy Reynolds number is large enough, then the local Richarsdon number can be smaller than the critical threshold 1/4 below which the KH instability is likely to develop [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Note on a paper of John W. Miles[END_REF]. Note that [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF] also report the presence of a KH instability in their numerical simulations. Aside from these papers, numerous complementary studies have testified of the occurrence of shear instabilities in strongly stratified flows, such as [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF] and [START_REF] Augier | Onset of secondary instabilities on the zigzag instability in stratified fluids[END_REF]. In these cases, the shear instability manifested itself as a secondary instability on the zigzag instability. These advances suggest that, apart from the zigzag instability previously described, the KH instability is a central mechanism of the direct energy cascade that is characteristic of strongly stratified turbulence, as pointed out by [START_REF] Waite | The vortex instability pathway in stratified turbulence[END_REF]. Indeed, it generates small scales by making the flow transition into three-dimensional turbulence locally.

Therefore, an open question is: does the process of vertical decorrelation imagined by Lilly always result in the trigger of the shear instability and ultimately in the generation of small scales via this mechanism ?

1.4 Goals of the dissertation 1.4.1 Dynamics of a vortex embedded in an environmental shear flow

The dynamics of the process of vertical decorrelation considered by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] has been the subject of a variety of other investigations [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF][START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF][START_REF] Demaria | The effect of vertical shear on tropical cyclone intensity change[END_REF][START_REF] Smith | On the dynamics of hurricane-like vortices in vertical-shear flows[END_REF][START_REF] Vandermeirsch | Resistance of a coherent vortex to a vertical shear[END_REF][START_REF] Jones | On the ability of dry tropical-cyclone-like vortices to withstand vertical shear[END_REF] that dealt with the case of a cyclone-like vortex subjected to an environmental shear flow. The underlying motivation of these works was to understand the different facets of the resiliency of the vortex to the external shear flow.

Although the conclusions of these papers are not directly transposable to the purely stratified case [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] because a background rotation is taken into account, they have not only shown that a vertical velocity exists due to the tilting of the vortex, but that it actually represents a major component of the global dynamics [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF][START_REF] Frank | Effects of environmental flow upon tropical cyclone structure[END_REF]Jones, 2000b). Hence, the fact that [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] ignores the vertical velocity in his analysis appears debatable. This is not the only objection that can be raised regarding the interpretation of Lilly. Indeed, [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF] have studied the dynamics of an initially aligned vortex embedded in a current difference in a two-layer quasigeostrophic flow. Their investigation has led to the conclusion that two distinct regimes exist: when the intensity of the current is weak enough compared to the vortex strength, the vortex is simply advected while remaining coherent along the vertical. This regime has been called non-tearing regime by the authors. In contrast, when the advection by the current is strong relative to the vortex, this one is progressively torn apart into two smaller vortices. This regime is called tearing regime. Hence, discriminating whether the vortex will be torn or not by the current is possible only by examining the competition between the external advection and the strength of the vortex, which is measured by its ability to remain coherent vertically. We stress that the analysis by [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF] is valid in the presence of a rapid background rotation. However, this competition still exists in purely stratified flows. Therefore, neglecting the vertical coupling as done by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] may not be appropriate in general.

Global problem

In light of the previous discussions, an important issue is to be able to state whether the vertical decorrelation of a flow structure by an environmental shear flow is a relevant and efficient mechanism of generation of small scales via the shear instability in stratified fluids or not. To say things more simply, the problem is to test Lilly's hypothesis. Hence, the concern of the present dissertation is the investigation of the dynamics of an initially vertical vortex in a stratified shear flow. Indeed, although [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] cleverly predicted the possibility of encountering shear instabilities in strongly stratified flows, recent findings [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF][START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] have demonstrated that his reasoning is insufficient to predict their occurrence in so far as it neglects important physical effects -including the vertical velocity and the vertical coupling of the velocity field -on one hand and that, on the other hand, the development of the instability is not unconditional, which directly contradicts Lilly's predictions.

Outline of the thesis

In the aim of testing Lilly's conjecture for the growth of the vertical shear, we will investigate the dynamics of a Lamb-Oseen vortex in a shear flow uniform in the horizontal and sinusoidally sheared along the vertical. We will rely on the full Navier-Stokes equations under the Boussinesq approximation to perform Direct Numerical Simulations (DNS) as well as theoretical asymptotic analyses of the dynamics of the vortex.

Before explaining the concrete results of the thesis, the chapter 2 is dedicated to the description of the numerical method. Some hints regarding the numerical simulations are also briefly given.

Then, in chapter 3, we investigate the volume-integrated budgets of kinetic energy and enstrophy. The goal of this study is to understand how the control parameters affect the saturation of the enstrophy. Indeed, when vertical gradients are large, the enstrophy is equivalent to the vertical shear of the horizontal velocity, which is a crucial quantity to predict the onset of the shear instability. It will be shown that the analysis of these global budgets is actually insufficient to account for the effect of the stratification on the variations of the enstrophy. This has therefore motivated local analyses.

These local analyses are carried out in chapter 4. The first analysis will consist in an asymptotic expansion for small times starting from the initial conditions. It will be shown that the initial response of the vortex is non-hydrostatic whatever the magnitude of the stratification. The second asymptotic analysis is conducted for long vertical wavelength following a similar approach as used by [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 1. General stability equations[END_REF] to describe the unstable interaction between columnar vortices in stratified-rotating fluids. This analysis will provide the governing equations for the displacements of the vortex center and its structure evolution due to the shear flow. We shall see that the vortex is not only advected in the direction of the shear flow but also in the orthogonal direction owing to its self-induced motion. In addition, it will be demonstrated that the vertical vorticity of the vortex decreases with time owing to dynamic and viscous effects. The initial non-hydrostatic regime evidenced in the first asymptotic analysis will be related to the transient excitation of internal waves at the start-up of the motion.

These theoretical results will eventually be discussed against numerical results from the DNS in chapter 5. A very good agreement will be observed. Hence, on the basis of these results, an estimation of the minimum local Richardson number will be derived in terms of the control parameters in the limit of strong stratification. This will allow concluding regarding the possibility of encountering the shear instability in strongly stratified fluids via the process of vertical decorrelation suggested by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF].

A general conclusion as well as perspectives will be given ultimately in chapter 6.

Chapter 2

Numerical method

The numerical program, called NS3D, that allowed us to perform direct numerical simulations (DNS) on parallel computers is presented in this chapter. Another postprocessing program has been developed in order to compute a variety of quantities that derive from the fields computed by NS3D. This Data Analysis Program (DAP) works on parallel architectures as well. It is described in a second time.

Governing equations

The dynamics of any flow in a stably stratified fluid is governed by the incompressible Navier-Stokes equations under the Boussinesq approximation. In dimensional form, they read ∇ • u = 0, (2.1)

∂u ∂t = u × ω -∇ p ρ 0 + u 2 2 + be z + ν∇ 2 u, (2.2) ∂b ∂t + u • ∇b + N 2 w = ν Sc ∇ 2 b, (2.3)
with u = (u, v, w) being the velocity field, ω its curl, p the pressure, b the buoyancy, e z the vertical unit vector oriented upwards, ν the viscosity, and Sc the Schmidt number. The total density field ρ t has been decomposed as ρ t (x, t) = ρ 0 + ρ(z) + ρ (x, t), ρ 0 being a constant reference density, ρ a linear mean density profile varying with the vertical coordinate z, and ρ (x, t) a perturbation density which is related to the buoyancy b by b = -gρ/ρ 0 . The Brunt-Väisälä frequency N measuring the ambient stable stratification is

N = - g ρ 0 dρ dz .
(2.4)

Direct Numerical Simulations

Global features

NS3D integrates the equations (2.1-2.3) in a parallelepipedic domain of dimensions l x × l y × l z with n x × n y × n z collocation points. Triply periodic boundary conditions are implemented, i.e.:

q(x + l x , y + l y , z + l z , t) = q(x, y, z, t), (2.5)

21
where the vector q = (u, v, w, p, b) gathers the unknowns of the problem. Since periodic boundary conditions are used, pseudo-spectral methods are particularly suited for resolving the problem numerically. Their principle is addressed in the following.

Formulation of the equations in spectral space

Let ϕ be a scalar periodic quantity, hence satisfying (2.5). Its three-dimensional Fourier transform F is given by: (2.6) where k = (k x , k y , k z ) is the total wavenumber and i 2 = -1. Thanks to the periodicity, a fundamental property of the Fourier transform is

F(ϕ)(k x , k y , k z , t) = lz 0 ly 0 lx 0 ϕ(x, y, z, t) e -i(kxx+kyy+kzz) dxdydz,
F ∂ϕ ∂x j = ik j F(ϕ), (2.7)
with x j being any space variable and k j the corresponding component in k. This property is essential for formulating the Boussinesq equations in spectral space. In the following, we derive the spectral form of the momentum equation (2.2) so as to illustrate some subtle aspects of pseudo-spectral methods. The spectral form of the thermodynamic equation will be given directly since it is based on a similar reasoning. First, applying the 3D Fourier transform to (2.1) yields:

k • F (u) = 0. (2.8)
Thus, in spectral space, the velocity field is orthogonal to the wave vector. Then, by taking the divergence of (2.2), the dynamic pressure per unit mass p * = p/ρ 0 + u 2 /2 is found to satisfy the following Poisson equation:

∇ 2 p * = ∇ • u × ω + be z + ν∇ 2 u .
(2.9)

As a result, we have

F(p * ) = -i k k 2 • F (a), (2.10)
where the acceleration a is defined as a = u × ω + be z . Finally, the Fourier transform is applied to the momentum equation (2.2). By introducing the unit vector e k such that

e k = k ||k|| , (2.11)
where ||k|| is the norm of k, and by taking advantage of the previous results, we find:

∂ ∂t + νk 2 F (u) = F (a) -e k • F (a) e k .
(2.12)

The right-hand side of (2.12) can be rewritten

Fm (a)e m -Fn (a) k n ||k|| k m ||k|| e m = δ mn - k m k n k 2 Fn (a)e n = P ⊥ ( F (a)), (2.13)
where the tensor P ⊥ (k) represents the projection on the space of solenoidal fields D(e k ) ⊥ :

P ⊥ (k) = δ mn - k m k n k 2 1 m,n 3 , (2.14)
where δ mn = 0 if m = n and δ mn = 1 if m = n. These results are not valid when the total wavenumber is zero. However, the incompressibility condition is automatically satisfied in this case, thus resolving the problem.

It eventually turns out that we can formulate (2.2) and (2.3) in spectral space in the form ∂ F (u) e νk 2 t ∂t = P ⊥ F (a) e νk 2 t , (2.15)

∂ F(b) e νk 2 t/Sc ∂t = -ik • F (bu) + N 2 F(w) e νk 2 t/Sc . (2.16)
The central features of the method are summarized below:

(i) In spectral space, the incompressibility imposes that the wave vector is orthogonal to the velocity field.

(ii) A projection operator, which appears in the form of a second-order tensor in the spectral formulation of the momentum equation, enforces the incompressibility.

(iii) The introduction of this projector is at the origin of the disappearance of the pressure field in the spectral formulation (2.15) of the momentum equation. Therefore, computing the pressure is actually not necessary for determining the remaining unknowns. However, it can be directly deduced from u and b if needed.

Time integration and pseudo-spectral method

Equations (2.15) and (2.16) are integrated in time using a fourth-order Runge-Kutta scheme. The core of the program is the algorithm implemented for computing the nonlinear terms found in these two equations. The procedure is the following:

1. The vorticity is calculated in spectral space, i.e. F (ω) = ik × F (u).

2. Backward Fourier transforms are applied to F (u), F (ω), and F(b). In this way, the fields u, ω, and b are found in physical space.

3. u × ω and bu are computed in the physical space.

4. F (u × ω) and F (bu) are deduced from the previous quantities by applying a forward Fourier transform. This yields the nonlinear terms found in (2.15) and (2.16).

The time advancement is performed.

This algorithm makes an extensive use of Discrete Fourier Transforms. They are calculated with a Fast-Fourier Transform (FFT) algorithm. The parallelization of NS3D and the assessment of its performances have been performed by [START_REF] Deloncle | Three-dimensional instabilities in stratified fluids[END_REF][START_REF] Deloncle | Three-dimensional instabilities in stratified fluids[END_REF], which offers a comprehensive description of the optimization of the whole program. The most limiting tasks in terms of computations are precisely the FTTs, which consume from 70 to 95% of the total computational time.

Aliasing

The Discrete Fourier Transform of a periodic field is responsible for the existence of socalled aliasing errors. They arise when nonlinear quantities are computed. To be more specific, the aliasing affects the accuracy of high-order modes by polluting their top one-third, as pointed out by [START_REF] Orszag | On the elimination of aliasing in finite-difference schemes by filtering highwavenumber components[END_REF]. [START_REF] Delbende | Réponse impulsionnelle dans les jets tournants et les sillages plans : analyse spatio-temporelle par simulation numérique directe[END_REF] illustrates how aliasing errors are produced in a single dimension space. Therefore, NS3D offers the possibility to truncate high order modes. Formally, we introduce the truncation radii R xs , R ys , and R zs and the maximum wavenumbers max (k x ), max (k y ), and max (k z ) along the three spectral directions. We apply an elliptic aliasing removal. Mathematically, this means that spectral quantities are set to zero for the wavenumbers (k x , k y , k z ) that satisfy the condition

k x R xs max (k x ) 2 + k y R ys max (k y ) 2 + k z R zs max (k z ) 2 > 1.
(2.17)

All the simulations presented in this dissertation have been performed using a classical two-third rule: R xs = R ys = R zs = 2/3.

Specificities of the present DNS

Initial state

The version of NS3D described in [START_REF] Deloncle | Three-dimensional instabilities in stratified fluids[END_REF] has been adapted in order to simulate the dynamics of an isolated columnar vortex subjected to a vertically sheared flow. More precisely, the initial state comprises a columnar Lamb-Oseen vortex embedded in a sinusoidal shear flow whose velocity field will be denoted U (z)e x with U (z) = U S sin(k z z), e x being the streamwise unit vector. The amplitude U S simply sets the intensity of the advection. The wavenumber k z of the profile U (z) is k z = 2π/l z so that a single wavelength is simulated. We have opted for a sinusoidal shear flow because it is straightforward to handle since periodic boundary conditions are implemented in NS3D. The vortex is purely vertical at t = 0. Initially, its vorticity field, in dimensional form, reads (2.18) where Γ is the circulation, a 0 the vortex radius, and r the radial coordinate. Using a Fast Fourier Transform, this vorticity field is inverted to obtain the corresponding velocity field in the computational domain. If the domain were infinite and non-periodic, this velocity field would be the Lamb-Oseen profile (2.19) e θ being the azimuthal unit vector. However, the actual profile differs from this Lamb-Oseen profile for large radius since the total circulation has to vanish in a periodic domain [START_REF] Pradeep | Effects of boundary condition in numerical simulations of vortex dynamics[END_REF]. Therefore, the total velocity field is purely horizontal initially. For this reason, the vertical velocity w and the buoyancy b are zero at t = 0. [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF].

ω = ζe z = Γ πa 2 0 e -r 2 /a 2 0 e z ,
u v = u θ (r)e θ = Γ 2πr 1 -e -r 2 /a 2 0 e θ ,
The Reynolds number Re and the Froude number F h are based on the initial conditions: .20) This implies that the time and length units are chosen as 2πa 2 0 /|Γ| and a 0 respectively. Accordingly, U S is non-dimensionalized by |Γ|/2πa 0 and k z by 1/a 0 :

Re = |Γ| 2πν , F h = |Γ| 2πa 2 0 N . ( 2
U * S = 2πa 0 |Γ| U S , k * z = a 0 k z .
(2.21)

In the following, the stars will be dropped for simplicity. We will no longer consider dimensional variables unless explicitly mentioned. All dimensionless quantities will be denoted like their dimensional counterparts for simplicity as well, as done for U S and k z . The Brunt-Väisälä frequency varies with altitude in the atmosphere and with depth in the oceans, as illustrated in figure 2.1, which displays a typical vertical profile of N in the ocean. To simplify things, it will be assumed that N is independent of the vertical coordinate z throughout this dissertation. Thus, the Froude number will be varied by simply setting N to the proper constant value.

It is worth stressing that we focus on flows such that the vortex is stronger than the external shear flow, i.e. |U S | < 1 in dimensionless form. We will also consider that the vortex circulation is positive, since there is no loss of generality in making this assumption.

Therefore, the control parameters of the problem are F h , Re, and the dimensionless velocity U S and wavenumber k z .

On the triply periodic boundary conditions

Periodic boundary conditions are commonly implemented in numerical programs that aim at investigating the dynamics of vortices beyond the range of parameters reachable in laboratory experiments. Nevertheless, such boundary conditions can disrupt the physics, leading to at least partially incorrect numerical results as pointed out by [START_REF] Pradeep | Effects of boundary condition in numerical simulations of vortex dynamics[END_REF].

For example, a pair of interacting vortices in a parallelepipedic box is actually embedded in a strain field that arises from the existence of an infinite number of pairs of vortices located in the virtual boxes that are adjacent to the domain which is simulated. This effect is precisely due to the periodicity constraint. Depending on the conditions, this strain may exert a strong influence on the pair of vortices under consideration, and, in turn, make its dynamics depart from the behaviour that should be observed normally, i.e. without any strain field. [START_REF] Pradeep | Effects of boundary condition in numerical simulations of vortex dynamics[END_REF] also stress that the net circulation of the flow within the numerical domain should be zero because of the periodic boundary conditions. The flow considered throughout this dissertation comprises a single vortex, which is at the origin of the existence of a non-zero circulation at the boundaries. In order to suppress this effect, [START_REF] Otheguy | Elliptic and zigzag instabilities on co-rotating vertical vortices in a stratified fluid[END_REF] have proposed to add the appropriate background rotation that compensates this net circulation. This background rotation can be estimated from the mean vertical vorticity uniformly generated by the vortex in the domain, which is equal to Γ/l x l y .

Rennich and Lele (1997) also provided a method to deal with problems related to periodic boundary conditions.

However, in order to mitigate the effect of periodicity, the simplest solution consists in attributing sufficiently large values to the domain dimensions, as in [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]. Obviously, this requires carrying out a careful preliminary analysis for identifying the proper values of the different parameters.

We have opted for the latest option for a variety of reasons, and notably because:

1. Implementing a method such as proposed by [START_REF] Rennich | Numerical method for incompressible vortical flows with two unbounded directions[END_REF] would have required to profoundly modify NS3D.

2. The strain due to the artificial images of the flow is inversely proportional to the square of the horizontal dimension of the numerical domain. As such, it is weak when l x and l y are large. This is valid also regarding the net circulation.

In the case of the model flow that we focus on in this dissertation, we have observed that the dimensions of the fluid domain should be adjusted depending on the intensity of the advection U S . This point will be detailed later since, to be explained satisfactorily, it requires introducing several intermediate notions.

Post-processing

NS3D computes the velocity field u = (u, v, w), its curl ω = (ω x , ω y , ζ), and the buoyancy b. A Data Analysis Program (DAP) has been therefore developed from NS3D in order to carry out further calculations from these fields. Its architecture is very similar to that of NS3D.

Since a wide variety of post-processing routines have been implemented, we will not describe the DAP in details. Instead, we will put the emphasis on the computation and processing of a crucial quantity for analysing the vortex deformations, namely the potential vorticity. The other tasks performed by the DAP will be listed and their subtleties will be briefly pointed out if needed.

Potential vorticity

The potential vorticity Π is an inviscid invariant of the governing equations (2.1-2.3), i.e.:

DΠ Dt = 0, D Dt = ∂ ∂t + u • ∇. (2.22)
In a stratified fluid, its expression, in dimensional form, is:

Π = ω • ∇b + N 2 e z .
(2.23)

The potential vorticity is one of the most challenging quantities to compute from the basic fields because it is the sum of three nonlinear terms. As such, it is a priori sensitive to aliasing errors. For this reason, we have computed this quantity without introducing any truncation of the high-order modes in the first instance. In a second time, it has been calculated by applying the same truncation as in NS3D, namely an elliptic dealiasing following the two-third rule by [START_REF] Orszag | On the elimination of aliasing in finite-difference schemes by filtering highwavenumber components[END_REF].

The comparison of the numerical results has revealed that differences are insignificant provided that the simulation is properly resolved. In contrast, they can grow quite significantly when the mesh is coarse, and especially when a shear instability is under-resolved. In order to illustrate that the computation of the potential vorticity is largely insensitive to the aliasing removal in properly resolved cases, figure 2.2 reports the evolution of the potential vorticity maxima computed in the horizontal planes z = l z /4, where the advection is the strongest and the vertical shear of the horizontal velocity is the weakest, and z = l z /2, where the advection is zero and the vertical shear of the horizontal velocity reaches its peak. The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000. The dimensions of the box are l x = l y = 18 and l z = 4. This run has been performed using n x = n y = 512 and n z = 256 collocation points. The fact that the curves are collapsed confirms the previous statement. For this reason, the DAP actually computes the potential vorticity without truncating the top one-third modes.

Tracking the vortex position

It will be useful to locate the vortex as time evolves. This is done by computing the position of the vortex center in each horizontal plane, which provides a representation of the deformations of the vortex axis.

A panel of different methods can be implemented, depending on how we define the center. In a similar study, [START_REF] Jones | On the ability of dry tropical-cyclone-like vortices to withstand vertical shear[END_REF] stressed that the outer regions of the vortex The elliptic dealiasing cancels the top one-third modes following [START_REF] Orszag | On the elimination of aliasing in finite-difference schemes by filtering highwavenumber components[END_REF]. The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000.

were more strongly tilted than its inner core, in agreement with former investigations such as [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF] and Jones (2000a). She emphasizes that, in such cases, defining the position of the vortex center as that of the maximum vertical vorticity does not provide reliable indications about the bending of the vortex because of the non-uniformity of the vortex tilt. This is the reason why we have opted for a similar definition as in [START_REF] Jones | On the ability of dry tropical-cyclone-like vortices to withstand vertical shear[END_REF], which is precisely based on the use of the potential vorticity. To be more specific, the position of the vortex (x c (z, t), y c (z, t)) in each horizontal plane has been obtained from the potential vorticity centroids:

x c (z, t) = xΠ h Π h , y c (z, t) = yΠ h Π h , (2.24)
where the brackets denote the horizontal averaging operator:

ϕ h = Π>Πc ϕ(x, y, z, t)dxdy.
(2.25)

The horizontal integration is carried out only in the regions where the potential vorticity is larger than a threshold Π c in order to avoid taking into account the small background vorticity due to the fact that the total vorticity is zero owing to the use of periodic boundary conditions. This threshold being defined as Π c = η max t=0 (Π), several tests have been done in order to set the coefficient η properly, that is for ensuring that the numerical results are independent of the size of the computational domain and of the particular value of Π c . This investigation has revealed that taking η = 0.05 was sufficient. Indeed, figure 2.3 displays four trajectories of the vortex center computed as previously explained, that correspond to four values of the parameter η, namely η = 0.05, η = 0.1, η = 0.5, and η = 0.9. The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000. The horizontal dimensions of the box are l x = l y = 18. The displacements ∆x and ∆y that are plotted are the displacements of the vortex with respect to its initial position, in the plane z = l z /4 where the advection is the strongest. They are reported for 0 t 35. It is seen that the trajectories are converged for η 0.1, hence confirming that η = 0.05 is an optimal choice. It is worth mentioning that the same method has been also implemented to locate the vortex on the basis of the vertical vorticity field. This simply consists in replacing Π by ζ in (2.24). This is different from the method that consists in assimilating the vortex center to the point where the vertical vorticity is maximum since the horizontal average will take the potential asymmetries into account and will therefore offer an adequate representation of the bending of the vortex. This is not the only advantage of such a method since it will also allow estimating the consistency of such a representation of the vortex deformations. Indeed, in chapter 5, we will show that apart from an initial transient stage, this computation of the vortex center provides nearly identical results whether Π or ζ is used.

Mitigation of the effect of periodic boundary conditions

As already mentioned, we have simulated the dynamics of the flow in a box whose horizontal dimensions are large enough to significantly weaken the effect of the periodic boundary conditions.

A careful analysis has revealed that a good indicator of how influential periodic boundary conditions are is the convergence of the vortex trajectory when l x and l y are varied, following the same idea as in figure 2.3.

More precisely, to any value U S of the advection corresponds a minimal threshold for l x and l y beyond which the vortex trajectories -as previously computed -are independent of the box dimensions, i.e. insensitive to the periodic boundary conditions. Without presenting the whole set of investigations that have been carried out, figure 2.4 shows that for U S = 0.1, the horizontal dimensions of the box must be l x,y 24 to ensure that the vortex trajectory is reasonably converged. Such a large domain requires introducing a large number of collocation points to properly resolve the finest scales. These requirements are actually too constraining in terms of computational resources. Instead, for U S = 0.2, setting l x,y 18 is sufficient and allows using an affordable number of collocation points to properly resolve the dynamics. The latest compromise being the best to reconcile numerical issues and the condition |U S | < 1, the vast majority of the simulations presented in this dissertation have been performed with l x = l y = 18 and U S comprised between 0.2 and 0.4.

Listing of the post-processing routines

To finish, we provide a list of the main routines that have been implemented in the DAP:

1. The potential vorticity field is computed without applying any dealiasing. Its maximum is also computed and located in each plane z = cte. Similarly, the position of its centroid is calculated in each of these planes.

2. The global budgets of kinetic and potential energies, enstrophy and potential enstrophy are computed from the basic fields determined by NS3D.

3. The total flow is decomposed into a mean component corresponding to a purely horizontal, vertically sheared flow, and a residual component which is in fact assimilable to the vortex flow. This is demonstrated in details in the next chapter.

4. The global and local budgets of kinetic energy and enstrophy associated with these two components are computed.

5. The vertical shear of horizontal velocity is computed for the total flow, the mean flow, and the vortex flow. The corresponding maxima are also determined in each plane z = cte. Similarly, the vertical shear of the total horizontal velocity on one hand, and of the vortex horizontal velocity on the other hand, are determined at the center of the vortex, in the plane z = l z /2 where the vertical shear is the most intense. Indeed, it will be shown that this quantity is essential to predict whether the shear instability can be triggered or not.

6. The vertical buoyancy gradient as well as the local Richardson number (see its definition further) are computed. A quantity of utmost importance for examining the onset of the shear instability is the minimum Richardson number, which is therefore determined in each horizontal plane. We do not apply any dealiasing to compute the Richardson number.

7. Globally, the DAP allows determining the vertical profiles of the minimum and maximum values of any field. In particular, this is done for the vertical velocity and the buoyancy, in addition to the fields already mentioned.

8. The profile of the mean flow velocity is also computed by the program.

Towards practical numerical simulations

We finish this chapter by giving some preliminary hints regarding the DNS.

As already mentioned, the present investigation aims at discriminating whether shear instabilities are likely to be triggered or not, depending on the values attributed to the control parameters. Following [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF] and [START_REF] Howard | Note on a paper of John W. Miles[END_REF], a necessary condition for the KH instability to develop in a stratified fluid is that the local Richardson number falls below 1/4 somewhere in the flow. Recent studies, such as [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF], [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] and [START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF], have shown that the buoyancy Reynolds number ReF 2 h (where Re is the Reynolds number and F h the Froude number) is of utmost importance in determining the dynamics of stratified turbulence. [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF] and [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF], besides reporting the existence of shear instabilities in strongly stratified flows, both proposed that the local Richardson number is inversely proportional to ReF 2 h . Thus, a large buoyancy Reynolds number is a priori conducive to the occurrence of the KH instability.

This is important for setting the horizontal Froude number as well as the Reynolds number in the numerical simulations. On one hand, the fact that the fluid is strongly stratified imposes F h < 1, and, if possible, F h 1. On the other hand however, since ReF 2

h must be at least reasonably large to favor the onset of shear instabilities, say

ReF 2 h 1, decreasing F h requires increasing Re if ReF 2
h is kept constant. Therefore, tending to the strongly stratified regime at a constant buoyancy Reynolds number requires increasing the resolution so as to resolve the finest scales, hence requiring larger computational resources. Thus, the setting of F h and Re is a matter of compromises. A variety of tests have shown that these issues are reconciled by choosing F h = 0.5 and increasing Re up to ∼ 7000 for example, whereas with F h = 0.1 the Reynolds number can be set up to Re ∼ 10000. With these settings, good but affordable resolutions are sufficient for properly resolving all the scales.

Note that, due to the Schmidt number Sc, the competition between inertial and viscous effects is not the same in the thermodynamic equation (2.3) as in the momentum equation (2.2). This is at the origin of another numerical difficulty. Indeed, for salt in the oceans for example, we have Sc ∼ 700. The "equivalent" Reynolds number that intervenes in (2.3) is therefore 700 times larger than in (2.2). Thus, in the DNS, if we opt for Re = 6000 by keeping Sc = 700 at the same time, then the actual Reynolds number in (2.3) is ∼ 4.2 × 10 6 . Even with Sc = 7, we would have an equivalent Reynolds number ∼ 4.2 × 10 4 . Both cases are out of reach numerically if simulated through DNS. For numerical convenience, we have therefore opted for keeping Sc = 1.

Chapter 3 Preliminary study: global energy and enstrophy budgets

In this preliminary investigation, we examine the global energy and enstrophy budgets in order to get some insights on the effects of the control parameters on the dynamics of the flow. Especially, the present analysis aims at understanding how the enstrophy is affected when the control parameters vary. Indeed, as this quantity is nearly equivalent to the vertical shear of the horizontal velocity when vertical gradients are large, this study is a first step to predict the development of the shear instability.

Example of a simulation

Before discussing the energy and enstrophy budgets, we show in figure 3.1 the evolution of the potential vorticity field and the vertical shear of the horizontal velocity for the parameters F h = 0.5, k z = π, U S = 0.2, and Re = 6000. The purpose is here to just give a brief idea of the flow evolution. This evolution will be discussed in more details in chapter 5.

Decomposition of the flow

In order to understand the mechanisms at work in the flow, it is interesting to decompose the total flow as

u = ū + u * , (3.1)
where the operator • denotes the horizontal average which for any quantity q is defined as q = 1 l x l y ly 0 lx 0 q(x, y, z, t)dxdy.

(3.2) Thus, ū is a mean flow varying only along the vertical and with time. The horizontal average of the complementary flow u * is zero by definition. At t = 0, we have ū = U (z)e x and u * = u v so that u * will be called "vortex flow". Applying the horizontal average operator to the dimensionless version of (2.1-2.

∂ ū ∂t + ∂w (ū + u * ) ∂z = - ∂ p ∂z e z + be z + 1 Re ∂ 2 ū ∂z 2 , (3.4) ∂ b ∂t + ∂wb ∂z = - w F 2 h + 1 ReSc ∂ 2 b ∂z 2 , (3.5)
where the periodicity of (u, p, b) has been used. Equation (3.3) implies w = 0. Indeed, this equation rigorously indicates that w depends on t only. However, an integration along the vertical then yields this result. Thus, the mean flow is purely horizontal and governed by

∂ ūh ∂t = - ∂w * u * h ∂z + 1 Re ∂ 2 ūh ∂z 2 . (3.6)
The equations for u * can be obtained by subtracting (3.6) from the dimensionless form of (2.2).

Interestingly, the analysis of the numerical results indicates that the mean flow exhibits a velocity field which closely resembles the sinusoidal profile imposed initially. This is illustrated in figure 3.2, where the vertical profiles of each component ū and v of the mean flow are reported at four different times t = 0, t = 10, t = 20 and t = 30, for F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000. It is seen that the departures from the initial sinusoidal profile are small. Hence, ūh = U (z)e x is an accurate approximation of the mean flow velocity at any time.

Kinetic and potential energies

Definitions

First, we define the volume-integrated kinetic energies of the mean flow on one hand, and of the vortex on the other hand, following the decomposition of the total flow introduced previously. The global kinetic energy of the mean flow simply reads:

Ēk = 1 l z ϑ ūh 2 2 dϑ, (3.7)
where ϑ stands for the numerical domain. The integral over ϑ is the triple integral over the fluid domain. The normalisation by l z allows performing comparisons between cases associated with distinct wavelengths [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]. We also define:

E * kh = 1 l z ϑ u * h 2 2 dϑ, E * kz = 1 l z ϑ w * 2 2 dϑ, E * k = E * kh + E * kz , (3.8)
where E * kh is the global horizontal kinetic energy of the vortex and E * kz its vertical counterpart. Hence, E * k is the total kinetic energy of the vortex per vertical length unit.

In addition, the potential energy per vertical length unit is:

E p = 1 l z ϑ F 2 h b 2 2 dϑ.
(3.9)

Volume-integrated budgets

The volume-integrated budgets that describe the variations of the kinetic energies are found to be d Ēk dt = -T + D (3.10)

for the mean flow, and

dE * k dt = C * + T + D * h + D * z + D w (3.11)
for the vortex. In these budgets, the transfer terms read

T = 1 l z ϑ ūh • ∂w * u * h ∂z dϑ, C * = 1 l z ϑ w * bdϑ, (3.12)
and the dissipation terms are

D = - 1 l z Re ϑ ∂ ūh ∂z 2 dϑ, (3.13) D * h = - 1 l z Re ϑ   ∂u * h ∂x 2 + ∂u * h ∂y 2   dϑ, D * z = - 1 l z Re ϑ ∂u * h ∂z 2 dϑ, (3.14) D w = - 1 l z Re ϑ   ∂w * ∂x 2 + ∂w * ∂y 2 + ∂w * ∂z 2   dϑ. (3.15)
The kinetic energy Ēk of the mean flow varies because of two terms: a transfer term T from the mean flow to the vortex and the dissipative term D which results from the vertical shear of the mean flow. The budget of the kinetic energy E * k of the vortex is more complex and contains:

1. the transfer term T and the conversion of the potential energy into kinetic energy C * .

2. two dissipative terms related to the horizontal shear and the vertical shear of the vortex horizontal velocity (3.14). We stress that these two contributions have been distinguished in the total dissipation for convenience in anticipation of further analyses.

3. the dissipative term (3.15) related to the vertical velocity.

Finally, the total potential energy satisfies:

dE p dt = -C * - F 2 h l z ReSc ϑ (∇b) 2 dϑ. (3.16)
Since this equation only involves the vertical velocity w = w * , it remains the same whether the total horizontal velocity field is decomposed into a mean component and a residual part or not.

Enstrophy

Definitions

From the velocity field decomposition defined previously, the vorticity ω = ∇ × u of the global flow can be expressed as:

ω = ω + ω * , ω = ∇ × ūh , ω * = ∇ × u * . (3.17)
As such, ω and ω * respectively stand for the vorticities associated with the mean flow and the vortex. Hence, we define:

Z = 1 l z ϑ ω2 2 dϑ, Z * = 1 l z ϑ ω * 2 2 dϑ, (3.18)
Z being the volume-integrated enstrophy associated with the mean flow and Z * the global enstrophy of the vortex. Note that ω is purely horizontal since ū depends on z only. The application of the horizontal average to the non-dimensional vorticity equation gives:

∂ ω ∂t + ∂w * ω * ∂z = ∂ω * z u * ∂z + 1 Re ∂ 2 ω ∂z 2 . (3.19)
The equation for ω * can be easily obtained by subtracting (3.19) from the full vorticity equation.

Volume-integrated budgets

The global enstrophy budgets can be written in the form

d Z dt = -T + S + D (3.20)
in the case of the mean flow, and

dZ * dt = T + S * + S * + B * + D * (3.21)
regarding the vortex. The different terms read:

T = 1 l z ϑ ω • ∂w * ω * ∂z dϑ, (3.22) S = 1 l z ϑ ω • ∂ω * z u * ∂z dϑ, (3.23) S * = 1 l z ϑ ω * • [(ω * • ∇) u * ] dϑ, S * = 1 l z ϑ ω * • [( ω • ∇)u * + (ω * • ∇) ūh ] dϑ, (3.24)
B * = 1 l z ϑ ω * • (∇b × e z ) dϑ, (3.25) D = - 1 l z Re ϑ ∂ 2 ūh ∂z 2 2 dϑ, D * = - 1 l z Re ϑ (∇ω * x ) 2 + (∇ω * y ) 2 + (∇ω * z ) 2 dϑ.
(3.26) The term T represents a transfer from the enstrophy of the mean flow to the enstrophy of the vortex. The term S represents the effect of the stretching and tilting of the vortex on the mean flow. Conversely, S * is the stretching and tilting of the vortex by the mean flow. S * is the stretching and tilting of the vortex on itself. B * is the effect of the baroclinic torque and, finally, D and D * are the dissipations of the mean flow and vortex enstrophies, respectively.

Energy and enstrophy evolutions

Figure 3.3 presents typical evolutions of the kinetic and potential energies in (a), and of the enstrophy in (b), for F h = 0.5, k z = π, U S = 0.2 and Re = 6000. The vertical kinetic energy is always small relative to its horizontal counterpart. This remains true throughout this dissertation. This is consistent with the fact that strong stratification inhibits vertical motions. The kinetic energy of the mean flow is nearly constant although it increases slightly owing to the transfer term T (see (3.12)). This is consistent with the fact that the velocity field of the mean flow remains close to the sinusoidal profile imposed initially (figure 3.2). In spite of the KH instability present at t 24 in this simulation, the global horizontal kinetic energy E * kh of the vortex decreases regularly and smoothly, following an approximately linear trend. Since the mean flow has an almost constant kinetic energy, the resulting variations of the total horizontal kinetic energy E kh are very similar to those of the vortex horizontal kinetic energy. The potential energy remains small compared to the horizontal kinetic energies. It exhibits a first regime where it grows regularly. Then, it saturates approximately when the KH instability develops. In the present case, the instability decays quite rapidly and disappears at t 30. Therefore, the plateau seen in the evolution of the potential energy persists after the decay of the instability. Ultimately, the potential energy slowly decays.

Regarding the enstrophy, that of the mean flow is almost constant as expected. In contrast, the global horizontal enstrophy of the vortex defined as

Z * h = 1 l z ϑ ω * h 2 2 dτ (3.27)
strongly increases. Since the vertical velocity is small compared to the horizontal velocity, the horizontal enstrophy ω * h 2 /2 corresponds in very good approximation to the vertical shear of the horizontal velocity of the vortex:

ω * h 2 ∂u * h ∂z 2 .
(3.28) Therefore, Z * h is a measure of the vertical shear of the vortex. It grows since the vortex is progressively bent by the external shear flow. The vortex horizontal enstrophy saturates at t 40, i.e. well after the time t 24 when the KH instability is the most developed. The maximum horizontal enstrophy of the vortex is approximately twice the horizontal enstrophy of the mean flow. Thus, there is an approximate 200% increase of the total horizontal enstrophy. In contrast, the vertical enstrophy defined as As seen in figure 3.4a, the kinetic energy of the vortex decreases because of three main effects. There are first a transfer to the potential energy (C * ) and also a transfer to the shear flow (T ). Later, these two terms decrease while the dissipation due to the vertical shear of the horizontal velocity (D * z ) increases and saturates at t 40. The other dissipative terms are very small (D * h + D w ). Figure 3.4b shows that the enstrophy of the vortex grows because of the stretching and tilting, first by the shear flow ( S * ), and then by the vortex on itself (S * ) also. In contrast, the baroclinic term (B * ), the transfer to the enstrophy of the vortex (T ), and dissipation (D * ) are all negative. All these terms saturate between t 20 and t 30. The dominant terms are the stretching S * and the dissipation D * while the stretching S * contributes also significantly to the budget but grows later and is at best of the order of S * . Therefore, in order to estimate the maximum enstrophy of the vortex (figure 3.3a), we can, as a first approach, approximate the enstrophy equation of the vortex in the form dZ * dt S * + D * .

Z z = Z * z = 1 l z ϑ ω *

Kinetic energy and enstrophy budgets of the vortex

(3.30)

The enstrophy of the vortex thus saturates when S * -D * . The order of magnitude of these two quantities can be estimated by introducing the scaling relations ||ω * || ∼ 1/δ where 1/δ ∼ ∂/∂z is the vertical length scale associated to the vertical shear (non-dimensionalized by a 0 ), ||∂u * /∂x|| ∼ ||∂v * /∂x|| ∼ ||∂u * /∂y|| ∼ ||∂v * /∂y|| ∼ O(1) since the velocities and lengths have been non-dimensionalized by |Γ|/2πa 0 and a 0 , and

|| ω|| ∼ k z U S . This gives S * ∼ l x l y k z U S δ , D * ∼ -l x l y 1 δ 4 Re , (3.31)
where the factor l x l y comes from the volume integration. Therefore, the relation S * -D * yields the order of magnitude

δ ∼ 1 (k z U S Re) 1/3 (3.32)
when the enstrophy of the vortex saturates. Hence, the maximum enstrophy of the vortex is max

(Z * ) ∼ l x l y δ 2 ∼ l x l y (k z U S Re) 2/3 . (3.33)
This scaling is independent of the Froude number F h as found by [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]. However, the scaling of [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF] is different and has been obtained by balancing S * and D * since there were no mean shear flow in their case. Since

S * ∼ l x l y /δ 2 , they have obtained δ ∼ 1/ √ Re and max (Z * ) ∝ Re instead of (3.32)- (3.33).
By assuming that the equation (3.30) holds approximately at any time, we can further obtain the typical time scale of the evolution of Z * . Indeed, (3.31) can be rewritten

S * ∼ l x l y k z U S √ Z * , D * ∼ - Z * 2 l x l y Re , ( 3.34)
since Z * ∼ l x l y /δ 2 . Thus, the simplified equation (3.30) for the enstrophy budget of the vortex can be written

dZ * dt l x l y k z U S √ Z * - Z * 2 l x l y Re . (3.35)
By rescaling the enstrophy and the time as (3.36) where

Z * = l x l y (k z U S Re) 2/3 Ẑ, t = τ v t,
τ v = Re 1/3 /(k z U S ) 2/3 , (3.35) becomes: d Ẑ d t Ẑ -Ẑ2 . (3.37)
This shows that the enstrophy of the vortex evolves as

Z * l x l y (k z U S Re) 2/3 = Ẑ = ϕ t τ v . (3.38)
If we neglect the transfer terms T and C * in the kinetic energy budget (3.11), i.e.

dE * k dt D * h + D * z + D w - 2Z * Re , (3.39)
we see that the kinetic energy should also follow the self-similar law

dE * k d t -2l x l y ϕ( t) (3.40) giving E * k = l x l y Φ t τ v . ( 3
k z = π/2 k z = π k z = 3π/2 0 1 2 3 1.5 2 2.5 3 3.5 4 4.5 5 t/τ v (b) k z = π/2 k z = π k z = 3π/2 Figure 3.5: Kinetic energy of the vortex E * k for F h = 0.5, U S = 0.2, Re = 6000, for k z = π/2 (dash-dotted line), k z =
π (dashed line) and k z = 3π/2 (solid line), as a function of (a) the time t and (b) the rescaled time t/τ v .

Parametric analysis

We now test the self-similar laws (3.38) and (3.41) by varying the control parameters. We first consider the kinetic energy of the vortex and then the corresponding enstrophy. Similarly, figure 3.6 displays the evolution of the vortex kinetic energy for various U S , for F h = 0.5, k z = π and Re = 6000. Again, the rescaling of t by τ v (figure 3.6b) leads to a self-similar decay in agreement with (3.41).

Kinetic energy of the vortex

The effect of the Reynolds number is tested in figure 3.7 for F h = 0.5, k z = π and U S = 0.2. It is seen that the evolution of E * k weakly depends on the Reynolds number, especially the initial part t 20. When t is rescaled by τ v (figure 3.7b), the collapse of the three curves remains satisfactory. This is because τ v depends on Re with the power law exponent 1/3. Moreover, the relative variations of Re in figure 3.7 are smaller than for the two previous parameters.

All these results are consistent with (3.41). Remarkably, this relation is not only valid for large times but also for short times. This suggests that the transfer terms T and C * , although dominant at the beginning, do not have a so important effect.

Enstrophy of the vortex

The self-similar law for the enstrophy of the vortex (3.38) is now checked when the parameters k z , U S , and Re are varied as done for the kinetic energy. Figures 3.8, 3.9, and 3.10 show that the scaling (3.38) gathers relatively well the different curves. The

k z = π/2 k z = π k z = 3π/2 0 1 2 3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t/τ v δ 2 Z * (b) k z = π/2 k z = π k z = 3π/2
Figure 3.8: Enstrophy of the vortex for F h = 0.5, U S = 0.2, Re = 6000, for k z = π/2 (dash-dotted line), k z = π (dashed line) and k z = 3π/2 (solid line). The unscaled enstrophy Z * is plotted against t in (a) and the rescaled enstrophy δ 2 Z * is plotted against the rescaled time t/τ v in (b).

maximum enstrophy of the vortex is approximately max (Z * ) 0.3/δ 2 and is attained at t 1.6τ v .

Effect of the Froude number

Finally, we examine the influence of F h on the evolutions of E * k and Z * . According to the self-similar laws (3.38) and (3.41), the Froude number should not have any effect. However, figure 3.11, which reports the evolution of E * k for different Froude numbers for k z = π, U S = 0.2, and Re = 6000, shows that the decay of E * k is more rapid at the beginning when F h increases. Nevertheless, we can notice that the three curves are approximately parallel for long times t 20. This shows that the transfer terms T and C * have, as a matter of fact, a non-negligible effect at the beginning while, after, the kinetic energy decay is mostly due to viscous effects. Figure 3.11b shows that the different kinetic energy evolutions collapse approximately if the time is rescaled by

τ v /F 1/3 h .
Similarly, figure 3.12 illustrates the effect of F h on Z * for the same parameters as previously. Again, the Froude number is found to exert a strong influence on the variations of Z * , thus invalidating (3.38). The different curves gather if Z * is rescaled by 1/F 1/3 h but not the time (figure 3.12b), suggesting the modified self-similar law:

F 1/3 h δ 2 Z * l x l y = ϕ t τ v . (3.42)
We have not succeeded in obtaining an interpretation for the factor F 1/3 h from the energy and enstrophy budgets. In particular, concerning the kinetic energy, the transfer towards potential energy weakens when the Froude number decreases. This explains qualitatively why the kinetic energy decay is less pronounced at the beginning for small Froude numbers. However, why this leads to the typical decay time τ v /F 1/3 h remains unknown from the previous analysis of the enstrophy and energy budgets. (a)

F h = 0.1 F h = 0.25 F h = 0.5 F 1/3 h t/τ v 0 1 2 3 1 2 3 4 5 (b) F h = 0.1 F h = 0.25 F h = 0.5 Figure 3.11: Kinetic energy E * k of the vortex for k z = π, U S = 0.
2, Re = 6000, for F h = 0.1 (dash-dotted line), F h = 0.25 (dashed line) and F h = 0.5 (solid line), as a function of (a) the time t and (b) the rescaled time (a)

F 1/3 h t/τ v .
F h = 0.1 F h = 0.25 F h = 0.5 t/τ v 0 1 2 3 4 F 1/3 h δ 2 Z * 0 0.05 0.1 0.15 0.2 0.25 (b) F h = 0.1 F h = 0.25 F h = 0.5 Figure 3
.12: Enstrophy of the vortex for k z = π, U S = 0.2, Re = 6000, for F h = 0.1 (dash-dotted line), F h = 0.25 (dashed line) and F h = 0.5 (solid line). The unscaled enstrophy Z * is plotted against t in (a) and the rescaled enstrophy

F 1/3 h δ 2 Z * is plotted against the rescaled time t/τ v in (b).

Conclusions

In the present chapter, we have investigated the global budgets of the volume-integrated kinetic energy and enstrophy in order to get some insight on the effect of the control parameters. The analysis of the budgets has provided an estimation of the minimum shear layer thickness δ which governs the maximum enstrophy of the vortex Z * ∝ 1/δ 2 . A scaling for the time scale that characterizes the viscous decay has also been obtained. These results have been shown to account for the influence of all the control parameters, except the horizontal Froude number. Indeed, the analysis of the budgets does not predict any dependency with respect to this parameter, although it strongly influences the variations of E * k and Z * . This paradox has therefore motivated an indepth investigation of the local dynamics in the next chapters.

Abstract

In this chapter, we investigate the dynamics of an initially vertical vortex embedded in a shear flow in a stratified fluid by means of two asymptotic analyses. The ultimate goal is to determine, whether or not, the Kelvin-Helmholtz instability can develop as speculated by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]. The analyses are mostly performed in the cases of the Lamb-Oseen vortex profile and a shear flow uniform in the horizontal and varying sinusoidally along the vertical.

The first asymptotic analysis is performed for small time: t 1, and shows that the dynamics is initially non-hydrostatic until t > 1/N , where N is the Brunt-Väisälä frequency.

The second asymptotic analysis is conducted by assuming a small vertical wavenumber k z rescaled by the horizontal Froude number

F h : k z F h 1.
The results show that the vortex axis is not only advected in the direction of the shear flow but also in the perpendicular direction owing to the self-induced motion of the vortex. In addition, internal waves are transiently excited at the beginning, explaining the non-hydrostatic regime evidenced in the first asymptotic analysis. Their effects on the displacements of the vortex axis are weak except at the beginning. The angular velocity of the vortex decays owing to dynamic and viscous effects. The former effect is due to the squeezing of isopycnals in the vortex core which implies a decrease of the vertical vorticity to satisfy potential vorticity conservation. In addition, a horizontal velocity field with an azimuthal wavenumber m = 2 is generated meaning that the shape of the vortex becomes slightly elliptical. These predictions will be tested against Direct Numerical Simulations in chapter 5.

Introduction

Stratified turbulence is strongly anisotropic with predominantly horizontal motions but high shear along the vertical [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]. Such layered structure inherently promotes the shear instability which breaks down the layers into smallerscale turbulence [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF][START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. This is one of the main mechanisms that accomplish the direct energy cascade towards small scales in stratified turbulence [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF][START_REF] Augier | Onset of secondary instabilities on the zigzag instability in stratified fluids[END_REF][START_REF] Augier | Spectral analysis of the transition to turbulence from a dipole in stratified fluid[END_REF][START_REF] Waite | The vortex instability pathway in stratified turbulence[END_REF].

In his pioneering paper on strongly stratified flows, [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] has conjectured a simple mechanism for the growth of the vertical shear. He considered a uniform horizontal flow directed along a given direction, let us say x, and vertically sheared: U = U (z)e x where z is the vertical coordinate. If a coherent structure with initial horizontal velocity u h0 (x, y, z) is embedded within this shear flow, it will be advected differentially depending on the vertical position:

∂u h ∂t = -U ∂u h ∂x . (4.1)
Lilly has further assumed that the coupling along the vertical is weak in the limit of strong stratification. When the vertical coupling is zero, the flow consists in a stack of horizontal two-dimensional flows evolving independently at each level z. Then, the solution at any time is simply: y, z). This shows that the vertical shear should grow algebraically with time according to

u h (x, y, z, t) = u h0 (x -U (z)t,
∂u h ∂z = - dU dz t ∂u h0 ∂x + ∂u h0 ∂z . (4.2)
On this basis, [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] argued that the Richardson number should reach sooner or later the critical value 1/4 for the onset of the shear instability [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Note on a paper of John W. Miles[END_REF]. However, a key assumption behind this conjecture is that there is no vertical coupling. This derives directly from Lilly's assumption that the vertical Froude number tends to zero in the limit of strong stratification. It is now acknowledged that the vertical Froude number remains of order unity because the typical vertical length scale scales like the buoyancy length scale in strongly stratified flows [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF].

It is therefore unknown whether or not this mechanism for the growth of the vertical shear is valid. In order to answer this question, we have considered the particular case of an initially vertical vortex in an ambient vertically sheared flow in a stratified fluid. This flow configuration is sketched in figure 4.1 together with the example of a direct numerical simulation illustrating the vortex evolution. In the present chapter, we shall study the dynamics of this flow by asymptotic analyses while direct numerical simulations will be performed in chapter 5.

The dynamics of an isolated vortex in an ambient vertical shear has been studied before in stratified rotating fluids in order to understand the dynamics of atmospheric tropical cyclones [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF](Jones, , 2000a(Jones, ,b, 2004;;[START_REF] Sutyrin | Intense vortex motion in a stratified fluid on the beta-plane: an analytical theory and its validation[END_REF][START_REF] Frank | Effects of environmental flow upon tropical cyclone structure[END_REF][START_REF] Reasor | A new look at the problem of tropical cyclones in vertical shear flow: vortex resiliency[END_REF][START_REF] Päschke | Motion and structure of atmospheric mesoscale baroclinic vortices: dry air and weak environmental shear[END_REF] or oceanic vortices [START_REF] Vandermeirsch | The net advective effect of a vertically sheared current on a coherent vortex[END_REF][START_REF] Vandermeirsch | Resistance of a coherent vortex to a vertical shear[END_REF]. In particular, [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF], using a three-dimensional continuous hydrostatic primitive equation model, has shown that an initially barotropic vortex in a linear shear flow tilts in the shear direction but also executes a cyclonic rotation about its mid-level center. [START_REF] Päschke | Motion and structure of atmospheric mesoscale baroclinic vortices: dry air and weak environmental shear[END_REF] have derived asymptotically the equations for the motion of the vortex center in the case of a weak shear flow.

The influence of a vertical shear on an initially aligned vortex has been also investigated by means of quasigeostrophic two-layer models [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF][START_REF] Demaria | The effect of vertical shear on tropical cyclone intensity change[END_REF][START_REF] Smith | On the dynamics of hurricane-like vortices in vertical-shear flows[END_REF][START_REF] Vandermeirsch | The net advective effect of a vertically sheared current on a coherent vortex[END_REF][START_REF] Vandermeirsch | Resistance of a coherent vortex to a vertical shear[END_REF]. [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF] have reported two distinct evolutions depending on the strength of the vortices compared to the shear. When the advection by the shear is below a critical value, the vortices in each layer are offset and rotate about each other. Conversely, when the advection is sufficiently strong, the two vortices are continuously advected away from each other. These two regimes have been called non-tearing and tearing regimes, respectively.

In contrast, the dynamics of an isolated vortex subjected to a vertical shear in stratified non-rotating fluids has been little investigated so far. [START_REF] Boulanger | Structure of a stratified tilted vortex[END_REF][START_REF] Boulanger | Tilt-induced instability of a stratified vortex[END_REF] have considered the dynamics of a vortex slightly tilted with respect to the direction of stratification. However, the tilt of the vortex is not created by an external shear flow but by the method of generation of the vortex. [START_REF] Majda | Model dynamics and vertical collapse in decaying strongly stratified flows[END_REF] have studied the effect of a vertically sheared mean horizontal flow on a vortex pair but only by means of the reduced equations of [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] which are valid in the limits of zero vertical and horizontal Froude numbers. Here and in chapter 5, we will study the dynamics of a Lamb-Oseen vortex in a sinusoidally sheared flow using the full Navier-Stokes equations under the Boussinesq approximation in the aim of testing Lilly's conjecture for the growth of vertical shear.

In the present chapter, two asymptotic analyses will be performed. The first analysis ( §4.3) is an asymptotic expansion for small times starting from the initial conditions. It will be shown that the initial response of the vortex is non-hydrostatic whatever the magnitude of the stratification. The second asymptotic analysis ( §4.4) is conducted for long vertical wavelength following a similar approach as used by [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 1. General stability equations[END_REF] to describe the unstable interaction between columnar vortices in stratified-rotating fluids. This analysis will provide the governing equations for the displacements of the vortex center and its structure evolution due to the shear flow. We shall see that the vortex is not only advected in the direction of the shear flow but also in the orthogonal direction owing to its self-induced motion. In addition, it will be demonstrated that the vertical vorticity of the vortex decreases with time owing to dynamic and viscous effects. The initial non-hydrostatic regime evidenced in the first asymptotic analysis will be related to the transient excitation of internal waves at the start-up of the motion.

These asymptotic predictions will be compared to Direct Numerical Simulations (DNS) in chapter 5.

Formulation of the problem

Governing equations

We consider the incompressible Navier-Stokes equations written within the Boussinesq approximation

∇ • u = 0, (4.3) ∂u ∂t + (u • ∇)u = -∇ p ρ 0 + be z + ν∇ 2 u, (4.4) ∂b ∂t + u • ∇b + N 2 w = κ∇ 2 b, (4.5)
with u = (u, v, w) being the velocity field in cartesian coordinates (x, y, z), p the pressure, b = -gρ/ρ 0 the buoyancy, g the gravity, e z the vertical unit vector oriented upwards, ν the viscosity, and κ the mass diffusivity. The total density field ρ t has been decomposed as ρ t (x, t) = ρ 0 + ρ(z) + ρ (x, t), ρ 0 being a constant reference density, ρ a linear mean density profile varying with the vertical coordinate z, and ρ (x, t) a perturbation density. The Brunt-Väisälä frequency measuring the ambient stratification N = -g/ρ 0 dρ/dz is assumed constant.

Initial condition

The initial flow is made of a columnar vortex embedded in a sinusoidal shear flow of the form U S = U (z)e x , where U (z) = ÛS sin( kz z), e x is the unit vector in the x direction, ÛS the amplitude and kz the vertical wavenumber. The vortex is chosen to have a Lamb-Oseen profile:

u v = Γ 2πr 1 -exp - r 2 a 2 0 e θ , (4.6)
where (r, θ, z) are cylindrical coordinates such that x = r cos(θ) and y = r sin(θ), Γ is the vortex circulation, a 0 its radius, and e θ the unit vector in the azimuthal direction. The total initial flow is therefore:

u = U S + u v . (4.7)

Non-dimensionalization

The equations (4.3-4.5) are non-dimensionalized by the characteristics of the vortex: the length and time units are taken as the vortex radius a 0 and the turnover time scale 2πa 2 0 /|Γ|. Keeping the same notation for the non-dimensional variables, the dimensionless equations read

∇ h • u h + ∂w ∂z = 0, (4.8) ∂u h ∂t + (u h • ∇ h ) u h + w ∂u h ∂z = -∇ h p + 1 Re ∇ 2 u h , (4.9) ∂w ∂t + u h • ∇ h w + w ∂w ∂z = - ∂p ∂z + b + 1 Re ∇ 2 w, (4.10) ∂b ∂t + u h • ∇ h b + w ∂b ∂z = - w F 2 h + 1 ReSc ∇ 2 b, (4.11)
where u h denotes the horizontal velocity. The Reynolds and Froude numbers are defined as

Re = |Γ| 2πν , F h = |Γ| 2πa 2 0 N , (4.12)
and Sc = ν/κ is the Schmidt number.

In non-dimensional form, the velocities of the shear flow and the vortex become (4.13) where U S = ÛS /(|Γ|/2πa 0 ) and k z = kz a 0 . Without loss of generality, the vortex circulation has been considered positive. Indeed, a vortex with a negative circulation can be transformed into a positive one by simply changing the cartesian frame (e x , e y , e z ) into (-e x , e y , -e z ).

U S = U S sin(k z z)e x , u v = 1 -exp(-r 2 ) r e θ ,

Initial dynamics of the vortex

Asymptotic problem

The initial dynamics of the vortex will be derived analytically from the equations (4.8-4.11) in the inviscid limit Re = +∞ by expanding all the variables q = (u, v, w, p, b) in power series of t:

q(x, y, z, t) = q 0 (x, y, z) + tq 1 (x, y, z) + t 2 q 2 (x, y, z) + t 3 q 3 (x, y, z) + O(t 4 ). (4.14)
Then, each power of t can be solved successively. In the following, the shear flow and the vortex will be assumed to have the arbitrary profiles U S = U (z)e x and u v = u θ (r)e θ as long as they do not need to be specified to solve the problem. The fluid domain is assumed unbounded in the horizontal direction.

Zeroth-order sub-problem

The zeroth-order problem reads

∇ • u 0 = 0, (4.15) u h1 + u 0 ∂u h0 ∂x + v 0 ∂u h0 ∂y + w 0 ∂u h0 ∂z = -∇ h p 0 , (4.16
)

w 1 + u 0 ∂w 0 ∂x + v 0 ∂w 0 ∂y + w 0 ∂w 0 ∂z = - ∂p 0 ∂z + b 0 , (4.17) b 1 + u 0 ∂b 0 ∂x + v 0 ∂b 0 ∂y + w 0 ∂b 0 ∂z = - w 0 F 2 h . (4.18)
As already explained in section 4.2.2, the initial flow is taken as the sum of the shear flow and the vortex velocity field, i.e.:

u 0 = -sin(θ)u θ (r) + U (z), v 0 = cos(θ)u θ (r), w 0 = 0, b 0 = 0. (4.19)
Although u 0 and v 0 are the horizontal velocity components in cartesian coordinates, it is more compact and convenient to express them in terms of the cylindrical coordinates (r, θ).

The solution (4.19) satisfies the divergence condition (4.15), and the density equation (4.18) immediately reduces to b 1 = 0. In order to find the remaining fields p 0 , u 1 , v 1 and w 1 , the divergence of (4.16)-(4.17) is taken. This yields

∇ 2 p 0 = 2Ω drΩ dr , (4.20)
since ∇ • u 1 = 0, where Ω = u θ /r is the vortex angular velocity. A solution is (4.21) which expresses the balance between the centrifugal force and the pressure gradient in the vortex. This balance is not affected by the presence of the shear flow at leading order because U is only a function of the vertical coordinate. Using (4.21), (4.16)-(4.17) yield the velocity field at first order:

p 0 = ξΩ 2 (ξ)dξ,
u 1 = sin(2θ) 2 r dΩ dr U (z), v 1 = - 1 2 cos(2θ)r dΩ dr + ζ U (z), w 1 = 0, (4.22)
with ζ the vertical vorticity:

ζ = 1 r dr 2 Ω dr . (4.23)
It is worth mentioning that (4.22) corresponds to a bending of the vortex as a whole without any deformation of its velocity distribution in the horizontal plane. Indeed, let us consider the translated vortex velocity

u v (x, t) = u v (x -∆x(z, t), y -∆y(z, t)), (4.24)
with ∆x and ∆y representing the cartesian displacements of the vortex center that are null at t = 0. For small ∆x and ∆y, we have: the velocity field at first order (4.22) is identical to the order O(t) of (4.25) with the displacements ∆x = U (z)t and ∆y = 0. It is important to stress that these results are valid for any vortex profile u θ (r) and any shear flow profile U (z).

u v (x, t) = u v (x, y) -∆x ∂u v ∂x -∆y ∂u v ∂y + ∆x 2 2 ∂ 2 u v ∂x 2 + ∆y 2 2 ∂ 2 u v ∂y 2 + ∆x∆y ∂ 2 u v ∂x∂y + ... ( 4

First-order sub-problem

The set of equations at first order reads:

∇ h • u h1 = 0, (4.27) 2u h2 + u 0 ∂u h1 ∂x + u 1 ∂u h0 ∂x + v 0 ∂u h1 ∂y + v 1 ∂u h0 ∂y = -∇ h p 1 , (4.28) 2w 2 = - ∂p 1 ∂z , (4.29) b 2 = 0, (4.30) r 0 1 2 3 4 β -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 (a) k z = 0 k z = π/6 k z = π/2 k z = π k z = 3π/2 r 0 1 2 3 4 k z β -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 (b) k z = π/2 k z = π k z = 3π/2
∇ 2 p 1 = -2 cos(θ) d dr Ω drΩ dr U (z). (4.31)
In order to solve this Poisson equation, we now have to specify the profile U (z). In the following, we consider the sinusoidal shear flow U (z) = U S sin(k z z). The case of a linear shear flow is investigated in appendix A. We look for p 1 in the form

p 1 = cos(θ)β(r)U (z). (4.32)
Then, (4.31) becomes:

d 2 β dr 2 + 1 r dβ dr -k 2 z + 1 r 2 β = - d dr 1 r dr 2 Ω 2 dr . (4.33)
The solution is:

β(r) = -rΩ 2 (r) + β b (r), (4.34) with β b (r) = k 2 z I 1 (k z r) +∞ r K 1 (k z ξ)ξ 2 Ω 2 (ξ)dξ + k 2 z K 1 (k z r) r 0 I 1 (k z ξ)ξ 2 Ω 2 (ξ)dξ, (4.35)
where it has been imposed that β is not singular at r = 0 and vanishes at infinity. I 1 and K 1 are the modified Bessel functions of order one of the first and second kind, respectively. Figure 4.2 displays the function β for different vertical wavenumbers k z in the case of the Lamb-Oseen profile (4.13). Starting from k z = 3π/2, which corresponds to a dimensional vertical wavelength slightly larger than the vortex radius a 0 , it is observed that β decreases as k z diminishes to 0 (figure 4.2a). In the limit k z = 0, we have β = -rΩ 2 . For finite wavenumbers π/2 k z 3π/2, figure 4.2b shows that k z β is almost independent of k z . These values cover wavelengths ranging approximately from the vortex radius to four radii. Finally, using (4.34), (4.28)-(4.29) give the second-order velocities:

u 2 = - 1 4 dβ b dr + β b (r) r + cos(2θ) dβ b dr - β b (r) r U S sin(k z z) + sin(θ) 2 cos(2θ) d dr Ω - r 2 dΩ dr - d dr r 2 dΩ dr U 2 S sin 2 (k z z), (4.36) v 2 = - sin(2θ) 4 dβ b dr - β b (r) r U S sin(k z z) + cos(θ) 2 cos(2θ) d dr r 2 dΩ dr -Ω + d dr Ω + ζ 2 U 2 S sin 2 (k z z), (4.37) w 2 = - cos(θ) 2 β(r)k z U S cos(k z z). (4.38)
In contrast to its counterpart at first order, this second order velocity field does not correspond in general to a bending of the vortex without internal deformation. This is the case only in the limit k z 1 for which β b can be neglected at leading order in (4.36)-(4.37) since

β b = k 2 z r 2 +∞ r ξΩ 2 (ξ)dξ + 1 2r r 0 ξ 3 Ω 2 (ξ)dξ + O(k 4 z ), (4.39)
i.e.

β b = O(k 2 z ).
The second order horizontal velocity field (4.36)-(4.37) is then equal to the order O(t 2 ) of (4.25) with the displacements ∆x = U (z)t and ∆y = 0 since

∂ 2 u v ∂x 2 = sin(θ) cos(2θ) d dr Ω -

Second-order sub-problem

At second order, we do not consider all the equations but only the buoyancy equation:

3b 3 = - w 2 F 2 h . (4.41)
This yields the first non-zero term in the expansion of the buoyancy:

b 3 = k z U S 6F 2 h cos(θ)β(r) cos(k z z). (4.42)
This shows that w = w 2 t 2 and b = b 3 t 3 at leading order. Using (4.41), the ratio of the vertical acceleration to the buoyancy is therefore

1 b Dw Dt = 6 F h t 2 . (4.43)
Hence, for small times such that t F h , the buoyancy is negligible relative to the vertical acceleration, meaning that the dynamics is initially non-hydrostatic. This transient non-hydrostatic phase lasts until t = O(F h ) whereas for t F h , a hydrostatic regime is expected to be recovered if F h < 1. This subsequent hydrostatic regime will occur for t 1 as assumed here, only if F h 1. We can further remark that the vertical velocity w 2 and the buoyancy b 3 are both proportional to k z β. From the behaviour of the function β displayed in figure 4.2, we can deduce that these two quantities are therefore proportional to k z for small k z whereas they are weakly sensitive to the variation of the vertical wavenumber for π/2 k z 3π/2.

Summary

The latter analysis shows that at first order in time, the vortex is simply advected as a whole in the x direction by ∆x = U (z)t. Internal deformation and generation of vertical motions occur only at order O(t 2 ). At this order, there is no displacement of the vortex implying that the higher orders in the displacements (∆x, ∆y) are at least O(t 3 ) or higher. This agrees with Jones' observation [START_REF] Jones | The evolution of vortices in vertical shear. I: Initially barotropic vortices[END_REF] in stratified rotating fluids that the initial effect of the vertical shear is a bending of the vortex in the streamwise direction only. We can further remark that the flow field at order O(t 2 ) does not depend on the Froude number but depends only on the vertical wavenumber k z of the sinusoidal shear flow. In other words, the stratification has no effect on the initial dynamics of the vortex. We have not determined the flow field at order O(t 3 ) but it is clear that it would also be independent of the Froude number. The stratification only intervenes in the buoyancy at order O(t 3 ). These results show the existence of an initial transient non-hydrostatic regime whose duration scales like the Froude number.

In the next section, we shall conduct an asymptotic analysis for long vertical wavelengths which is valid for longer times. This investigation will show that the initial non-hydrostatic phase is due to the generation of internal waves and that the vortex is also displaced in the y direction for finite times.

Long-wavelength dynamics of the vortex

Formulation of the problem

In section §4.3, we have studied the initial dynamics of the vortex by means of an asymptotic expansion in powers of t for any vertical wavenumber. The dynamics of the vortex can be also described asymptotically for finite time in the limit of small vertical wavenumber. Indeed, when k z = 0, i.e. in the two-dimensional limit, an exact solution is (4.44) which means that the vortex is simply advected at constant speed U . This solution is expected to remain valid at leading order when ε = k 2 z F 2 h is such that ε 1, i.e. when the vertical variations ∂/∂z of U (z) rescaled by the Froude number F h are weak:

u h = u h0 (x -U t, y) + U e x ,
√ ε = O F h ∂ ∂z 1. (4.45)
We consider that k z F h is small and not simply that k z 1 because of the self-similarity of strongly stratified fluids [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF]. Indeed, this self-similarity implies that the horizontal velocity depends on F h and k z only through their product

k z F h when F h
1. This means that the condition ε 1 can be fulfilled either by having k z 1 or F h 1. In each case, the other parameter F h or k z can be large provided that k z F h 1. We will keep the asymptotic analysis as general as possible in order to encompass both cases. Nevertheless, some parts of the analysis will assume F h 1 in order to simplify the calculations. Consistently with (4.45), we rescale the vertical coordinate and wavenumber as follows:

z = √ ε F h z, kz = k z F h √ ε , (4.46)
so that kz = O(1) and ∂/∂ z = O(1). Then, we can look for the horizontal velocity in the form of an asymptotic expansion:

u h = u h0 (x -U (z)t -δx, y -δy, η, τ ) + U (z)e x + εu h1 (x, y, z, t, η, τ ) + O(ε 2 ), (4.47)
where δx(z, η, τ ) and δy(z, η, τ ) are corrections to the displacements of the vortex center varying along the vertical and with the slow time scale τ = εt. In contrast to classical multiscale expansions, the fast and slow times t and τ should not be viewed as completely independent variables. We shall see indeed that the slow time evolution depends on the fast time. The introduction of τ is simply a convenient way to separate the dynamics at different orders1 . Similarly, a slow vertical length scale η = εz is introduced. We also define local horizontal cartesian coordinates (x, ỹ) centered on the vortex center at the level z:

x = x -U (z)t -δx(z, η, τ ), ỹ = y -δy(z, η, τ ). (4.48)
The cylindrical coordinates associated to these new coordinates are denoted with a tilde (x, ỹ) = (r cos( θ), r sin( θ)). In addition to the horizontal velocity, the other fields are expanded in the form:

p = p 0 (x, ỹ, η, τ ) + εp 1 (x, ỹ, z, t, η, τ ) + O(ε 2 ), (4.49) w = √ εF h w 0 (x, ỹ, z, t, η, τ ) + εw 1 (x, ỹ, z, t, η, τ ) + O(ε 2 ) , (4.50) b = √ ε F h b 0 (x, ỹ, z, t, η, τ ) + εb 1 (x, ỹ, z, t, η, τ ) + O(ε 2 ) , (4.51)
where the vertical velocity and the buoyancy have been scaled by √ εF h and √ ε/F h , respectively, for convenience. These scalings are equivalent to those of [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] for strongly stratified flows. Weak viscous effects will be also taken into account by assuming 1/Re = O(ε). Thus, we write:

1 Re = ε Re , (4.52)
where Re = O(1).

The expansions (4.47, 4.49, 4.50, 4.51) will be introduced in the governing equations (4.8-4.11). It will be useful to consider also the equation for the vertical vorticity ζ = ∂v/∂x -∂u/∂y:

∂ζ ∂t + u h • ∇ h ζ + w ∂ζ ∂z = ω h • ∇ h w + ζ ∂w ∂z + 1 Re ∇ 2 ζ, (4.53)
where ω h is the horizontal vorticity.

Because of the introduction of the local coordinates (x, ỹ) that depend on t, τ , z and η (see (4.48)), the time and vertical derivatives of any quantity Λ in (4.8-4.11) and (4.53) become:

∂Λ ∂t = ∂Λ ∂t x,ỹ,τ -U (z) ∂Λ ∂ x + ε   ∂Λ ∂τ x,ỹ,t - ∂δx ∂τ ∂Λ ∂ x - ∂δy ∂τ ∂Λ ∂ ỹ   , (4.54) ∂Λ ∂z = √ ε F h   ∂Λ ∂ z x,ỹ,η -t dU dz ∂Λ ∂ x - ∂δx ∂ z η ∂Λ ∂ x - ∂δy ∂ z η ∂Λ ∂ ỹ   + ε √ ε F h   ∂Λ ∂ η x,ỹ,z - ∂δx ∂ η z ∂Λ ∂ x - ∂δy ∂ η z ∂Λ ∂ ỹ   . (4.55)
It is important to remark the presence of the term tdU/dz in (4.55). Because of this term, the long-wavelength assumption (4.45) is expected to be no longer valid when √ εtdU/dz > 1, i.e. when k z F h U S t > 1 for the shear flow profile U = U S sin(k z z).

Leading order problem

At zeroth-order, equations (4.8-4.11) reduce to:

∇ h • u h0 = 0, (4.56) (u h0 • ∇ h ) u h0 = -∇ h p 0 , (4.57) F 2 h   ∂w 0 ∂t x,ỹ,τ + u h0 • ∇ h w 0   = - ∂p 0 ∂ z + b 0 , (4.58) ∂b 0 ∂t x,ỹ,τ + u h0 • ∇ h b 0 = -w 0 , (4.59)
where the operator ∇ h is with respect to x and ỹ, or equivalently r and θ. Note that the advection by the shear flow U (z)∂/∂ x is absent in (4.57-4.59) since it cancels with the term coming from the time derivative (see (4.54)). The equations (4.56)-(4.57) are the steady two-dimensional Euler equations. Consistently with the initial conditions chosen in §4.2.2, we take an axisymmetric vortex at leading order, namely u h0 = rΩ(r, η, τ )e θ. Then, (4.57) reduces to:

∂p 0 ∂r = rΩ 2 (r, η, τ ), (4.60)
which is the balance between the centrifugal force and the radial pressure gradient as if the shear flow were absent. Therefore, the vertical pressure gradient is given by

∂p 0 ∂ z = ∂r ∂ z ∂p 0 ∂r = -cos( θ) dU dz t + ∂δx ∂ z + sin( θ) ∂δy ∂ z rΩ 2 (4.61)
since we have at leading order:

∂r ∂ z = -cos( θ) dU dz t + ∂δx ∂ z -sin( θ) ∂δy ∂ z . (4.62)
Using (4.61), the equations (4.58) and (4.59) can be written as

F 2 h ∂w 0 ∂t + Ω ∂w 0 ∂ θ = b 0 + cos( θ) dU dz t + ∂δx ∂ z + sin( θ) ∂δy ∂ z rΩ 2 , (4.63) ∂b 0 ∂t + Ω ∂b 0 ∂ θ = -w 0 , (4.64)
where the fact that the time derivative ∂/∂t is taken at constant x, ỹ and τ and the vertical derivative ∂/∂ z at constant x, ỹ and η will be no longer mentioned to simplify the notation. To solve (4.63)-(4.64), w 0 and b 0 are decomposed into two components:

w 0 = w f (r, θ, z, t, η, τ ) + w s (r, θ, z, η, τ ), b 0 = b f (r, θ, z, t, η, τ ) + b s (r, θ, z, η, τ ), (4.65) that are solution of F 2 h ∂w f ∂t + Ω ∂w f ∂ θ = b f + cos( θ)rΩ 2 dU dz t, (4.66) ∂b f ∂t + Ω ∂b f ∂ θ = -w f , (4.67)
and

F 2 h Ω ∂w s ∂ θ = b s + cos( θ) ∂δx ∂ z + sin( θ) ∂δy ∂ z rΩ 2 , (4.68) Ω ∂b s ∂ θ = -w s . (4.69)
Deriving the analytic expressions of the slowly evolving components (w s , b s ) is straightforward:

w s = -sin( θ) ∂δx ∂ z + cos( θ) ∂δy ∂ z CΩ, (4.70) b s = -cos( θ) ∂δx ∂ z + sin( θ) ∂δy ∂ z C, (4.71)
where

C(r, η, τ ) = rΩ 2 1 -F 2 h Ω 2 (4.72)
has no critical point and is well-defined when F h < 1. Determining the fast components w f and b f is more difficult. They can be found in the form

w f = w p + w w , b f = b p + b w , (4.73)
where (w p , b p ) is a particular solution of (4.66)-(4.67):

w p = cos( θ)F -t sin( θ)Ω C dU dz , (4.74) b p = -2F h sin( θ)G + t cos( θ) C dU dz , (4.75) with F(r, η, τ ) = 1 + F 2 h Ω 2 1 -F 2 h Ω 2 , G(r, η, τ ) = F h Ω 1 -F 2 h Ω 2 , (4.76)
while (w w , b w ) is a homogeneous solution of (4.66)-(4.67) corresponding to internal waves oscillating at the non-dimensional frequency 1/F h :

w w = 1 F 2 h cos(Ωt -θ) C 1 (r, z, η, τ ) cos t F h + C 2 (r, z, η, τ ) sin t F h + 1 F 2 h sin(Ωt -θ) S 1 (r, z, η, τ ) cos t F h + S 2 (r, z, η, τ ) sin t F h , (4.77) b w = 1 F h cos(Ωt -θ) C 2 (r, z, η, τ ) cos t F h -C 1 (r, z, η, τ ) sin t F h + 1 F h sin(Ωt -θ) S 2 (r, z, η, τ ) cos t F h -S 1 (r, z, η, τ ) sin t F h , (4.78)
where C 1 , C 2 , S 1 and S 2 are arbitrary functions. They are deduced by imposing that the vertical velocity and the buoyancy are null at t = τ = 0 since the vortex is strictly vertical at the start-up of the evolution. This implies that δx(z, η, 0) = δy(z, η, 0) = 0 so that w s and b s are zero at τ = 0. Therefore, we have to impose w f = b f = 0 at t = τ = 0, yielding:

C 1 = -F 2 h CF dU dz , S 1 = 0, C 2 = 0, S 2 = -2F 2 h CG dU dz . (4.79)
Strictly speaking, C 1 and S 2 are fixed at τ = 0 so that we should set C 1 = -F 2 h C(r, η, τ = 0)F(r, η, τ = 0)dU/dz and similarly for S 2 . However, it is also legitimate to impose (4.79) whatever τ since the evolution over the slow time scale τ is a higher order effect. Finally, collecting (

w 0 = cos( θ)F -t sin( θ)Ω -F cos(Ωt -θ) cos t F h -2G sin(Ωt -θ) sin t F h C dU dz + -sin( θ) ∂δx ∂ z + cos( θ) ∂δy ∂ z CΩ, (4.80) b 0 = -F h 2 sin( θ)G + t F h cos( θ) + 2G sin(Ωt -θ) cos t F h -F cos(Ωt -θ) sin t F h C dU dz -cos( θ) ∂δx ∂ z + sin( θ) ∂δy ∂ z C. (4.81)
The vertical velocity and buoyancy (4.80)-(4.81) contain three types of terms evolving on three different time scales:

1. the terms corresponding to w p and b p that grow linearly with the fast time t in response to the shear flow U (z), 2. internal waves oscillating at the frequency 1/F h (i.e. the Brunt-Väisälä frequency in dimensional units) that are generated at the onset of the motion, 3. the terms coming from w s and b s due to the slow evolution of the displacement corrections δx and δy.

In

√ ε F h b 0 = - 1 6F 2 h cos(θ)rΩ 2 dU dz t 3 + O(t 4 , τ ), (4.83)
since δx and δy are at most O(τ ) = O(εt) for τ 1 because δx(z, η, 0) = δy(z, η, 0) = 0. When U = U S sin( kz z), (4.82) and (4.83) are identical to the expressions (4.38) and (4.42) obtained for small time in §4.3 in the limit k z 1 since then β = -rΩ 2 + O(k 2 z ). It can be concluded that the initial transient non-hydrostatic regime found in §4.3 comes from the generation of internal waves of frequency 1/F h . Indeed, if the waves (4.77)-(4.78) were absent in (4.80)-(4.81), the vertical velocity and buoyancy would mostly grow linearly with time w ∝ F 2 h t and b ∝ t according to (4.74)-(4.75). These scaling laws fulfill the typical scaling associated to the hydrostatic balance, w = O(F 2 h b) in contrast to the full solution (4.80)-(4.81). For large time, the amplitude of the waves w w , b w will become small relative to w p , b p that grow linearly with time. Hence, the hydrostatic balance will be recovered. In addition, (4.80)-(4.81) show that the waves are advected differentially by the azimuthal motion in the vortex (terms cos(Ωtθ) and sin(Ωtθ)). We will see in §4.4.4 that their influence quickly dies out for this reason.

First-order problem

At first order, it is more convenient to consider the divergence equation and the equation for the vertical vorticity

ζ 1 = ∂v 1 /∂ x -∂u 1 /∂ ỹ: ∇ h • u h1 = - ∂w 0 ∂ z + ∂r ∂ z ∂w 0 ∂r + ∂ θ ∂ z ∂w 0 ∂ θ , (4.84) ∂ζ 1 ∂t + u h0 • ∇ h ζ 1 + u h1 • ∇ h ζ 0 = - ∂ζ 0 ∂τ + cos( θ) ∂δx ∂τ ∂ζ 0 ∂r + sin( θ) ∂δy ∂τ ∂ζ 0 ∂r + dU dz + ∂r ∂ z ∂u 0 ∂r + ∂ θ ∂ z ∂u 0 ∂ θ ∂w 0 ∂ ỹ - ∂r ∂ z ∂v 0 ∂r + ∂ θ ∂ z ∂v 0 ∂ θ ∂w 0 ∂ x + ζ 0 ∂w 0 ∂ z + ∂r ∂ z ∂w 0 ∂r + ∂ θ ∂ z ∂w 0 ∂ θ -w 0 ∂r ∂ z ∂ζ 0 ∂r + ∇ 2 ζ 0 Re , ( 4.85)
where ζ 0 = 2Ω + r∂Ω/∂r and

∂ θ ∂ z = sin( θ) r dU dz t + ∂δx ∂ z - cos( θ) r ∂δy ∂ z (4.86)
at leading order while ∂r/∂ z has been given previously in (4.62). The forcing terms in the left-hand side of (4.85) correspond to the slow evolution of ζ 0 and the displacements δx, δy (first row), the tilting of horizontal vorticity into vertical vorticity (second row), the stretching and vertical advection of vertical vorticity (third row), and viscous effects (last row). The viscous term in (4.85) reads

∇ 2 ζ 0 = ∂ 2 ζ 0 ∂r 2 + 1 r ∂ζ 0 ∂r + ε F 2 h ∂ 2 ζ 0 ∂ z2 . (4.87)
This expression suggests that we should neglect the second term since ε 1. However, we will keep it for two reasons: first, the condition ε 1 can be fulfilled either by having k z 1 or F h 1. In the latter case, the parameter ε/F2 h = k 2 z can be large since a wavenumber k z of order unity or larger can satisfy the condition k z F h 1. Secondly, because the vertical shear grows algebraically with time as shown by (4.62) and (4.86), the second term can become large when ε/F 2 h (tdU/dz) 2 > 1 2 . Actually, this is the main reason why this term may have an important effect. Hence, we will take into account only the terms proportional to tdU/dz in (4.62)-(4.86) and neglect the others to express the viscous term (4.87). This yields:

∇ 2 ζ 0 = ∂ 2 ζ 0 ∂r 2 + 1 r ∂ζ 0 ∂r + ε 2F 2 h ∂ 2 ζ 0 ∂r 2 + 1 r ∂ζ 0 ∂r dU dz 2 t 2 - ε F 2 h cos( θ) ∂ζ 0 ∂r d 2 U dz 2 t + ε 2F 2 h cos(2 θ) ∂ 2 ζ 0 ∂r 2 - 1 r ∂ζ 0 ∂r dU dz 2 t 2 . (4.88)
We emphasize again that the term O(ε/F 2 h ) is taken into account in order to encompass both cases k z 1 or F h 1. However, in the former case, this term is expected to be negligible as long as k z tdU/dz 1. We see that the viscous term contains azimuthal modes m = 0, 1, and 2. In fact, the other forcing terms in (4.84) and (4.85) can be also decomposed into these three azimuthal modes. Looking for example at the divergence equation (4.84), we see two categories of forcing terms. First, the term ∂w 0 /∂ z will produce terms proportional to cos( θ) or sin( θ). Secondly, the terms ∂r/∂ z∂w 0 /∂r and ∂ θ/∂ z∂w 0 /∂ θ will give azimuthal modes of the form cos(m θ) and sin(m θ) with m = 0 or m = 2. The righthand side of (4.85) can be decomposed similarly into the azimuthal modes m = 0, m = 1, and m = 2. For this reason, (4.84) and (4.85) can be re-written as follows:

∇ h • u h1 = - ∂w 0 ∂ z + δD 0 + δD 2 , (4.89) ∂ζ 1 ∂t + Ω ∂ζ 1 ∂ θ + u h1 • ∇ h ζ 0 = cos( θ) ∂δx ∂τ + sin( θ) ∂δy ∂τ ∂ζ 0 ∂r + ζ 0 ∂w 0 ∂ z - ε ReF 2 h cos( θ) ∂ζ 0 ∂r d 2 U dz 2 t + δZ 0 + δZ 2 , (4.90)
where δZ 0 and δD 0 correspond to the axisymmetric forcing terms, while δZ 2 and δD 2 are the forcing terms proportional to cos(2 θ) or sin(2 θ). In contrast, the forcing terms written explicitly in the right-hand sides of (4.89) and (4.90) are all proportional to cos( θ) or sin( θ).

In order to solve (4.89) and (4.90), the horizontal velocity is decomposed into potential and rotational components by introducing a potential Φ 1 and a streamfunction ψ 1 :

u h1 = ∇ h Φ 1 -∇ × ψ 1 e z . (4.91)
The potential and streamfunction are further decomposed into three components

Φ 1 = Φ 10 + Φ 11 + Φ 12 , ψ 1 = ψ 10 + ψ 11 + ψ 12 , (4.92)
corresponding to the azimuthal modes m = 0, m = 1, and m = 2, respectively. They can be solved separately. We will begin by the components (Φ 11 , ψ 11 ) which will provide the governing equations for the displacements (δx, δy). The other components will be next determined in section 4.4.5.

Determination of Φ 11 and ψ 11

We first determine the potential from (4.89):

∇ 2 h Φ 11 = - ∂w 0 ∂ z . (4.93)
Like the vertical velocity (4.80), Φ 1 is separated into three potentials evolving with the different time scales τ , t, and t/F h (corresponding to internal waves):

Φ 11 = Φ s (r, θ, z, η, τ ) + Φ p (r, θ, z, t, η, τ ) + Φ w (r, θ, z, t, η, τ ), (4.94) with ∇ 2 h Φ s = - ∂w s ∂ z , ∇ 2 h Φ p = - ∂w p ∂ z , ∇ 2 h Φ w = - ∂w w ∂ z , (4.95)
where w s is defined in (4.70), w p in (4.74), and w w in (4.77). The potential varying on the slow time scale reads

Φ s = 1 2 sin( θ) ∂ 2 δx ∂ z2 -cos( θ) ∂ 2 δy ∂ z2 Φ ps , (4.96) with Φ ps (r, η, τ ) = r r +∞ Ω(ξ, η, τ )C(ξ, η, τ )dξ - 1 r r 0 ξ 2 Ω(ξ, η, τ )C(ξ, η, τ )dξ, (4.97)
where it has been imposed that Φ s is not singular at r = 0 and vanishes as r tends to infinity. The two other potentials Φ p and Φ w can be derived similarly:

Φ p = 1 2 -cos( θ)Φ pc + t sin( θ)Φ ps d 2 U dz 2 , (4.98) Φ w = 1 2 r r +∞ C(ξ, η, τ )F(ξ, η, τ ) cos(Ω(ξ, η, τ )t -θ)dξ - 1 r r 0 ξ 2 C(ξ, η, τ )F(ξ, η, τ ) cos(Ω(ξ, η, τ )t -θ)dξ cos t F h d 2 U dz 2 + r r +∞ C(ξ, η, τ )G(ξ, η, τ ) sin(Ω(ξ, η, τ )t -θ)dξ - 1 r r 0 ξ 2 C(ξ, η, τ )G(ξ, η, τ ) sin(Ω(ξ, η, τ )t -θ)dξ sin t F h d 2 U dz 2 , (4.99)
where

Φ pc (r, η, τ ) = r r +∞ C(ξ, η, τ )F(ξ, η, τ )dξ - 1 r r 0 ξ 2 C(ξ, η, τ )F(ξ, η, τ )dξ. (4.100)
The potential Φ 11 being determined, an equation for the streamfunction ψ 11 is obtained from the equation (4.90) for the vertical vorticity:

∂ ∂t + Ω ∂ ∂ θ ∇ 2 h ψ 11 - 1 r ∂ζ 0 ∂r ∂ψ 11 ∂ θ = cos( θ) ∂δx ∂τ + sin( θ) ∂δy ∂τ - ∂Φ 11 ∂r ∂ζ 0 ∂r + ζ 0 ∂w 0 ∂ z - ε ReF 2 h cos( θ) ∂ζ 0 ∂r d 2 U dz 2 t. (4.101)
The first two terms in the left-hand side of (4.101) come from the slow evolution of δx and δy. The third term is the advection of ζ 0 by the first order potential flow. The fourth term is the stretching of ζ 0 . The last term corresponds to the dominant viscous effect. In order to solve (4.101), the streamfunction is decomposed as

ψ 11 = ψ sd (r, θ, z, τ ) + ψ s (r, θ, z, η, τ ) + ψ p (r, θ, z, t, η, τ ) + ψ w (r, θ, z, t, η, τ ), (4.102)
where the different streamfunctions are solution of:

Ω ∂∇ 2 h ψ sd ∂ θ - 1 r ∂ζ 0 ∂r ∂ψ sd ∂ θ = cos( θ) ∂δx ∂τ - εt ReF 2 h d 2 U dz 2 + sin( θ) ∂δy ∂τ ∂ζ 0 ∂r , (4.103) Ω ∂∇ 2 h ψ s ∂ θ - 1 r ∂ζ 0 ∂r ∂ψ s ∂ θ = - ∂ζ 0 ∂r ∂Φ s ∂r -ζ 0 ∇ 2 h Φ s , (4.104) ∂ ∂t + Ω ∂ ∂ θ ∇ 2 h ψ p - 1 r ∂ζ 0 ∂r ∂ψ p ∂ θ = - ∂ζ 0 ∂r ∂Φ p ∂r -ζ 0 ∇ 2 h Φ p , (4.105) ∂ ∂t + Ω ∂ ∂ θ ∇ 2 h ψ w - 1 r ∂ζ 0 ∂r ∂ψ w ∂ θ = - ∂ζ 0 ∂r ∂Φ w ∂r -ζ 0 ∇ 2 h Φ w , (4.106)
where Φ s , Φ p and Φ w are given in (4.96), (4.98), and (4.99), respectively. The solutions of (4.103) and (4.104) which are non-singular are

ψ sd = -sin( θ) ∂δx ∂τ - εt ReF 2 h d 2 U dz 2 + cos( θ)
∂δy ∂τ (r + µrΩ) , (4.107)

ψ s = cos( θ) ∂ 2 δx ∂ z2 + sin( θ) ∂ 2 δy ∂ z2 ψs + αrΩ , (4.108)
where µ and α are arbitrary functions of η and τ and ψs (r, η, τ ) = rΩ(r, η, τ )

r 0 1 η 3 Ω 2 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξdη + rΩ(r, η, τ ) r +∞ C(ξ, η, τ )dξ - r 2 ∂Φ ps ∂r .
(4.109)

The functions µ, α are arbitrary because rΩ cos( θ) and rΩ sin( θ) are homogeneous solutions of (4.103) and (4.104). These solutions derive from the translational invariance [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 1. General stability equations[END_REF]. In order to set µ and α, we need to use a normalisation condition.

Here, we have chosen to impose that the mean displacement associated to the first order vorticity ∇ 2 h ψ 11 is zero, i.e. In this way, only the displacements U (z)t + δx and δy of the leading order vorticity ζ 0 will contribute to the mean displacement of the vortex. This will ease the comparison with the numerical simulations performed in chapter 5 where the vortex displacements will be computed from vertical vorticity centroids. Hence, to enforce (4.110), we impose

∇ 2 h ψ sd , r = 0, (4.111) ∇ 2 h ψ s , r = 0, (4.112)
where the scalar product is defined as

< f, g >= +∞ 0 f grdr (4.113)
for any functions f and g. Equations (4.111) and (4.112) lead to

µ = 0, α = - ∇ 2 r ψs , r ∂ζ 0 ∂r , r , (4.114)
where the operator ∇ 2 r is given by

∇ 2 r Λ = ∂ 2 Λ ∂r 2 + 1 r ∂Λ ∂r - Λ r2 (4.115)
for any function Λ. Determining ψ p and ψ w is more difficult because of the dependence of the forcings on the fast time t in (4.105)-(4.106). We first derive ψ p analytically. Since the potential Φ p evolves linearly with t (see (4.98)), we seek a solution in the form

ψ p = -sin( θ) (λ 0 rΩ + λ 1 r + ψ ps (r, η, τ )) + t cos( θ) (λ 1 rΩ + ψ pc (r, η, τ )) d 2 U dz 2 , ( 4
.116) where λ 0 and λ 1 are arbitrary functions of η and τ . This leads to the following equations for ψ ps and ψ pc :

Ω∇ 2 r ψ pc - 1 r ∂ζ 0 ∂r ψ pc = - 1 2 ∂ζ 0 ∂r ∂Φ ps ∂r + ζ 0 ∇ 2 r Φ ps , (4.117) Ω∇ 2 r ψ ps - 1 r ∂ζ 0 ∂r ψ ps = - 1 2 ∂ζ 0 ∂r ∂Φ pc ∂r + ζ 0 ∇ 2 r Φ pc -∇ 2 r ψ pc , (4.118)
The solutions that are non-singular at r = 0 are:

ψ pc = -ψs , (4.119) ψ ps (r, η, τ ) = r 2 ∂Φ pc ∂r -r r +∞ C(ξ, η, τ )dξ -2F h rΩ r 0 C(ξ, η, τ )G(ξ, η, τ )dξ -2F h rΩ r 0 1 η 3 Ω 2 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )G(ξ, η, τ )dξdη + 2rΩ r 0 1 η 3 Ω 3 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξdη -r r 0 1 η 3 Ω 2 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξdη, (4.120)
where ψs is defined in (4.109). The two arbitrary functions λ 0 and λ 1 in (4.116) come from the homogeneous solution of (4.105) which is of the form

ψ h = rΩ cos( θ)f + sin( θ)g + r sin( θ) ∂f ∂t -cos( θ) ∂g ∂t , (4.121)
where f and g are arbitrary functions of z, t, η and τ . In (4.116), we have taken f (z, t, η, τ ) = -tλ 1 (η, τ )d 2 U/dz 2 and g(z, η, τ ) = -λ 0 (η, τ )d 2 U/dz 2 in accordance with the time dependence of this equation. As before, we impose the normalisation conditions

∇ 2 r ψ ps , r + λ 0 (η, τ ) ∂ζ 0 ∂r , r = 0, (4.122) ∇ 2 r ψ pc , r + λ 1 (η, τ ) ∂ζ 0 ∂r , r = 0, (4.123)
which enforce that the mean displacement associated to ∇ 2 h ψ p is zero. This sets λ 0 and λ 1 to:

λ 0 = - ∇ 2 r ψ ps , r ∂ζ 0 ∂r , r , λ 1 = - ∇ 2 r ψ pc , r ∂ζ 0 ∂r , r . (4.124)
Similarly, the solution of (4.106) is sought in the form where the potentials Φ wc and Φ ws are given by:

ψ w = -cos( θ) ψ wc (
Φ wc (r, t, η, τ ) = 1 2 1 r r 0 ξ 2 C(ξ, η, τ )F(ξ, η, τ ) cos(Ω(ξ, η, τ )t)dξ -r r +∞ C(ξ, η, τ )F(ξ, η, τ ) cos(Ω(ξ, η, τ )t)dξ cos t F h + 1 r r 0 ξ 2 C(ξ, η, τ )G(ξ, η, τ ) sin(Ω(ξ, η, τ )t)dξ -r r +∞ C(ξ, η, τ )G(ξ, η, τ ) sin(Ω(ξ, η, τ )t)dξ sin t F h , (4.128) Φ ws (r, t, η, τ ) = 1 2 1 r r 0 ξ 2 C(ξ, η, τ )F(ξ, η, τ ) sin(Ω(ξ, η, τ )t)dξ -r r +∞ C(ξ, η, τ )F(ξ, η, τ ) sin(Ω(ξ, η, τ )t)dξ cos t F h - 1 r r 0 ξ 2 C(ξ, η, τ )G(ξ, η, τ ) cos(Ω(ξ, η, τ )t)dξ -r r +∞ C(ξ, η, τ )G(ξ, η, τ ) cos(Ω(ξ, η, τ )t)dξ sin t F h . (4.129)
Because the time and radial dependencies in Φ wc and Φ ws are not separated, it is not possible to solve (4.126) and (4.127) analytically. We have therefore implemented a numerical resolution of (4.126)-(4.127). The time advancement is performed by a forward Euler scheme, spatial derivatives are computed by centered finite differences, and the Laplacian operator is inverted by means of the tridiagonal algorithm. The imposed boundary conditions are ψ ws = ψ wc = 0 at r = 0 and that ψ ws and ψ wc vanish as r → +∞. Since the total streamfunction ψ 11 should be zero at t = τ = 0, the initial condition should be ψ w = -ψ p because ψ s (r, θ, z, η, τ = 0) = 0 and ψ sd (r, θ, z, τ = 0) = 0 since δx = δy = ∂δx/∂τ = ∂δy/∂τ = 0 at τ = 0. Thus, we impose at t = τ = 0:

ψ wc = 0, ψ ws = -[ψ ps + λ 0 rΩ + λ 1 r -ϑ(η, τ = 0)r] , (4.130)
where ϑ(η, τ = 0) is set such that ψ ws → 0 as r → +∞ initially, as required to integrate (4.126) and (4.127) numerically. From the expression (4.120) of ψ ps , we have:

ϑ(η, τ ) = λ 1 (η, τ ) -F h +∞ 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )G(ξ, η, τ )dξ + lim R→+∞ 2Ω( R, η, τ ) R 0 1 η 3 Ω 3 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξdη -lim R→+∞ R 0 1 η 3 Ω 2 (η, η, τ ) η 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξdη. (4.131)
There is no loss of generality in imposing that ψ ws vanishes as r → +∞ initially because the function ϑr belongs to the homogeneous solution of (4.106):

ψ wh = ϑt cos( θ)rΩ + ϑr sin( θ). (4.132)
In other words, introducing the term ϑr is equivalent to substituting C w by C w + ϑt in (4.125). After the numerical integration, we enforce that the mean displacement associated to ∇ 2 h ψ w is zero, yielding:

C w (t, η, τ ) = - ∇ 2 r ψ wc , r ∂ζ 0 ∂r , r , S w (t, η, τ ) = - ∇ 2 r ψ ws , r ∂ζ 0 ∂r , r . (4.133)
Figures 4.3 and 4.4 show the streamfunctions ψ wc and ψ ws plotted against r at different times, for F h = 0.1 and F h = 0.5 respectively, at z = l z /4. In the next section, we will show that the angular velocity Ω does not evolve with τ in this plane. After an initial transient phase, the streamfunctions globally decay with time in the vortex core.

We have now fully determined the four components of the streamfunction at first order ψ 11 . They all satisfy the boundary condition at r = 0. However, for large r, we see directly in (4.107) that ψ sd grows like r. Similarly, it can be shown that the three other components behave as where A is given by:

ψ s ∼ A 2 cos( θ) ∂ 2 δx ∂ z2 + sin( θ) ∂ 2 δy ∂ z2
A(η, τ ) = +∞ 0 ξ 2 Ω 2 (ξ, η, τ )C(ξ, η, τ )dξ. (4.137)
Note that we have used the fact that ψ wc and ψ ws tend to zero as r → +∞. In order to enforce the boundary condition ψ 11 → 0 as r → +∞, the sum of the four components should be zero, giving the governing equations for (δx, δy):

∂δx ∂τ = A(η, τ ) 2 ∂ 2 δy ∂ z2 -ϑ(η, τ ) + ∂C w ∂t (t, η, τ ) d 2 U dz 2 + εt ReF 2 h d 2 U dz 2 , (4.138) ∂δy ∂τ = - A(η, τ ) 2 t d 2 U dz 2 + ∂ 2 δx ∂ z2 - ∂S w ∂t (t, η, τ ) d 2 U dz 2 . (4.139)
When U = 0, these equations are identical to those found by [START_REF] Billant | Zigzag instability of vortex pairs in stratified and rotating fluids. Part 1. General stability equations[END_REF] for the self-induction of an isolated vortex in stratified non-rotating fluids. Figure 4.5 shows A and ϑ as functions of F h at z = l z /4, where the angular velocity remains the same Ω = [1 -exp(-r 2 )] /r 2 whatever τ (see §4.4.6). They are almost independent of F h for F h 0.5. The wave forcing terms C w and S w are also shown in figure 4.6 for two Froude numbers F h = 0.1 and F h = 0.5, at the same location. In appendix B, it is shown that they behave for small time as

C w = -ϑt + O(t 3 ), S w = - A 4 t 2 + O(t 4 ). (4.140)
Then, they go through oscillations as seen in figure 4.6. These oscillations exhibit two typical periods: the period of rotation in the vortex core T ∼ 2π and the buoyancy period T ∼ 2πF h . For F h = 0.5 (figure 4.6b), these two periods are comparable so that C w and S w exhibit irregular oscillations while for F h = 0.1 (figure 4.6a), fast oscillations superimposed on slow oscillations can be clearly distinguished. For large times, C w and S w tend to constant values because the waves are increasingly sheared by the vortex rotation through the terms cos(Ωtθ) and sin(Ωtθ) in Φ w (see (4.99)). As a result, their contributions ∂C w /∂t and ∂S w /∂t to the vortex displacements vanish as t → +∞.

In order to be able to solve (4.138)-(4.139) in general, we first need to know how the angular velocity evolves on the slow time τ . This will be determined in the next two sections while further analyses of (4.138)-(4.139) will be conducted in §4.4.7. 4.4.5 Determination of (Φ 10 , ψ 10 ) and (Φ 12 , ψ 12 )

We now determine (Φ 10 , ψ 10 ) and (Φ 12 , ψ 12 ) corresponding to the azimuthal modes m = 0 and m = 2 of the potential and streamfunction at first order (4.92).

In order to keep the calculations relatively simple, we will assume here that the Froude number is small so that the hydrostatic approximation is satisfied. The comparison with the DNS in chapter 5 will show that this assumption provides accurate predictions for F h 1. In the limit F h → 0, the vertical velocity and buoyancy

w 0 = cos( θ) -t sin( θ)Ω -cos(Ωt -θ) cos t F h rΩ 2 dU dz + -sin( θ) ∂δx ∂ z + cos( θ) ∂δy ∂ z rΩ 3 , (4.141) b 0 = -t cos( θ)rΩ 2 dU dz -cos( θ) ∂δx ∂ z + sin( θ) ∂δy ∂ z rΩ 2 . (4.142)
The only remaining term due to internal waves is the third term in (4.141). However, if the Froude number F h were set to zero directly in (4.63)-(4.64), this term would be absent in (4.141). Furthermore, the previous section has shown that the effects due to these internal waves are restricted to the beginning of the evolution and are rather weak (see also chapter 5). For this reason, we will neglect the third term in (4.141) in the following. This will avoid lengthy calculations. In addition, the equations (4.138)-(4.139) for the displacement perturbations derived in the previous section show that (δx, δy) remain small. Therefore, these displacements will be also neglected in (4.141)-(4.142). Hence, (4.141)-(4.142) are simplified as follows: In (4.147), we have neglected the horizontal dissipation (first two terms in the righthand side of (4.88)) and retained only the terms coming from the vertical dissipation (third term in the right-hand side of (4.88)) because the latter term grows like t 2 in contrast to the former.

w 0 = cos( θ) -t sin( θ)Ω rΩ 2 dU dz , (4.143) b 0 = -t cos( θ)rΩ 2 dU dz . ( 4
We first determine the potentials Φ 10 and Φ 12 . They satisfy (see (4.89)):

∇ 2 h Φ 10 = δD 0 , ∇ 2 h Φ 12 = δD 2 . (4.149)
The solutions that are not singular at r = 0 and that vanish at infinity are

Φ 10 = t 2 r +∞ ξΩ 2 (ξ, η, τ )dξ dU dz 2 , (4.150) Φ 12 = t cos(2 θ)Φ c2 - 3 2 t 2 sin(2 θ)Φ s2 dU dz 2 , ( 4.151)
where:

Φ c2 (r, η, τ ) = 1 2r 2 r 0 ξ 3 Ω 2 (ξ, η, τ )dξ, Φ s2 (r, η, τ ) = 1 3r 2 r 0 ξ 3 Ω 3 (ξ, η, τ )dξ. (4.152)
Then, we look for the streamfunctions ψ 10 and ψ 12 . Their equations are directly obtained from (4.90). Using (4.147), (4.148), (4.150) and (4.151): in order that the vertical vorticity at order 1 is orthogonal to the vertical vorticity at leading order. Taking the scalar product of (4.153) with

∂ ∂t + Ω ∂ ∂ θ ∇ 2 h ψ 10 - 1 r ∂ζ 0 ∂r ∂ψ 10 ∂ θ = - ∂ζ 0 ∂τ -2tΩ 3 + 3trΩ 2 ∂Ω ∂r dU dz 2 + εt 2 2 ReF 2 h ∂ 2 ζ 0 ∂r 2 + 1 r ∂ζ 0 ∂r dU dz 2 , (4.153) ∂ ∂t + Ω ∂ ∂ θ ∇ 2 h ψ 12 - 1 r ∂ζ 0 ∂r ∂ψ 12 ∂ θ = εt 2 2 ReF 2 h cos(2 θ) ∂ 2 ζ 0 ∂r 2 - 1 r
ζ 0 shows that ∂ ∇ 2 h ψ 10 , ζ 0 /∂t = 0 only if ∂ζ 0 ∂τ = -2tΩ 3 + 3trΩ 2 ∂Ω ∂r dU dz 2 + εt 2 2 ReF 2 h ∂ 2 ζ 0 ∂r 2 + 1 r ∂ζ 0 ∂r dU dz 2 . (4.156)
The streamfunction ψ 12 will be sought in the form

ψ 12 = ψ * 12 + ψ ν 12 , (4.157)
where ψ ν 12 is the solution forced by the viscous term in the right-hand side of (4.154) (first term) while ψ * 12 is the solution forced by all the other terms. We will first determine ψ * 12 in the form

ψ * 12 = 1 2 -cos(2 θ)ψ 2b + t sin(2 θ)ψ 2s -t 2 cos(2 θ)ψ 2c dU dz 2 , (4.158)
where the streamfunctions ψ 2b (r, η, τ ), ψ 2s (r, η, τ ), and ψ 2c (r, η, τ ) are solution of

∇ 2 * r ψ 2s 2 + Ω∇ 2 * r ψ 2b - 1 r ∂ζ 0 ∂r ψ 2b = rΩ ∂Ω ∂r , (4.159) -∇ 2 * r ψ 2c + Ω∇ 2 * r ψ 2s - 1 r ∂ζ 0 ∂r ψ 2s = rΩ 2 ∂Ω ∂r + 1 r3 ∂ζ 0 ∂r r 0 ξ 3 Ω 2 (ξ, η, τ )dξ, (4.160) Ω∇ 2 * r ψ 2c - 1 r ∂ζ 0 ∂r ψ 2c = rΩ 3 ∂Ω ∂r - 1 r3 ∂ζ 0 ∂r r 0 ξ 3 Ω 3 (ξ, η, τ )dξ, (4.161)
where ∇ 2 * r is defined as:

∇ 2 * r ϕ = ∂ 2 ϕ ∂r 2 + 1 r ∂ϕ ∂r - 4ϕ r2 . (4.162)
It can be remarked that the decomposition (4.158) contains a constant term implying that ψ * 12 does not vanish at t = 0. Yet, the total streamfunction at first order ψ 1 (4.92) should be zero at t = 0 since there is no perturbation initially. Because ψ 10 and ψ 11 already fulfill this condition, it should be the case also for ψ * 12 . However, the simplified expression of the vertical velocity (4.143), which is at the origin of the inviscid forcing terms in (4.154) (see (4.85)), does not vanish either at t = 0. This is a direct consequence of the fact that the term cos(Ωtθ) cos(t/F h ), corresponding to internal waves, has been neglected in (4.143). Nevertheless, if this term had been kept, the vertical velocity would be zero at t = 0 and ψ 1 would contain an additional contribution ψ w 12 due to this wave term. This contribution could be computed so that it cancels ψ * 12 at t = 0, i.e. ψ * 12 + ψ w 12 = 0 initially. This is the procedure that has been followed to compute ψ 11 in §4.4.4. However, the contribution due to the internal waves rapidly decays as seen in figure 4.6. This is the reason why we have chosen to neglect this term here so as to simplify the calculations. However, this implies that w 0 and ψ * 12 are not exactly zero at t = 0.

In (4.159-4.161), the angular velocity Ω and vertical vorticity ζ 0 evolve on the slow time according to (4.156). In the next section, it will be shown that this evolution can be very well accounted for by assuming that the vortex keeps a Lamb-Oseen profile but with time varying maximum angular velocity Ω c (η, τ ) and radius ãe (η, τ ):

Ω = Ω c (η, τ ) Ω(r), Ω = 1 -exp(-r 2 ) r2 , ( 4.163)
where r = r/ã e . Then, by rescaling ψ 2b , ψ 2s , and ψ 2c as follows: Note that they all rapidly vanish outside the vortex core.

ψ 2b = ã2 e Ω c ψ2b (r), ψ 2s = ã2 e Ω 2 c ψ2s (r), ψ 2c = ã2 e Ω 3 c ψ2c (r), ( 4
These three equations share the same linear operator which admits a homogeneous solution behaving as ψ 2h = r2 for r → +∞ and ψ 2h = 2.524r 2 for r → 0 [START_REF] Moore | The instability of a straight vortex filament in a strain field[END_REF][START_REF] Le Dizès | Theoretical predictions for the elliptical instability in a two-vortex flow[END_REF]. By using this homogeneous solution, it is possible to determine numerically complete solutions that vanish at infinity and are non-singular at r = 0. They are shown in figure 4.7. For r 1, they behave as

ψ2b = β 2 r2 + O(r 4 ), ψ2s = σ 2 r2 + O(r 4 ), ψ2c = χ 2 r2 + O(r 4 ), (4.168)
with:

β 2 = -0.1506, σ 2 = 0.3786, χ 2 = -9.830 × 10 -3 . (4.169)
By using (4.164), this implies that the streamfunction ψ * 12 behaves in the vicinity of the vortex core as

ψ * 12 = r2 2 -β 2 Ω c cos(2 θ) + σ 2 Ω 2 c sin(2 θ)t -χ 2 Ω 3 c cos(2 θ)t 2 dU dz 2 + O(r 4 ). (4.170)
In chapter 5, this will allow us to determine the evolution of the vertical shear of the horizontal velocity on the vortex axis. A similar calculation is conducted in appendix C to determine ψ ν 12 . For r 1, it is of the form:

ψ ν 12 = εr 2 2 ReF 2 h β ν 2 Ω c sin(2 θ) + χ ν 2 t cos(2 θ) + σ ν 2 Ω c t 2 sin(2 θ) dU dz 2 + O(

Evolution of the angular velocity of the vortex

By rescaling the slow time, the vertical coordinate, and the Reynolds number, the evolution equation (4.156) can be re-written after some manipulations:

∂ζ 0 ∂t = - F 2 h r ∂r 2 Ω 3 ∂r t dU dz 2 + t 2 2Re 1 r ∂ ∂r r ∂ζ 0 ∂r dU dz 2 . (4.173) x z z = l z /2 x = 0 h h x = 0 δh - Figure 4
.8: Lines of constant total buoyancy b t near the point x = 0, z = l z /2 at t = 0 (grey dashed lines) and t = t 0 (grey solid lines). The bold black solid line represents the vortex axis at t = t 0 . δh is the height variation at x = 0 of the iso-buoyancy line located at z = h + l z /2 initially.

3 Integrating in r leads to an equation for Ω, namely

∂Ω ∂t = -F 2 h tΩ 3 + t 2 2Rer ∂ζ 0 ∂r dU dz 2 .
(4.174)

This equation shows that the angular velocity of the vortex decays because of dynamic and viscous effects (first and second terms of the right-hand side, respectively). A physical interpretation of the origin of the dynamic effect can be gained by considering the conservation of potential vorticity near the vortex axis. First, let us consider the buoyancy at leading order b = √ εb 0 /F h = -txΩ 2 dU/dz, where b 0 is given in (4.144) and x = x-U (z)t. Figure 4.8 sketches the lines of constant total buoyancy b t = b+z/F 2 h near the point x = 0, z = l z /2 at two times: t = 0 (grey dashed lines) and t = t 0 (grey solid lines) where t 0 is small. At t = t 0 , the vortex axis x = 0 is inclined towards the left since x = U (z)t 0 -k z U S t 0 (z -l z /2) near z = l z /2. The iso-buoyancy lines are also no longer horizontal but slanted downward towards the positive x direction. On the vortex axis x = 0, they remain at the same vertical level while on the line x = 0, they are displaced by

δh = -F 2 h Ω 2 ht 2 0 dU dz 2 (4.175)
where h = z 0 -l z /2 with z 0 the initial vertical level of the iso-buoyancy line. δh is therefore negative as clearly seen in figure 4.8. Since the vortex is inclined, the 3 From this equation, it can be remarked that an alternative choice for the scaling of the slow time scale and Reynolds number could be τ

= √ εt, 1/ReF 2 h = √ ε/ Re b where Re b = O(1)
. Introducing these scalings in (4.173) leads to an equation where only the slow time scale τ appears instead of both τ and t as in (4.156). However, the asymptotic analysis with these scalings would give identical results as the present one. iso-buoyancy lines are indeed displaced towards the vortex axis where the pressure is minimum. This means that the vortex is squeezed near x = 0, z = l z /2. The conservation of potential vorticity between the iso-density lines at z = z 0 and z = l z /2 (figure 4.8)

ζ h = ζ + δζ h + δh , (4.176)
implies that the vertical vorticity at t = t 0 will decrease by δζ = ζδh/h. The variation of vertical vorticity between t = 0 and t = t 0 is therefore

∂ζ ∂t = δζ t 0 = -F 2 h ζΩ 2 t 0 dU dz 2 . (4.177)
This corresponds exactly to (4.174) for Re = ∞ since ζ = 2Ω on the vortex axis. However, we stress that (4.174) is valid not only on the vortex axis near z = l z /2 but everywhere. Both the dynamic and viscous decay of Ω are proportional to (dU/dz) 2 . Hence, the decay will be maximum at z = 0, l z /2 and will vanish at z = l z /4, 3l z /4 for the sinusoidal profile U = U S sin(k z z).

The equation (4.174) can be solved by a multiple scale analysis when the viscous effects are small. To do so, we introduce the time t = (F h t|dU/dz|) 3 . Then, (4.174) becomes

∂Ω ∂ t = - Ω 3 3 t1/3 + γ 6r ∂ζ 0 ∂r , γ = 1 ReF 3 h dU dz , (4.178)
where γ is the only remaining control parameter. Assuming γ 1, the angular velocity can be expanded as

Ω = Ω 0 + γΩ 1 + O(γ 2 ) (4.179)
and we introduce a slow time scale υ = γ t. At leading order, (4.178) becomes

∂Ω 0 ∂ t = - Ω 3 0 3 t1/3 (4.180)
whose solution is

Ω 0 = Ω i (r, υ) 1 + t2/3 Ω 2 i (r, υ) , (4.181)
where Ω i is the angular velocity profile at t = 0. It can depend also on the slow time υ. At next order, we obtain:

∂Ω 1 ∂ t + ∂Ω 0 ∂υ = - Ω 2 0 Ω 1 t1/3 + 1 6r ∂ζ 0 ∂r .
(4.182)

After substituting (4.181) in this equation, its general solution is found to be:

Ω 1 = Θ - ∂Ω i ∂υ + 1 6r ∂ζ i ∂r - 1 2 ∂Ω i ∂r 2 3 arctan(Ω i t1/3 ) Ω 4 i - 3 t1/3 Ω 3 i + t Ω i 1 + t2/3 Ω 2 i 3/2 , (4.183)
where Θ is an arbitrary constant. This constant should be set to zero to ensure that Ω 1 → 0 as r → +∞. In addition, we have to impose that Ω = Ω i initially, implying that Ω 1 = 0 at t = 0. This leads to an equation for Ω i over the slow time υ:

∂Ω i ∂υ = 1 6r ∂ζ i ∂r .
(4.184) Since Ω i (r, υ = 0) = [1 -exp(-r 2 )] /r 2 , the solution of (4.184) is:

Ω i = 1 -exp(-r 2 /ã 2 ) r2 , ( 4.185)
with ã2 = 1 + 2υ/3. Gathering the two orders and rescaling the slow time scale υ give the total expression of Ω:

Ω = Ω i 1 + t2/3 Ω 2 i - γ 2 ∂Ω i ∂r 2 3 arctan(Ω i t1/3 ) -3 t1/3 Ω i + tΩ 3 i Ω 4 i 1 + t2/3 Ω 2 i 3/2 + O(γ 2 ), (4.186)
with

Ω i = 1 r2 1 -exp - r2 1 + 2γ t/3 . (4.187)
Equation ( 4.186) implies that the angular velocity at the vortex center evolves as

Ω c = 1 1 + 2γ t/3 2 + t2/3 + O(γ 2 ), (4.188)
i.e. in terms of t and the original parameters:

Ω c = 1   1 + 2t 3 3Re dU dz 2   2 + F 2 h t 2 dU dz 2 .
(4.189)

Figure 4.9 compares the asymptotic solution (4.186) to the exact solution of (4.178) computed numerically at different times t. We see that the agreement is excellent for ). There are only some slight discrepancies in the vortex core r 1 when t 1. In practice, the asymptotic solution (4.186) turns out to be very accurate not only for small values of γ but also for large values such as γ = 10 (not shown). In figure 4.9, the Lamb-Oseen profile .190) with ã2 e = 1/Ω c , where Ω c is given by (4.188), has been also plotted with black dashdotted lines. This profile is close to the asymptotic solution (4.186). This is the reason why the approximate profile (4.190) has been used to ease the numerical resolution of the streamfunction ψ 12 (4.157).

Ω e = 1 -exp(-r 2 /ã 2 e ) r2 ( 4
Figure 4.10 shows the asymptotic angular velocity at the vortex center (4.188) and the exact solution obtained numerically, for different values of the control parameter γ. The agreement is very good not only for γ = 0.1 but also for γ = 1 and γ = 10. Therefore, (4.188) will be used to describe the evolution of Ω c regardless of the value of γ.

Analysis of the vortex deformations

Having determined the evolution of Ω, we now come back to the study of the equations (4.138)-(4.139) for the vortex displacements.

Since U = U S sin( kz z), we will look for the perturbation displacements in the form δx(z, η, τ ) = X(η, τ ) sin( kz z), δy(z, η, τ ) = Ŷ (η, τ ) sin( kz z), (4.191) where X and Ŷ are amplitudes that vary slowly along the vertical and with time.

Then, (4.138) and (4.139) become

∂ X ∂τ = - k2 z A(η, τ ) 2 Ŷ (η, τ ) + k2 z U S ϑ(η, τ ) + ∂C w ∂t (t, η, τ ) -ε k2 z U S t ReF 2 h , (4.192) ∂ Ŷ ∂τ = k2 z A(η, τ ) 2 U S t + X(η, τ ) + k2 z U S ∂S w ∂t (t, η, τ ). (4.193)
It is convenient to further rescale the slow time τ and the vertical coordinate η in terms of the time t and the vertical coordinate z: (4.195) where ω = (k z F h ) 2 A/2. The solution of (4.194)-( 4.195) can be decomposed as (4.197) giving (4.198)

∂ X ∂t = -ω(z, t) Ŷ (z, t) + (k z F h ) 2 U S ϑ(z, t) + ∂C w ∂t (z, t) - k 2 z U S t Re , (4.194) ∂ Ŷ ∂t = ω(z, t) U S t + X(z, t) + (k z F h ) 2 U S ∂S w ∂t (z, t),
X = X * + (k z F h ) 2 U S C w + t 0 ϑ(z, υ)dυ - k 2 z U S t 2 2Re , (4.196) Ŷ = Y * + (k z F h ) 2 U S S w + U S t 0 ω(z, υ)υdυ,
∂X * ∂t = -ωY * + O (k z F h ) 4 ,
∂Y * ∂t = ωX * + O (k z F h ) 4 , k 4 z F 2 h Re . ( 4
X * = X * c (z) cos t 0 ω(z, υ)dυ + X * s (z) sin t 0 ω(z, υ)dυ , (4.200) Y * = -X * s (z) cos t 0 ω(z, υ)dυ + X * c (z) sin t 0 ω(z, υ)dυ , (4.201)
where X * c and X * s are arbitrary functions of z only. The initial conditions X = Ŷ = 0 at (4.203) 1. the displacement due to internal waves (first terms in the right-hand side of (4.202)-(4.203))

t = 0 impose X * c = X * s = 0. Therefore, a solution of (4.194)-(4.195) valid up to order O [(k z F h ) 2 ] and O [k 2 z /Re] is X = (k z F h ) 2 U S C w + t 0 ϑ(z, υ)dυ - k 2 z U S t 2 2Re , (4.202) Ŷ = (k z F h ) 2 U S S w + U S t 0 ω(z, υ)υdυ,
2. the second term in the right-hand side of (4.202) describes a small delay (since ϑ is negative as seen in figure 4.5b) in the advection by the shear flow due to its three-dimensional character.

3. viscous effects (last term in the right-hand side of (4.202))

4. the self-induction (last term in the right-hand side of (4.203)). The shear flow bends the vortex in the x direction and, in turn, this creates a self-induced motion in the y direction.

Behaviour for small time

By substituting C w and S w in (4.194)-( 4.195) by their asymptotic expressions (B.21)-(B.22) derived for small time (see appendix B), we find that the amplitudes of the perturbation displacements behave for t 1 and in the inviscid limit as:

X = k 2 z U S ςt 3 + O(t 4 ), Ŷ = k 2 z U S σt 4 + O(t 5 ). (4.204)
Remarkably, the initial evolution of ( X, Ŷ ) is independent of the Froude number at leading order. At this order, X and Ŷ vary like t 3 and t 4 respectively. This implies that the leading order buoyancy b 0 (4.81) varies like t 3 for small times consistently with the short time expansion of b derived in §4.3. However, it is not possible to compare quantitatively the order k 2 z t 3 of the short time expansion to its counterpart in the longwavelength expansion. Indeed, (4.204) shows that X and Ŷ are of order k 2 z for small times. Thus, the contribution of the perturbation displacements δx and δy in (4.81) is one order smaller in k 2 z than the other terms. Hence, these contributions are actually of the same order as the next order term εb 1 in the expansion (4.51). Since b 1 will not be computed herein, the determination of the order k 2 z of b is therefore incomplete. It is interesting to compare the predictions (4.204) to those of a simple model where the advection by the shear flow is directly added to the equations of the self-induced motion of the single vortex, i.e.: [START_REF] Marshall | Tearing of an aligned vortex by a current difference in two-layer quasigeostrophic flow[END_REF] for two layer quasi-geostrophic flows, with constant coefficients. For small times t 1, the solution of (4.205)-(4.206) is:

∂∆ X ∂t = -ω Ŷ + U S , (4.205) ∂ Ŷ ∂t = ω∆ X, ( 4
X = - ω2 U S 6 t 3 + O(t 4 ), Ŷ = ωU S 2 t 2 + O(t 3 ). (4.207) Thus, X is initially of order k 4 z F 4 h t 3 instead of k 2 z t 3 and Ŷ is of order k 2 z F 2 h t 2 instead of k 2 z t 4
. This shows that internal waves have a strong effect on the initial dynamics of the vortex.

Evolution in the planes

z = l z /4, 3l z /4
In the planes z = l z /4, z = 3l z /4, A and ϑ are independent of t because the initial angular velocity profile does not evolve at these locations since dU/dz = 0 (see (4.174)). Thus, (4.202)-(4.203) become

X = (k z F h ) 2 U S [C w + ϑt] - k 2 z U S t 2 2Re , (4.208) Ŷ = U S ωt 2 2 + (k z F h ) 2 U S S
X = (k z F h ) 2 U S [C ∞ w + ϑt] - k 2 z U S t 2 2Re , (4.210) Ŷ = U S ωt 2 2 + (k z F h ) 2 U S S ∞ w .
(4.211)

Conclusion

We have carried out two asymptotic analyses of the evolution of a columnar vortex in an ambient stratified shear flow. The vortex has initially a Lamb-Oseen profile while the shear flow is horizontally uniform and varies sinusoidally along the vertical:

U (z) = U S sin(k z z).
The Froude number F h has been assumed to be lower than unity so that there is no critical layer where the angular velocity of the vortex is equal to the buoyancy frequency. The Reynolds number has been considered infinite in the first analysis and large in the second analysis. These asymptotic analyses highlight different aspects of the evolution of the vortex. First, the asymptotic analysis for small time shows that the horizontal flow is independent of the stratification up to order t 3 . At order t, the flow consists simply in a horizontal displacement of the vortex without deformation while at order t 2 , internal deformations of the vortex occur for finite vertical wavenumber k z . On the vertical, the dynamics is non-hydrostatic as long as t < 1/N .

The second asymptotic analysis is a long-wavelength analysis assuming that the rescaled vertical wavenumber k z F h is small. This particular small parameter derives from the self-similarity of strongly stratified flows [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] according to which the vertical wavenumber always appears multiplied by the Froude number F h in the Euler equations under the Boussinesq approximation.

This asymptotic analysis provides governing equations for the displacement of the vortex axis and for its angular velocity. The vortex axis is advected by the shear flow at leading order ∆x = U (z)t but four different effects influence its trajectory at order (k z F h ) 2 . First, the vortex axis deviates perpendicularly to the direction of the shear flow owing to the self-induction of the vortex when it is bent. Secondly, internal waves, which are generated at the start-up of the motion, influence the initial evolution of the vortex. These waves oscillate at the buoyancy pulsation N and explain why the dynamics is initially non-hydrostatic. Later, their effects on the vortex trajectory quickly decay off. The third and fourth effects are three-dimensional and viscous effects that slightly reduce the advection in the direction of the shear flow.

The angular velocity of the vortex decreases with time because of dynamic and viscous effects. The dynamic effect is due to the squeezing of the isopycnal surfaces when the vortex is inclined. To conserve potential vorticity, the core vertical vorticity has therefore to decrease. The viscous decay of the angular velocity is also quickly enhanced since the vertical shear increases algebraically with time. These two effects are not uniform along the vertical: they are maximum in the regions of high shear and vanish at the vertical levels where the shear is zero.

In addition, a horizontal velocity field with an azimuthal wavenumber m = 2 arises at order (k z F h ) 2 implying that the vortex becomes slightly elliptical.

Chapter 5

Evolution of a vortex in a strongly stratified shear flow. Part 2. Numerical simulations.

Article in preparation for submission to the Journal of Fluid Mechanics.

85

Introduction

In this chapter, we continue the analysis of the evolution of a vortex embedded in a vertically sheared flow in a strongly stratified fluid. The main purpose is to determine if the vertical shear can grow indefinitely as speculated by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] and thus leads to the unconditional development of the shear instability.

This instability is thought to be an important process for the generation of small scales in stratified flows [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]. In the case of a columnar counter-rotating vortex pair, [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF] and [START_REF] Waite | Instability and breakdown of a vertical vortex pair in a strongly stratified fluid[END_REF] have reported that the vertical shear generated by the zigzag instability can lead to the development of the shear instability. This occurs when the buoyancy Reynolds number Re b = ReF 2 h (Re is the classical Reynolds number and F h the Froude number) is above a threshold since the minimum Richardson number is inversely proportional to Re b [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF][START_REF] Augier | Onset of secondary instabilities on the zigzag instability in stratified fluids[END_REF]. The subsequent destabilization of the Kelvin-Helmholtz billows leads to small-scale turbulence with spectral characteristics similar to those of randomly forced stratified turbulence [START_REF] Augier | Spectral analysis of the transition to turbulence from a dipole in stratified fluid[END_REF][START_REF] Waite | The vortex instability pathway in stratified turbulence[END_REF].

However, a counter-rotating vortex pair is a particular flow not frequently encountered in turbulent flows. Here, we consider the more generic configuration of a single vortex in an ambient shear flow. In chapter 4, we have studied such flow by means of two asymptotic analyses. First, a short-time analysis has evidenced the existence of an initial non-hydrostatic regime. Secondly, a long-wavelength analysis has been carried out for k z F h 1, where k z is the vertical wavenumber of the sinusoidal shear flow. This analysis provides a complete description of the vortex dynamics: the evolution of the vortex axis and angular velocity as well as secondary flows created as the vortex is bent. In the present chapter, we will conduct DNS of this flow and analyse its dynamics in the light of the asymptotic analyses.

The chapter is organized as follows. The initial conditions, control parameters, and numerical method are described in §5.2. An overview of two typical simulations is first given in §5.3. Then, the initial evolution of the flow is compared to the results of the short-time asymptotic analysis in §5.4. In §5.5, the long-wavelength analysis is first briefly summarized in §5.5.1 and its predictions for the deformations of the vortex axis are compared to the numerical simulations in §5.5.2. We then focus on the evolution of the flow at the vortex center and the mid-vertical level where the vertical shear is maximum ( §5.5.3). Again, the asymptotic analysis is used to rationalize the numerical results. Finally, section §5.5.4 concentrates on the Richardson number predicted by the long-wavelength analysis. Section §5.7 summarizes and discusses the results.

Formulation of the problem

Initial conditions and governing equations

As in chapter 4, the flow at t = 0 is chosen as

u = U S + u v , (5.1)
where U S is a sinusoidal shear flow and u v a columnar vortex with a Lamb-Oseen profile: (5.2) where (x, y, z) and (r, θ, z) are cartesian and cylindrical coordinates, respectively. The associated unit vectors are (e x , e y , e z ) and (e r , e θ , e z ). In (5.2), the length and time have been non-dimensionalized by the vortex radius a 0 and the turnover time 2πa 2 0 /|Γ| of the vortex. The shear amplitude U S and wavenumber k z are therefore non-dimensional: U S = ÛS /(|Γ|/2πa 0 ), k z = kz a 0 , where ÛS and kz are the corresponding dimensional quantities.

U S = U S sin(k z z)e x , u v = 1 -exp(-r 2 ) r e θ ,
The governing equations are the incompressible Navier-Stokes equations under the Boussinesq approximation. They read in dimensionless form: 6000 3π/2 0.2 1.333 832 832 448 0.005 Table 5.1: Overview of the physical and numerical parameters of some typical simulations. For all simulations, the horizontal dimensions of the domain are l x = l y = 18.

∇ h • u h + ∂w ∂z = 0, (5.3) ∂u h ∂t + (u h • ∇ h ) u h + w ∂u h ∂z = -∇ h p + 1 Re ∇ 2 u h , (5.4) ∂w ∂t + u h • ∇ h w + w ∂w ∂z = - ∂p ∂z + b + 1 Re ∇ 2 w, (5.5) ∂b ∂t + u h • ∇ h b + w ∂b ∂z = - w F 2 h + 1 ReSc ∇ 2 b, (5.6) F h Re k z U S l z n x n y n z δt 0.
where u h = (u, v) and w are the horizontal and vertical velocities in cartesian coordinates, p the pressure, b the buoyancy, and

Re = |Γ| 2πν , F h = |Γ| 2πa 2 0 N , Sc = ν κ (5.7)
are the Reynolds, Froude and Schmidt numbers, with ν the viscosity, κ the diffusivity and N the Brunt-Väisälä frequency which is assumed constant.

Numerical method

Equations (5.3-5.6) are integrated numerically by means of a pseudo-spectral method with periodic boundary conditions and a fourth-order Runge-Kutta time advancement scheme [START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]. An elliptic truncation of the top one-third of the modes is applied. The viscous and diffusive terms are integrated exactly.

The horizontal size of the computational domain is taken large l x = l y = 18 in order to minimize the effect of the periodic boundary conditions. The vertical size is set to l z = 2π/k z , so that a single wavelength of the shear flow is simulated.

Table 5.1 lists the parameters of some typical simulations. The number of grid points in the x and y directions have been varied from n x = n y = 384 to n x = n y = 832 depending on the values of the Reynolds and Froude numbers. The number of grid points in the vertical direction ranges from n z = 256 to n z = 448 depending on the values of k z , Re and F h . Typically, a high resolution is required for the parameters where the shear instability develops because it generates small billows while a moderate resolution is sufficient for the other cases. When k z increases, the horizontal and vertical resolutions have to be increased also since the vertical gradients are larger. The accuracy of the results has been checked by increasing the resolution or the domain horizontal sizes in several runs. The time step varies from δt = 0.0025 to δt = 0.01. All the numerical simulations have been carried out for Sc = 1. The Froude number has been always kept below unity so as to remain in the strongly stratified regime. The shear amplitude U S is also always kept below unity meaning that the vortex is stronger than the shear flow. The vertical wavenumber has been varied in the range 0.1 k z 3π/2. The Reynolds number has been varied from Re = 3000 to Re = 10000.

Overview of the dynamics

Qualitative description

We first begin by a description of two different simulations in order to give a global overview of the flow dynamics.

Figures 5.1 and 5.2 display the potential vorticity at different times for k z = π, F h = 0.1 and k z = 3π/2, F h = 0.5, respectively, whereas U S and Re are fixed to U S = 0.2 and Re = 6000. The first column shows three-dimensional contours while the second column represents a corresponding horizontal cross-section at the vertical level z = l z /4. The vortex is mostly displaced in the direction of the shear flow, but also slightly in the perpendicular direction as seen in the horizontal cross-sections. Hence, the vertical plane containing the vortex axis is actually oblique relative to the (x, z) plane. The displacement in the y direction is weaker in figure 5.1 than in figure 5.2.

A common feature of both simulations is that the potential vorticity decreases faster in the regions of high shear z = 0, l z /2 than in the regions of weak shear z = l z /4, 3l z /4. Thus, the vortex seems to be torn apart into two separate pancake vortices [START_REF] Beckers | Dynamics of pancake-like vortices in a stratified fluid: experiments, model and numerical simulations[END_REF] at large times.

Figure 5.3 displays the corresponding total vertical shear of the horizontal velocity √ S z = (∂u/∂z) 2 + (∂v/∂z) 2 (color) in the vertical cross-section at y = 9, i.e. passing through the vortex center at t = 0. The superimposed black lines show the total density ρ t . For k z = π, F h = 0.1 (left column of figure 5.3), the shear is maximum in the vortex core at the point x c = 9, z c = 0, l z /2 (note that these coordinates correspond to those of the computational domain where the vortex center is initially in the middle x = 9, y = 9). As the vortex is progressively bent, S z grows monotonically with time and becomes rapidly much higher than the maximum ambient shear max ( Sz ) = k 2 z U 2 S (figure 5.3a,c,e,g). The iso-density lines remain nearly flat since the stratification is strong for this case. Figure 5.4 shows that the minimum of the Richardson number (black solid line)

Ri = 1 F 2 h + ∂b ∂z S z
(5.8) decreases with time down to min (Ri) = 3.7 at t = 22 and then slowly re-increases.

The quantity min (Ri) thus remains well above the critical value 1/4 necessary for the development of the shear instability of a steady parallel inviscid shear flow [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Note on a paper of John W. Miles[END_REF].

For k z = 3π/2, F h = 0.5 (right column of figure 5.3), the growth of the maximum shear S z is not monotonic. There is a first stage where the shear is very weak within the vortex core (see figure 5.3b at t = 4), i.e. the response of the vortex tends to cancel the ambient shear. Then, the vortex becomes tilted as for k z = π, F h = 0.1, and the maximum shear is encountered in the vicinity of x c = 9, z c = 0, l z /2 (figure 5.3d). The regions of high shear are remarkably thin. Later on, the flow strongly dissipates in these regions and the shear becomes maximum at points away from (x c , z c ) (figure 5.3f,h). During this evolution, the iso-density lines are strongly deformed in contrast to k z = π, F h = 0.1. Some overturns can even be seen in the regions of high shear at t = 26 (figure 5.3f ). As seen in figure 5.4, the minimum Richardson number for this simulation (grey solid line) decreases below Ri c = 0.25 for 11 t 37. The (y, z) cross-section at x = 9 of the buoyancy at t = 26 (figure 5.5b) confirms the presence of Kelvin-Helmholtz billows near z = 0, l z /2. In contrast, no billows can be seen in these regions in the corresponding (x, z) cross-section at y = 9 (figure 5.5a). This means that the axes of the Kelvin-Helmholtz billows are mostly oriented in the x direction, i.e. they are parallel to the direction of the ambient shear flow. The black contours in figure 5.5 delineate the regions where Ri < 0.25. In addition to the unstable regions near (x c , z c ), there exist also other unstable regions above and below each pancake vortex at z = l z /4 and z = 3l z /4. We can also see some billows in these regions (figures 5.5a and 5.3f ) but in this case, their axes are perpendicular to the direction of the ambient shear. These unstable regions appear only in a second stage after those near (x c , z c ).

Time evolution of the vertical shear at the center

For k z = π, F h = 0.1, the Richardson number is always minimum at the center point (x c = 9, y c = 9, z c = l z /2) and at the symmetric point (x c = 9, y c = 9, z c = 0). For k z = 3π/2, F h = 0.5, the Kelvin-Helmholtz instability also develops first at these points. It is therefore interesting to investigate the evolution of the vertical shear at these locations.

To this end, we first decompose the flow into a mean flow varying only along the vertical and with time ū(z, t) and a complementary flow u * :

u = ū + u * ,
(5.9)

where the overbar denotes the horizontal average over the computational domain, which for any quantity q is defined as q = 1 l x l y ly 0 lx 0 q(x, y, z, t)dxdy.

(5.10) At t = 0, we have ū = U S and u * = u v so that ū and u * will be called "shear flow" and "vortex flow", respectively. Figure 5.6a shows the evolution of the different shear components ∂ ū/∂z, ∂u * /∂z, and ∂v * /∂z at the center (x c = 9, y c = 9, z c = l z /2) for F h = 0.1, k z = π, U S = 0.2 and Re = 6000. The quantity ∂v/∂z is always equal to zero at the center and is not plotted. More generally, ∂v/∂z always remains very small at any vertical position compared to ∂ ū/∂z. We see that the mean shear ∂ ū/∂z (dash-dotted line) remains almost constant and equal to -k z U S = -0.6283.

In contrast, the shear component ∂v * /∂z (solid line) grows first linearly and then saturates at t 22 at the value ∂v * /∂z = 6.3, i.e. ten times the maximum ambient shear |k z U S |. The other component ∂u * /∂z remains very weak up to t = 10 and then increases up to ∂u * /∂z 1 at t = 30. This quantity therefore saturates at a lower value and later than its counterpart ∂v * /∂z.

The initial behaviour of the vertical shear of the vortex ∂u * /∂z and ∂v * /∂z can be simply understood by considering that the vortex is displaced at the velocity U (z) in the x direction, i.e. u v (x -U t, y), giving:

∂u v ∂z = - dU dz t ∂u v ∂x .
(5.11)

Since u v = (-Ωy, Ωx), where Ω(r) is the angular velocity of the vortex, we have ∂v v /∂x = Ω = 1 and ∂u v /∂x = 0 at the center r = 0. Thus, (5.11) yields

∂u v ∂z = 0, ∂v v ∂z = k z U S t
(5.12) at z = l z /2. The straight line in figure 5.6a confirms that ∂v * /∂z increases initially at the rate k z U S t. This also explains why ∂u * /∂z remains very small initially. The subsequent evolutions will be explained later thanks to the asymptotic analysis performed in chapter 4. In contrast to the case F h = 0.1, k z = π (figure 5.6a), ∂v * /∂z follows the law (5.12) only at the very beginning t 2. Instead, both shear components ∂u * /∂z and ∂v * /∂z first oscillate with a phase lag and with a period around 2π, i.e. the period corresponding to the angular velocity on the vortex axis Ω = 1. Because of these oscillations, we can notice that ∂v * /∂z goes back to zero around t 4 -5 while ∂u * /∂z is approximately opposite to ∂ ū/∂z. Thus, the total shear S z is weak at the center as already observed in figure 5.3b at t = 4.

Then, ∂u * /∂z and ∂v * /∂z are both abruptly amplified up to an absolute value around 10. Remarkably, ∂u * /∂z becomes now negative and saturates earlier than ∂v * /∂z. Later on, |∂u * /∂z| decreases very quickly while |∂v * /∂z| decays more slowly. The vorticity components -ω *

x and ω * y have been also plotted in figure 5.6b. They are again almost identical to ∂v * /∂z and ∂u * /∂z except -ω *

x for 21 t 35. This corresponds to the time interval when the Kelvin-Helmholtz billows exist. They produce a finite vertical velocity w * , making the term ∂w * /∂y in ω *

x no longer negligible. In contrast, the term ∂w * /∂x is still negligible in ω * y , most probably because the axes of the Kelvin-Helmholtz billows are aligned with the x direction. previously described. The kinetic energies have been decomposed into a mean part and a vortex part using the decomposition (5.9):

Global energy and enstrophy evolutions

Ēk = 1 l z V ū2 2 dV , E * kh = 1 l z V u * h 2 2 dV , E * kz = 1 l z V w * 2 2 dV .
(5.13)

The integral over the computational domain V is divided by l z in order to enable the comparisons between simulations carried out with distinct vertical wavelengths.

Similarly, the global potential energy per vertical length unit is:

E p = 1 l z V F 2 h b 2 2 dV .
(5.14)

The kinetic energy of the mean flow Ēk (black dash-dotted lines) remains approximately constant even if it increases slightly at the beginning mostly for k z = 3π/2, F h = 0.5 (figure 5.7b). In contrast, the horizontal kinetic energy of the vortex E * kh (black dashed lines) decreases regularly following an approximately linear trend. The vertical kinetic energy E * kz (grey dashed lines) and the potential energy E p (grey solid lines) remain always very weak compared to the horizontal kinetic energy.

Likewise, figure 5.8 displays the evolutions of the global enstrophies per vertical length unit, decomposed using (5.9):

Z = 1 l z V ω2 2 dV , Z * h = 1 l z V ω * h 2 2 dV , Z * z = 1 l z V ζ * 2 2 dV .
(5.15)

The global horizontal enstrophy of the vortex Z * h (black dashed lines) increases and then decreases while its vertical counterpart Z * z (grey dashed lines) continuously decays. Unexpectedly, the growth of Z * h is much more pronounced for k z = π, F h = 0.1 than for k z = 3π/2, F h = 0.5 although the vertical wavenumber k z is higher in this second case. The enstrophy of the mean shear Z remains approximately constant like the The contour interval is 0.002 for w and 0.001 for b. The contour at the level zero is not drawn. The negative iso-contours are colored in grey. The parameters are F h = 0.5, k z = π, U S = 0.2, and Re = 6000.

mean kinetic energy Ēk . Since the mean enstrophy Z for k z = 3π/2 is more than twice the one for k z = π, the maximum of the total enstrophy (solid lines) is comparable in the two simulations even if max (Z * h ) is lower for k z = 3π/2.

Comparison to the short-time asymptotic analysis

After this global description of two simulations, we turn to a detailed comparison of the initial dynamics with the predictions of the asymptotic analysis for small time performed in chapter 4 in the inviscid limit. The main result of this analysis is that the buoyancy and vertical velocity should behave initially as

b = t 3 k z U S 6F 2 h β(r) cos(θ) cos(k z z), (5.16) w = -t 2 k z U S 2 β(r) cos(θ) cos(k z z), (5.17)
where β is a function which depends on r and k z (see (4.34)-(4.35) and figure 4.2 in chapter 4). The expressions (5.16)-(5.17) imply that the dynamics is non-hydrostatic, i.e. b < w, as long as t < F h . These asymptotic predictions will be first checked against a simulation for the parameters F h = 0.5, k z = π, U S = 0.2, Re = 6000, and then a parametric study will be carried out.

A detailed comparison

Figures 5.9 and 5.10 compare the asymptotic predictions for w and b to the observed fields in the DNS at z = l z /2, at two times t = 0.4 and t = 1.6, respectively. As seen in (5.16)-(5.17), z = l z /2 is the vertical level where w and b are maximum. There is a good quantitative agreement between asymptotics and numerics at t = 0.4 both for the vertical velocity (figure 5.9a,b) and the buoyancy (figure 5.9c,d). We can notice that the buoyancy is in phase opposition with the vertical velocity. As expected, the asymptotics depart from the numerical results for t > 1 as can be seen in figure 5.10 at t = 1.6. Indeed, the wavenumber one pattern of the buoyancy has started to rotate in the numerical simulation (figure 5.10d), while such rotation is not reproduced by the asymptotic solution (figure 5.10c). Moreover, the extremum values predicted by the asymptotics no longer perfectly match those in the DNS. Nevertheless, the structure of the vertical velocity (figure 5.10a,b) is qualitatively similar in the asymptotics and numerics.

In order to estimate the time range of validity of the asymptotics, the evolutions of the peak values w m and F 2 h b m of the vertical velocity and buoyancy are plotted in figure 5.11. An excellent agreement between theory (black solid lines) and numerics for F h = 0.5 (black dashed lines) is observed for t 0.8 approximately. Beyond this time, the growth of the maximum vertical velocity and buoyancy is slowlier in the numerics than predicted by the asymptotics.

Effect of the Froude number

Figure 5.11 also shows w m and F 2 h b m from a DNS for the larger Froude number F h = 0.8 (grey dashed lines). In agreement with the asymptotics, w m and F 2 h b m are independent of the Froude number for t < 1 whereas they become dependent for larger time.

According to (5.16)-(5.17), the buoyancy becomes larger than the vertical acceleration when t > √ 6F h , i.e. sooner as the Froude number diminishes. In other words, the duration of the initial non-hydrostatic regime scales with the Froude number.

Effect of the vertical wavenumber

The effect of the vertical wavenumber k z on the maximum vertical velocity w m and buoyancy b m is investigated in figure 5.12. Both the asymptotic and the numerical results are almost distributed on a single curve. This is due to the fact that the factor k z β appearing in (5.16)-(5.17) is approximately independent of k z in the range π/2 k z 3π/2 as illustrated in figure 4.2 of chapter 4. In contrast, for small k z , w m and b m are proportional to k z since β is almost independent of the vertical wavenumber in that limit (not shown).

Effect of the shear flow amplitude U S

The influence of the shear amplitude U S is considered in figure 5.13. The maximum vertical velocity and buoyancy are both proportional to U S in agreement with the asymptotic predictions. Even if there is a departure between the asymptotics and the numerics, it is remarkable that the scaling relations w ∝ U S and b ∝ U S hold very well in the DNS over the entire time range investigated.

Comparison to the long-wavelength asymptotic analysis

Reminder

In chapter 4, we have performed a long-wavelength asymptotic analysis for k z F h 1, i.e. for small vertical Froude number

F v = k z F h = |Γ|/(a 0 l z N ).
Leading order viscous effects have been also taken into account. This analysis has provided evolution equations for the position of the vortex center at each level z:

∆x = U (z)t + δx(z, t),
(5.18) ∆y = δy(z, t),

(5.19) where

δx = X(z, t) sin(k z z), (5.20) δy = Ŷ (z, t) sin(k z z).
(5.21)

The amplitudes ( X, Ŷ ) follow the equations (5.23) where C w and S w are the effects of internal waves excited at t = 0 and that decay quickly afterwards. The parameters ω and ϑ correspond to the self-induction of the vortex and an advection correction, respectively. The last term in the right-hand side of (5.22) is the leading viscous effect. The expressions of all the parameters in (5.22)-(5.23) are given in chapter 4.

X = (k z F h ) 2 U S C w + t 0 ϑ(z, υ)dυ - k 2 z U S t 2 2Re , (5.22) Ŷ = (k z F h ) 2 U S S w + U S t 0 ω(z, υ)υdυ,
In addition, the asymptotic analysis has shown that the vortex remains axisymmetric at leading order but its angular velocity evolves according to

∂Ω ∂t = -F 2 h tΩ 3 + t 2 2Rer ∂ζ 0 ∂r dU dz 2 ,
(5.24)

where ζ 0 = (1/r)∂r 2 Ω/∂r is the vertical vorticity and r is the local radius with respect to the center of the vortex at the level z:

r2 = (x -U (z)t -δx) 2 + (y -δy) 2 .
(5.25)

The associated local angle is denoted θ. The first term in the right-hand side of (5.24) ensures the conservation of potential vorticity while the second term describes the leading viscous effect. This effect is proportional to t 2 because the vertical shear grows like tdU/dz at leading order. In chapter 4, the equation (5.24) has been solved asymptotically when

γ = 1 ReF 3 h dU dz
(5.26) is small. In particular, an analytic expression for the angular velocity on the vortex axis has been obtained:

Ω c = 1 (1 + 2γT 3 /3) 2 + T 2 , ( 5.27)
where T = F h t|dU/dz|. This expression is in very good agreement with the exact result obtained by numerical integration of (5.24). It will be therefore used in the following.

A last prediction of the asymptotic analysis is that the horizontal velocity of the vortex at order O[(k z F h ) 2 ] is of the form

u h1 = ∇ h Φ 1 -∇ × ψ 1 e z ,
(5.28) with Φ 1 = Φ 10 + Φ 11 + Φ 12 , ψ 1 = ψ 11 + ψ 12 .

(5.29)

The potential Φ 10 does not depend on θ and thus corresponds to an axisymmetric velocity field in the local reference frame with origin at the vortex center. In contrast, (Φ 11 , ψ 11 ) and (Φ 12 , ψ 12 ) are contributions with an azimuthal wavenumber m = 1 or m = 2, respectively. They correspond to a displacement and an elliptic deformation of the vortex, respectively. All these functions have been determined in chapter 4. These results will enable us to obtain predictions for the evolution of the vertical shear of the horizontal velocity in section §5.5.3. Before, we begin by presenting in the next section a comparison between the asymptotic predictions for the vortex center (5.18)-(5.19) and the DNS.

Deformations of the vortex axis

In order to estimate the position of the vortex center in the numerical simulations, we have used two methods: one based on the potential vorticity Π = ω • [∇b + N 2 e z ] and the other on the vertical vorticity ζ. In each case, the displacements of the vortex center have been estimated from vorticity centroids: (5.31) where the brackets denote ϕ h = ϕ>ϕc ϕdxdy.

∆x Π c (z, t) = xΠ h Π h , ∆y Π c (z, t) = yΠ h Π h , (5.30) or ∆x ζ c (z, t) = xζ h ζ h , ∆y ζ c (z, t) = yζ h ζ h ,
(5.32)

The horizontal integration is carried out only in the regions where the vorticity is larger than a critical value Π c , ζ c in order to avoid taking into account the small background vorticity due to the fact that the total vorticity is zero owing to the use of periodic boundary conditions. The values Π c = 0.05 max t=0 (Π) and ζ c = 0.05 max t=0 (ζ) have been chosen as they provide results almost independent of the size of the computational domain and the particular values of the thresholds.

The tracking method based on the potential vorticity seems more natural since it is a transported quantity in the inviscid limit. However, we shall see that the method based on the vertical vorticity will enable a closer comparison to the asymptotic results. This is because the condition used to normalize the streamfunction at first order implies (5.33) and, similarly ∆y ζ c = ∆y, where (∆x, ∆y) are the asymptotic displacements (5.18)-(5.19).

∆x ζ c = xζ h ζ h = xζ 0 (r) h ζ 0 h = ∆x,
In the next sections, we compare the asymptotic and numerical results for different parameters.

In-depth analysis of a simulation Figure 5.14 compares the total displacements ∆x = U S t + δx and δy as predicted by the asymptotics to the positions ∆x c and ∆y c of the vortex estimated from the vertical vorticity and potential vorticity centroids at the level z = l z /4 for F h = 0.5, k z = 0.3, U S = 0.2, and Re = 6000. The agreement is excellent even at t = 20. The dominant displacement is in the x direction (figure 5.14a) and given by U S t while the deviations δx and δy are much smaller.

Figures 5.14b,c display the y displacement of the vortex estimated from the vertical vorticity, ∆y ζ c , and from the potential vorticity, ∆y Π c , respectively. They are both in very good agreement with the asymptotic prediction δy even if a small discrepancy arises after t 12. However, if we focus on the initial evolution of ∆y ζ c and ∆y Π c (figure 5.14d), we see that they are actually different for t 4. The displacement ∆y ζ c (crosses) is first slightly negative for t 2 in excellent agreement with the asymptotic prediction δy (black solid line). In contrast, the y displacement estimated from the potential vorticity ∆y Π c (open circles) increases monotonically. The estimation of the vortex center from the vertical vorticity is in much better agreement with the asymptotics than the estimation from the potential vorticity because the normalisation condition used in the asymptotic analysis is based on the vertical vorticity.

Nevertheless, if the effects of the internal waves are neglected (i.e. C w = S w = 0), the asymptotic prediction (black dashed line) is then close to ∆y Π c . This confirms that internal waves play a key role at the start-up of the flow evolution. Because of these waves, the initial evolution of δy is of the form δy = k 2 z U S σt 4 , where σ = -3.826 × 10 -3 is a constant, as shown by the grey solid line in figure 5.14d. In contrast, when internal waves are neglected, δy evolves initially as δy ∝ t 2 (black dashed line).

We have also compared the vertical velocity and buoyancy fields predicted by the asymptotics against their numerical counterparts. Figures 5.15 and 5.16 display horizontal cross-sections of w and b, respectively, in the plane z = l z /2 where they are maximum. A very good qualitative and quantitative agreement is observed even at t = 18 (figures 5.15c,f and 5.16c,f ) apart from the existence of small wave-like disturbances in the DNS that are absent in the asymptotics. We can remark that the vertical velocity and buoyancy are approximately in phase opposition at t = 2 (figures 5.15a,d and 5.16a,d) as observed for short times (figures 5.9, 5.10) while for t = 6 and t = 18, they have a π/2 phase difference. The close agreement between the asymptotic and numerical results can be further seen in the temporal variations of the maximum values of w and b (figure 5.17). The asymptotic and numerical results begin to slightly depart from each other only at t 10. The slowly growing discrepancy as time increases is probably due to the fact that the vertical derivative grows like k z U S t at leading order owing to the shear flow. Hence, the long wavelength assumption F h ∂/∂z 1 is expected to be less and less valid as time increases. The maximum values of w and b increase approximately linearly with oscillations superimposed. These oscillations are due to the internal waves excited at t = 0. Two periods T = 2πF h = π and T = Ω(r = 0) = 1 are mixed explaining why the oscillations look somewhat irregular, especially for w m (figure 5.17a).

Parametric study

We now examine the effects of the parameters F h , k z , and U S .

Figure 5.18a presents the evolution of the y displacement for a different Froude number and wavenumber F h = 0.1, k z = 1.5 such that the product k z F h = 0.15 is the same as in figure 5.14. The shear flow amplitude and Reynolds number have been kept to the same values U S = 0.2 and Re = 6000. As expected, the evolution of δy (black solid line) is similar to the one for k z = 0.3, F h = 0.5 (shown again by a black dashed line in figure 5.18a for reference). This confirms that the relevant main control parameter is k z F h and not k z and F h , separately. There is a small difference between the curves for F h = 0.1, k z = 1.5 and F h = 0.5, k z = 0.3 because the functions ω, ϑ, C w and S w depend on the Froude number. However, this dependence is weak when F h 0.5.

In contrast, if only F h (figure 5.18b) or only k z (figure 5.18c) are varied relative to the reference case, the y displacement differs since k z F h has changed. We can notice that the agreement of the asymptotic predictions δy with the measured y displacement of the vertical vorticity centroids is slightly better than for the reference case since k z F h is lower. Finally, lowering the shear flow amplitude to U S = 0.1 while keeping the other parameters constant (figure 5.18d) reduces the y displacement as expected. The agreement between the theoretical prediction and the numerical results is as good as for the reference case (figure 5.14). Figure 5.19 shows the evolution of the maximum buoyancy for the same four set of parameters as in figure 5.18. Once again, the asymptotics closely follow the numerics for all times investigated in the four cases all the more than k z F h is low. We can notice that the oscillations are strongly reduced when F h = 0.1 (figure 5.19a,b) compared to F h = 0.5 (figure 5.19c,d). This is because the amplitude of the internal waves in the buoyancy scales as the Froude number. Hence, they vanish in the hydrostatic limit F h = 0.

Evolution of the vertical shear

We now investigate the evolution of the flow at the vortex center r = 0 and at the levels z c = 0, l z /2 where the ambient shear is maximum. As seen in section §5.3, these are the locations where the Richardson number reaches its minimum or where the Kelvin-Helmholtz instability first appears.

Using the asymptotic expressions of the velocity and buoyancy fields, it is possible to obtain simple predictions for the vertical shear of the horizontal velocity and the buoyancy at these locations. This should enable us to predict the evolution of the Richardson number and the occurrence of the Kelvin-Helmholtz instability.

The horizontal velocity of the vortex has been written in the asymptotic analysis in the general form u h ≡ u h (x, ỹ, z, t), (5.34)

where we recall that (x, ỹ) are the local cartesian horizontal coordinates with origin at the vortex center at the level z. Hence, differentiating with respect to z leads to

∂u h ∂z = - ∂∆x ∂z ∂u h ∂ x - ∂∆y ∂z ∂u h ∂ ỹ + ∂u h ∂z x,ỹ .
(5.35)

At the vortex center x = ỹ = 0 and the level z = l z /2, we have (∂u h /∂z) x,ỹ = 0.

In addition, we can assume ∂∆x/∂z tdU/dz, ∂∆y/∂z 0 since the displacement corrections (δx, δy) are small compared to U (z)t as seen in the previous section. Hence, (5.35) reduces to

∂u h ∂z (r = 0, z = l z /2) = -t dU dz ∂u h ∂ x = k z U S t ∂u h ∂ x .
(5.36) Therefore, only the x derivatives of the horizontal velocity field are needed to compute the vertical shear of the vortex at r = 0, z = l z /2. In the asymptotic analysis, its horizontal velocity has been written in the form

u h = u h0 + (k z F h ) 2 u h1 + O (k z F h ) 4 , (5.37) with u 0 = -Ω(r, z, t)ỹ, v 0 = Ω(r, z, t)x (5.38)
while u h1 is defined in (5.28). This leads to

∂u ∂ x (r = 0, z = l z /2) = (k z F h U S ) 2 3 4 -σ 2 tΩ 2 c - 1 ReF 2 h β ν 2 Ω c + σ ν 2 t 2 Ω c , (5.39) ∂v ∂ x(r = 0, z = l z /2) = Ω c -(k z F h U S ) 2 β 2 Ω c + 1 4 + χ 2 t 2 Ω 3 c - χ ν 2 t ReF 2 h ,
(5.40)

where Ω c (t) is the angular velocity at the vortex center and σ 2 , β 2 , χ 2 , β ν 2 , σ ν 2 and χ ν 2 are constants given in chapter 4. Using (5.36), we therefore obtain

∂u ∂z (r = 0, z = l z /2) = (k z U S ) 3 F 2 h t 3 4 -σ 2 tΩ 2 c - 1 ReF 2 h β ν 2 Ω c + σ ν 2 t 2 Ω c , (5.41) ∂v ∂z (r = 0, z = l z /2) = k z U S tΩ c -(k z U S ) 3 F 2 h t β 2 Ω c + 1 4 + χ 2 t 2 Ω 3 c - χ ν 2 t ReF 2 h .
(5.42) Similarly, the vertical gradient of the buoyancy can be obtained as

∂b ∂z (r = 0, z = l z /2) = - ∂∆x ∂z ∂b ∂ x (5.43)
by neglecting (δx, δy) and because (∂b/∂z) x,ỹ = 0 at r = 0. Using the leading order expression of the buoyancy in the hydrostatic limit, we have

∂b ∂ x(r = 0, z = l z /2) = k z U S tΩ 2 c (5.44) giving therefore ∂b ∂z = (k z U S tΩ c ) 2 .
(5.45)

We now compare these predictions to the results of the DNS presented in section §5.3 for F h = 0.1, k z = π, U S = 0.2, and Re = 6000.

Detailed comparisons for a reference case

Figure 5.20 compares ∂u/∂x, ∂v/∂x, and ∂b/∂x at the vortex center and at z c = l z /2 predicted by (5.39), (5.40), and (5.44), respectively, to their counterparts from the DNS. Apart from ∂u/∂x, the agreement is excellent over the entire time range investigated.

In particular, the saturation of the buoyancy gradient is properly accounted for (figure 5.20c). The beginning of the evolution of ∂u/∂x is well predicted by the asymptotics until t 5 but not later. This discrepancy might be due to the fact that the asymptotic model is not uniformly asymptotic in time, since the long-wavelength assumption is expected to be no longer valid when k z F h U S t > 1. Figure 5.21 compares ∂u/∂z, ∂v/∂z and ∂b/∂z at the vortex center and at z c = l z /2 from the DNS and predicted by (5.41), (5.42) and (5.45). Except for ∂u/∂z once again, a very good agreement is observed at all times. In particular, the times and levels of saturations of ∂v/∂z and ∂b/∂z are both accurately estimated. Note that ∂v/∂z and ∂b/∂z saturate simultaneously because these two quantities are functions of a unique arrangement, namely tΩ c , except the small viscous term in (5.42). As already seen in figure 5.6a, the level of saturation of ∂u/∂z is much lower than for ∂v/∂z. From (5.41)-(5.42), we can see that this is because ∂u/∂z is one order smaller in (k z F h ) 2 than ∂v/∂z. This comes from the fact that the leading order velocity is such that ∂u 0 /∂ x = 0 while ∂v 0 /∂ x = Ω c .

In figure 5.21, we have also plotted with black dashed lines the evolution of each asymptotic quantity in the inviscid limit, i.e. for Re = ∞. In this limit, the angular velocity at the vortex center follows the law

Ω c = 1 1 + (k z F h U S t) 2 .
(5.46)

Remarkably, we see in figure 5.21 that the three quantities grow monotonically but tend to finite values in the limit t → ∞:

lim Re,t→∞ ∂u ∂z = 3 4 -σ 2 k z U S , (5.47) lim Re,t→∞ ∂v ∂z = 1 F h 1 -β 2 (k z F h U S ) 2 - 1 4 + χ 2 , (5.48) lim Re,t→∞ ∂b ∂z = 1 F 2 h .
(5.49)

This clearly demonstrates that the vertical shear cannot grow indefinitely in strongly stratified flows even for Re = ∞, contradicting the conjecture of [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]. We will see later whether or not the bounds (5.47-5.49) are sufficiently high to reach the threshold Ri = 1/4 necessary for the development of the shear instability. Before, we will conduct a parametric study in order to check if the maxima of the vertical gradients of u, v and b are always well predicted by the asymptotic formulas (5.41), (5.42) and (5.45).

As already remarked, ∂v/∂z and ∂b/∂z are one order higher than ∂u/∂z (see figure 5.21 and (5.41), (5.42), (5.45)). Equations (5.42) and (5.45) show that the quantities ∂v/∂z and ∂b/∂z both depend only on the product Λ = tΩ c except the viscous term in (5.42) but which is small. Therefore, they both reach their maxima when dΛ dt = 0.

(5.50)

The time of saturation is therefore given by T sat = k z F h U S t sat = (3/4γ) 1/3 where γ is defined in (5.26), i.e.

t sat = 1 k z F h U S 3 4γ 1/3 = 3Re 4k 2 z U 2 S 1/3 . (5.51)
At this time, the angular velocity is

Ω c = Ω sat = 2 3 1 + 2 9γ 2/3 . (5.52)
It is then straightforward to obtain the maxima of ∂v/∂z and ∂b/∂z:

max ∂v ∂z = M v = η F h 1 -β 2 (k z F h U S ) 2 - 1 4 + χ 2 η 2 + k z U S ReF 2 h χ ν 2 3 4γ 2/3 , (5.53) max ∂b ∂z = M b = η 2 F 2 h , (5.54)
where η is a function of the single parameter γ, namely

η = 2 3 1 + 2 9γ 2/3 3 4γ 1/3 .
(5.55)

In contrast, the viscous term in (5.41) governs the level of saturation of ∂u/∂z as seen in figure 5.21a. The maximum M u of ∂u/∂z is reached at t * sat = t sat . Since it is not possible to derive simple analytic expressions for t * sat and M u , these two quantities have been determined numerically using the asymptotic formulas (5.27) and (5.41).

Effect of k z

Figure 5.22 shows the evolution of ∂u/∂z, ∂v/∂z, and ∂b/∂z at the vortex center and z = l z /2 in numerical simulations for different vertical wavenumbers k z , for F h = 0.1, U S = 0.2, and Re = 6000. As k z increases, these quantities saturate earlier and at a higher level (figure 5.22a,c,e). The right plots (figure 5.22b,d,f ) show the same curves when each quantity is rescaled by its asymptotic maximum M u , M v (5.53), M b (5.54), the time t being rescaled by the corresponding saturation time t * sat or t sat . The different curves of ∂v/∂z and ∂b/∂z collapse very well with rescaled maxima close to unity and reached when t/t sat = 1. The time of saturation is thus proportional to k -2/3 z as seen in (5.51). Regarding ∂u/∂z (figure 5.22b), the level of saturation is correctly predicted by the asymptotics and the three rescaled curves reach their respective peaks at the same time t/t * sat 0.6 different from unity. The fact that the theory does not predict the exact saturation time is not surprising given the discrepancy observed in figure 5.21a. However, it is noticeable that the scaling of the saturation time as well as the maximum of ∂u/∂z are properly estimated.

Effect of U S

As for k z , figure 5.23 displays the evolution of ∂u/∂z, ∂v/∂z, and ∂b/∂z at r = 0, z = l z /2 in numerical simulations for different values of U S . Enhancing the shear amplitude again leads to higher saturation levels that are reached earlier. When rescaled (figure 5.23d,f ), (∂v/∂z)/M v and (∂b/∂z)/M b saturate around unity at t/t sat = 1. As seen in (5.51), the time of saturation thus scales like U -2/3 S . Regarding ∂u/∂z, the rescaling aligns the peaks around t/t * sat 0.6 and their levels are again close to unity as observed for k z .

Effect of Re

Figure 5.24 shows that increasing the Reynolds number delays the saturation of ∂u/∂z, ∂v/∂z and ∂b/∂z at r = 0, z = l z /2 and slightly increases the peak values. This is consistent with (5.51), since this equation shows that t sat is weakly affected by the Reynolds number because Re intervenes only to the power 1/3 instead of -2/3 for k z and U S .

The rescaling of ∂v/∂z and ∂b/∂z by their asymptotic maxima and of the time by t sat (figure 5.24d,f ) accurately collapses the evolution and saturation of these quantities. The prediction of the saturation of ∂u/∂z is less satisfying this time. Indeed, the curve corresponding to the highest Reynolds number Re = 10000 departs from the other ones (figure 5.24b). Nevertheless, the two other curves still saturate at t/t * sat 0.6 and their peaks are both close to unity as for k z and U S .

Effect of F h

Eventually, we examine the influence of the horizontal Froude number in figure 5.25. When F h is increased, from F h = 0.1 to F h = 0.5, ∂v/∂z and ∂b/∂z at the vortex center at z = l z /2 decrease (figure 5.25c,e) while ∂u/∂z increases (figure 5.25a). In addition, large oscillations are present for F h = 0.5. When the curves are rescaled as before, the collapse is less satisfactory than for the other parameters investigated previously, especially for the case F h = 0.5. Two reasons can explain these discrepancies. First, as F h increases, the parameter k z F h goes from 0.2 to 1. Thus, the hypothesis of the long-wavelength analysis, i.e. k z F h 1, is less and less valid as F h increases. Secondly, the law for the evolution of the angular velocity has been derived under the hydrostatic approximation, i.e. equivalently F h 1, in order to simplify the calculations. This assumption implies, in particular, that the internal waves excited at t = 0 are neglected. However, it is likely that these internal waves are responsible for the oscillations seen for F h = 0.5 in figure 5.25. To test this claim, it could be therefore interesting to derive the evolution of the angular velocity without making the hydrostatic approximation.

Evolution of the Richardson number

Having validated the asymptotics for the evolution of the vertical shear of the vortex and the vertical gradient of the buoyancy at the vortex center at z = l z /2, we now consider the evolution of the asymptotic Richardson number at this location: (5.56) where ∂u/∂z, ∂v/∂z and ∂b/∂z are given by (5.41), (5.42), and (5.45).

Ri = 1 F 2 h + ∂b ∂z -k z U S + ∂u ∂z 2 + ∂v ∂z 2 ,
In figure 5.4, we have plotted (5.56) for the two sets of parameters considered in section §5.3. The agreement with the numerical results is good for k z = π, F h = 0.1, while for k z = 3π/2, F h = 0.5, the initial decrease of the asymptotic Richardson number is well captured but subsequently, it re-increases instead of decreasing. As discussed above, this discrepancy comes probably from the oscillations due to internal waves that are present at F h = 0.5 in the numerical simulations. These oscillations are absent in the asymptotics because of the use of the hydrostatic approximation. Hence, the asymptotics are expected to be valid only if the Froude number is below F h 0.25.

To study the minimum value of (5.56) depending on the control parameters k z , U S , F h , Re, it can be first remarked that (5.41), (5.42), and (5.45) can be written in the form (5.61) This relation implies that min (Ri) < 0.25 only if Ri 0 0.12. Therefore, the Richardson number will go below the threshold Ri = 0.25 in the inviscid limit only if the initial Richardson number is already below this critical value, i.e. only if the ambient shear flow is unstable by itself.

∂u ∂z = k z U S F u (k z F h U S , Re b , t), (5.57) ∂v ∂z = 1 F h F v (k z F h U S , Re b , t), (5.58) ∂b ∂z = 1 F 2 h F b (k z F h U S ,
= 1 + F b (k z F h U S ) 2 (F u -1) 2 + F 2 v , ( 5
(Ri) = 2Ri 0 σ 2 + 1 4 2 + Ri 0 3 4 -χ 2 - β 2 Ri 0 2 .

Occurrence of the shear instability at F h = 0.5

Since the shear instability has been only observed in the numerical simulations performed at F h = 0.5, it is worth clarifying when it occurs at this particular Froude number, depending on the other control parameters k z , U S , Re.

To do so, and following the same idea as in §5.5.4, we have plotted in figure 5.27 a map discriminating stable cases from their unstable counterparts in the plane (Ri 0 , Re b ). For a given U S , the instability develops only if the buoyancy Reynolds number is sufficiently high and Ri 0 sufficiently low (see the circles for U S = 0.2). However, this marginal stability line varies with U S : it moves towards lower Ri 0 as U S increases. Indeed, for Re b = 1500, the instability develops for U S = 0.1 even if Ri 0 is as high as Ri 0 = 40, while for U S = 0.2 and U S = 0.3, it occurs only for Ri 0 10 and Ri 0 2, respectively. For U S = 0.4, no instability has been observed even at Ri 0 1. Therefore, at this Froude number, the occurrence of the Kelvin-Helmholtz instability cannot be predicted on the basis of Ri 0 and Re b only as suggested by the asymptotics. Nevertheless, there is no contradiction with the asymptotic predictions because of the effect of the internal waves in the moderately stratified regime.

Conclusion

We have performed direct numerical simulations of the evolution of an initially columnar vortex in an ambient shear flow in a strongly stratified fluid. The numerical results have been compared to the asymptotic analyses carried out in chapter 4.

The DNS show that the vortex is progressively bent in the direction of the shear flow but the vortex core also deviates in the orthogonal direction. The decay of its potential vorticity is enhanced in the regions of high shear. For F h = 0.5 and sufficiently high Reynolds number Re and vertical wavenumber k z , the Kelvin-Helmholtz instability is first triggered in the center of the vortex at the vertical levels where the ambient vertical shear is maximum. For the lower Froude numbers investigated F h = 0.1 and F h = 0.25, the Kelvin-Helmholtz instability has never been detected but the Richardson number reaches its minimum at the same location. We have therefore concentrated our efforts on the understanding of the evolution of the flow at this location thanks to the longwavelength asymptotic analysis performed for k z F h 1 in chapter 4. Besides, we have first shown that the short-time and the long-wavelength analyses account very well for the initial non-hydrostatic regime and the deformation of the vortex axis, respectively.

For F h 0.25, the long-wavelength asymptotic analysis turns out to predict very well the evolution of the vertical shear of the vortex in the spanwise direction ∂v/∂z and of the vertical gradient of the buoyancy ∂b/∂z at the vortex center and z = l z /2 for all the parameters investigated. The prediction for the streamwise shear of the vortex ∂u/∂z is less accurate but this is not too dramatic since it is one order smaller in (k z F h ) 2 than ∂v/∂z and ∂b/∂z.

The long-wavelength asymptotic analysis reveals that the growth of these three quantities is bounded owing to the conservation of potential vorticity and by viscous effects. Thereby, the maxima of ∂v/∂z and ∂b/∂z at leading order are of the form h is the buoyancy Reynolds number and Ri 0 = 1/(k z F h U S ) 2 the initial Richardson number. For small γ (i.e. near the inviscid limit), the functions F v and F b tend to constants close to unity, while for large γ, they behave as (5.63) at leading order.

M v = F v (γ)/F h and M b = F b (γ)/F 2 h , where γ = 1/(ReF 3 h k z U S ) = √ Ri 0 /
F v (γ) ∼ 1 γ 1/3 ∼ F h (k z U S Re) 1/3 , (5.62) F b (γ) ∼ 1 γ 2/3 ∼ F 2 h (k z U S Re) 2/3 ,
The Richardson number at the vortex center and z = l z /2 based on the asymptotic expressions of ∂u/∂z, ∂v/∂z, and ∂b/∂z is in good agreement with the DNS provided that the Froude number is small F h 0.25. This limitation comes from the fact that the hydrostatic approximation has been assumed in some parts of the long-wavelength analysis.

The minimum asymptotic Richardson number depends only on the two parameters Ri 0 and Re b . For large Re b , it can be lower than the critical value 1/4 necessary for the development of the shear instability only if the initial Richardson number is already smaller than this value.

This result puts into question the idea that the shear instability is easily triggered in strongly stratified shear flows. It also demonstrates that the mechanism proposed by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF], according to which the vertical shear should grow algebraically with time without bounds, is not valid for the flow studied herein. Because of the vertical coupling, neglected by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF], the flow actually decays in the regions of high shear, thereby limiting the maximum vertical shear attainable. This is due to potential vorticity conservation and enhanced viscous effects.

The shear instability has been only observed for F h = 0.5 and sufficiently high Re and k z . The role played by oscillations due to the transient internal waves excited at the start-up of the motion seems crucial for the development of the shear instability in this case. These oscillations are indeed of particularly large amplitudes for F h = 0.5 and high k z . Hence, the development of the shear instability for F h = 0.5 might be specific to the particular configuration and initial conditions chosen here. However, in stratified turbulent flows, although the structures can be continuously vertically sheared, it is unlikely that a coherent structure will be impulsively sheared at a given time. As seen in section §5.4, this impulsive start-up leads to a transient non-hydrostatic regime and to the generation of large internal waves when F h is not small. In contrast, the shearing process is expected to be much more continuous and smooth in stratified turbulent flows, limiting the generation of internal waves.

Chapter 6

Conclusions and perspectives

The purpose of the present dissertation has been to understand the different facets of the response of an initially vertical vortex embedded in an environmental vertically sheared flow, in presence of a strong stratification. To do so, we have considered the case of a model flow comprising a straight Lamb-Oseen vortex at t = 0 subjected to a flow uniform in the horizontal but varying sinusoidally along the vertical. Under the effect of the differential advection due to the shear flow, the vortex is progressively decorrelated along the vertical, this process being at the origin of the existence of localized regions of intense vertical shear. This investigation therefore enables testing Lilly's conjecture [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]) that the vertical shear should grow algebraically with time, leading to the generation of small scales through the shear instability for large enough times. More generally, the present study has aimed at clarifying many facets of the dynamics of the vortex: its displacements, its structure evolution, and notably the variations of the density as well as the vertical shear of the horizontal velocity in the vortex core.

In this regard, a first approach has been based on the analysis of the global budgets of energy and enstrophy (see chapter 3). We have been able to determine an estimation of the minimum shear layer thickness that is associated with the peak value attained by the global enstrophy of the vortex.

Remarkably, the scaling which has been obtained for the maximum enstrophy in terms of the control parameters indicates that it is proportional to the Reynolds number to the power 2/3, suggesting that the minimum Richardson number is in turn inversely proportional to Re 2/3 . This result questions former estimations stating that the minimum Richardson number scales as the inverse of the buoyancy Reynolds number ReF 2

h [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Deloncle | Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation[END_REF]. In fact, the exponent 2/3 comes from the presence of the ambient shear flow: the dynamics is not the same as in freely decaying flows. This is not the only remarkable result of this investigation. Indeed, the theoretical scalings derived in chapter 3 account for the effect of all the control parameters except the Froude number: the interpretations presented in this chapter do not predict any dependency with respect to this parameter although it clearly influences the decay of the kinetic energy as well as the variations of enstrophy.

In order to overcome this paradox, we have carried out local investigations of the dynamics based on asymptotic analyses in chapters 4 and 5. First, the study of the short-time dynamics has shown that the initial response of the vortex to the ambient shear flow is non-hydrostatic regardless of the level of stratification. After this transient stage however, the hydrostatic balance is progressively recovered. The initial nonhydrostatic regime arises because of the excitation of internal waves at the beginning of the flow evolution.

Another central result of these analyses is that the core angular velocity of the vortex decays because of dynamic and viscous effects. Notably, this decay is not uniform in the vertical: it is maximum in the regions of maximum vertical shear whereas it cancels at the locations where the vertical shear is close to zero. The decrease of the angular velocity actually modulates the growth of the vertical shear of the horizontal velocity which, in turn, cannot grow indefinitely as speculated by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]: it has an upper bound accurately estimated by the asymptotics. A key factor in this process of modulation is the vertical velocity, which had been ignored by Lilly. The asymptotics have not only enabled us to obtain analytic predictions regarding the saturation of the different quantities but also for their global temporal variations, thus providing a comprehensive description of the dynamics.

Nevertheless, these analytic results have been obtained in the asymptotic limit of long vertical wavelength and strong stratification, hence allowing to use the hydrostatic approximation for the sake of simplicity. Their use in the moderately stratified regime has shown that they still account for the observed dynamics qualitatively, but are no longer valuable quantitatively. In particular, we have shown that the dependence of the vertical shear of horizontal velocity with respect to the Froude number is quite complex. Nevertheless, the effect of this parameter has been understood thanks to the local analyses presented in chapters 4 and 5, since the formulas derived in these chapters capture quite well the physics provided that the stratification is strong. Hence, these advances offer the possibility of explaining the effect of the Froude number on the variations of the volume-integrated energy and enstrophy, which remained unexplained in chapter 3, by integrating the local equations obtained via the asymptotics over the entire fluid domain.

Beyond the effect of the Froude number, the analytic predictions of chapter 4 yield an accurate expression for the local Richardson number in chapter 5. A major result is that, for strong stratification, its minimum is a function of time only, and of two well-identified parameters: the Richardson number associated with the ambient shear flow only, and the buoyancy Reynolds number. In this case, we have shown that the shear instability can be triggered only if the ambient shear flow is itself unstable in the sense of the Miles & Howard criterion [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF][START_REF] Howard | Note on a paper of John W. Miles[END_REF]. Thus, our results prove that the process of vertical decorrelation itself is insufficient to trigger the shear instability in strongly stratified fluids, in strong contrast with Lilly's hypothesis. However, we have shown that the shear instability develops for moderate stratification and large enough vertical wavenumber even if the advection by the ambient shear flow is not that strong.

Globally, the present investigations offer a variety of perspectives for future works on related issues. First, the present long-wavelength analysis could be improved and probably extended to the moderately stratified case by taking into account the wave terms that have been neglected herein when tackling the problem of the evolution of the vortex structure. More generally, the asymptotic method presented in chapter 4 could be replicated to other model cases, for example by assuming that the ambient shear flow varies linearly along the vertical, or by considering that a rapid background rotation is present over the stratification. The latter configuration should be interesting because there would no longer be any excitation of internal waves in this case. Thus, the dynamics should be significantly different from that observed in the purely stratified case. In addition, our results prove that the decay of the core angular velocity is of utmost importance in order to model the evolution of the vertical shear of horizontal velocity properly. This could have important implications for two-layer models (or multi-layer models more generally). Indeed, they usually ignore this effect and, as such, may elude central physical processes, making them insufficiently representative of the real physics.

Appendix B

Derivation of the streamfunction ψ w for small time

In this appendix, we derive analytically the streamfunction ψ w associated to internal waves for small time. This will allow us to obtain the asymptotic expressions of C w and S w for t 1. It is first useful to rewrite the right-hand side of the equation for ψ p (4.105) in the form Accordingly, the solution ψ w of (4.106) is expressed as: where ψ pc is defined in (4.119) and where it has been also imposed that ψ w1 → 0 as r → +∞. In turn, (B.9) and (B.10) yield where the constants of integration have been chosen so that ψ w3 vanishes at infinity. At fourth order, only ∇ 2 r ψ w4 will be needed: it is not necessary to derive ψ w4 . From (B.11) and (B.14), we have directly:

ψ w = -
ψ w2 = A 4 rΩ, ψ w3 = -
∇ 2 r ψ w4 = ζ w4 48F 2

Résumé détaillé de la thèse

Contexte Cette thèse étudie un processus de décorrélation d'un tourbillon initialement rectiligne par un écoulement externe cisaillé sinusoïdalement selon la verticale, dans un fluide stratifié. Il a été proposé qu'un tel mécanisme devrait déclencher des instabilités de cisaillement et ainsi contribuer à la production de petites échelles en turbulence fortement stratifiée, que l'on rencontre dans l'atmosphère et les océans dans une gamme d'échelles où la force de Coriolis est négligeable.

Plus précisément, [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] a émis l'idée qu'en présence d'une forte stratification, toute structure cohérente soumise à un champ de cisaillement de la forme U (z), avec z la coordonnée verticale, devrait présenter un cisaillement vertical de vitesse horizontale augmentant avec le temps en raison de l'advection différentielle selon la verticale. En d'autres termes, la forte stratification tend à "scinder" l'écoulement en un empilement de couches dont les dynamiques sont indépendantes. En effet, la conjecture de Lilly suppose que la structure est advectée sans génération de vitesse verticale et que tout éventuel couplage vertical de son champ de vitesse est négligé. L'équation mathématique modélisant cette dynamique s'écrit :

∂u h ∂t = -U (z)
∂u h ∂x où l'on a supposé que l'écoulement ambiant est porté par un axe e x associé à la coordonnée x, avec u h le champ de vitesse (horizontal) de la structure considérée. La solution est u h (x, y, z, t) = u h0 (x -U (z)t, y, z) avec u h0 le champ de vitesse u h initial, si bien que le cisaillement vertical de vitesse horizontale ∂u h ∂z = ∂u h0 ∂z -t dU dz ∂u h0 ∂x croît algébriquement avec le temps t. On s'attend donc à une décroissance continue du nombre de Richardson, celui-ci étant inversement proportionnel au carré du cisaillement vertical de vitesse horizontale. La transition sous le seuil critique 1/4, condition nécessaire au déclenchement d'une instabilité de cisaillement, devrait ainsi être inconditionnelle et l'instabilité devrait donc se développer après un certain laps de temps dépendant notamment de l'intensité de l'advection U .

Ce mécanisme, bien que fréquemment invoqué dans les études de turbulence fortement stratifiée, n'a pourtant jamais été testé en pratique. Cette thèse vise à le caractériser dans le cas d'un écoulement modèle comprenant un tourbillon de Lamb-Oseen soumis à un écoulement ambiant de profil de vitesse U (z) = U S sin(k z z).

Résultats La première partie est dédiée à l'analyse des évolutions des énergie cinétique et enstrophie totales du tourbillon calculées au moyen de Simulations Numériques 143 Directes (DNS), en fonction des paramètres de contrôle. Cette étude montre que la dynamique est différente de celle des écoulements libres non cisaillés : du fait de la présence de l'écoulement ambiant, la condition de saturation entre les termes d'étirement et de dissipation dans le bilan d'enstrophie globale implique que l'enstrophie maximale du tourbillon sature comme Re 2/3 , Re étant le nombre de Reynolds, et non pas proportionnellement à Re comme pour les écoulements sans cisaillement ambiant. Néanmoins, cette condition de saturation ne rend pas compte de l'effet observé de la stratification.

Afin de comprendre cet effet, la dynamique locale du tourbillon a fait l'objet de deux analyses asymptotiques présentées dans une deuxième partie. Une première étude pour temps courts prouve que la réponse initiale du tourbillon est non-hydrostatique quelle que soit la stratification. Il est également démontré que, de manière remarquable, la stratification n'influence en réalité que très peu la dynamique initiale. Une deuxième analyse pour grande longueur d'onde fournit les équations qui décrivent l'évolution de la vitesse angulaire du tourbillon et des déformations de son axe. L'excitation d'ondes internes au début de l'évolution est à l'origine du régime non-hydrostatique initial. Le tourbillon est principalement advecté dans la direction du cisaillement ambiant mais aussi perpendiculairement du fait de son auto-induction. Sa vitesse angulaire décroît à cause d'effets dynamique et visqueux. La décroissance dynamique est dûe à un resserrement des lignes isopycnes dans le coeur du tourbillon, qui implique une atténuation de la vorticité verticale afin que la vorticité potentielle se conserve. Enfin, un champ de vitesse associé à un mode azimuthal m = 2 est généré, indiquant que le tourbillon devient légèrement elliptique.

Dans une troisième et dernière partie, des DNS révèlent que l'instabilité de cisaillement se développe seulement lorsque la stratification est modérée et la longueur d'onde du cisaillement ambiant suffisamment petite. Les résultats numériques sont comparés aux prédictions asymptotiques. Les déformations du tourbillon sont prédites avec précision. Les ondes internes influencent fortement le mouvement du tourbillon au début de l'évolution mais leur effet s'estompe rapidement. L'axe du tourbillon ne se déforme pas de façon purement sinusoïdale : il subit une torsion verticale. Par ailleurs, les évolutions du cisaillement vertical de vitesse horizontale et du gradient vertical de flottabilité sont prédites de manière fine et exhaustive par l'analyse asymptotique pour grande longueur d'onde lorsque le nombre de Froude est petit, c'est-à-dire dans la limite fortement stratifiée. Le nombre de Richardson asymptotique admet un minimum qui n'est pratiquement jamais inférieur au seuil critique 1/4 nécessaire au déclenchement de l'instabilité de cisaillement. La saturation du cisaillement vertical est dûe au déclin accéléré du tourbillon dans les régions où le cisaillement ambiant est maximal.

Conclusions Les résultats asymptotiques dans la limite des grandes longueurs d'onde permettent d'aboutir à une expression analytique du nombre de Richardson local, qui permet de le prédire finement lorsque la stratification est forte. On démontre, grâce à cette expression, que ce nombre ne peut être inférieur au seuil critique 1/4 nécessaire pour que l'instabilité de cisaillement se développe, à moins que l'écoulement ambiant ne soit lui-même instable. Ce résultat est cohérent avec les simulations numériques, cellesci n'ayant jamais permis d'observer l'instabilité pour de petits nombres de Froude.

Ce comportement est lié à la saturation du cisaillement vertical de vitesse horizontale, qui résulte du déclin accéléré du tourbillon dans les régions où |dU/dz| est maximal. Cette saturation entre en contradiction directe avec l'hypothèse formulée par Lilly. On soulignera que la décroissance dynamique de la vitesse angulaire du tourbillon suffit à provoquer une saturation des cisaillements verticaux.

Au vu de ces résultats pour un tourbillon, il est conclu que les processus de décorrélation d'une structure cohérente par un champ de cisaillement ambiant semblent insuffisants, en fluide fortement stratifié, pour déclencher des instabilités de cisaillement, et donc pour participer à la génération de petites échelles, contrairement à la conjecture de Lilly (1983). vortex by a strongly stratified shear flow decorrelation, vortex, shear, stratification s the vertical decorrelation of an initially straight vortex by an ambient sinusoidal shear been conjectured that such process should trigger shear instabilities and, as such, contribute to the generation of small scales in strongly stratified turbulence, which is encountered in the atmosphere and the oceans in an intermediate range of re Coriolis effects are negligible. The first part analyses the evolutions of the total energies and enstrophy of the vortex in Direct Numerical ics differs from free because of the presence of the ambient shear flow, the balance between stretching and dissipation terms in the global enstrophy budget implies that the enstrophy of the vortex is proportional to the to the power 2/3. However, such lance does not account for the observed effect , the local dynamics has been investigated for short time or longanalysis provides the governing equations for the velocity of the vortex and the deformations of its axis. mainly show that the angular velocity decays because dynamic and viscous effects. from the squeezing of the isopycnals in the vortex core which implies a decrease of the vertical vorticity to conserve potential vorticity. In a third part, the DNS show that the shear inst only triggered for a moderate of the vertical shear of horizontal velocity and of the vertical buoyancy gradient asymptotics in the strongly stratified limit Richardson number always threshold 1/4 necessary for instability. The saturation of the vertical shear is due to the accelerated decay of the vortex in the regions of maximum ambient shear. These results suggest that the shear instability is not easily triggered by decorrelation processes in strongly stratified flows in contradiction with previous conjectures. l'instabilité de cisaillement est difficilement déclenchée par des processus de décorrélation dans les écoulements fortement stratifiés, contredisant ainsi les conjectures précédemment formulées dans la littérature. the evolution of the angular the vortex and the deformations of its axis. We mainly show that the angular velocity decays because of dynamic and viscous effects. The dynamic decay results from the squeezing of the isopycnals in the vortex core a decrease of the vertical vorticity to

In a third part, the DNS show that the shear instability is moderate stratification. The evolutions of the vertical shear of horizontal velocity and of the gradient are finely captured by the in the strongly stratified limit. The asymptotic Richardson number always remains larger than the critical necessary for the development of the shear The saturation of the vertical shear is due to the accelerated decay of the vortex in the regions of maximum These results suggest that the shear instability is not easily processes in strongly stratified flows in contradiction with previous conjectures.

Figure 1

 1 Figure 1.2: Generation of turbulence by a grid. The turbulent part of the flow, visible in the downstream region, is chaotic, disordered, without any clear trend in the motion of the fluid parcels in strong contrast with the far field laminar flow of the upstream region. Reproduced from Van Dyke (1982).

Figure 1 . 4 :

 14 Figure 1.4: Illustration of the development of the zigzag instability in the case of a pair of counter-rotating vortices. Reproduced from the experiments presented in Billant and Chomaz (2000a).

Figure 1

 1 Figure 1.5: Example of experiment showing the KH instability of an interface between two fluids of equal depth and different densities in relative motion. Reproduced from Thorpe (1968).

Figure 2

 2 Figure 2.1: A typical profile of the Brunt-Väisälä frequency in the ocean. The measurement was performed in the North Atlantic ocean, near 28 • North, 70 • West. It is reproduced from Gill (1982).

z

  = l z /4: no trunc. z = lz/4: trunc. z = l z /2: no trunc. z = l z /2: trunc.

Figure 2 . 2 :

 22 Figure2.2: Evolution of the potential vorticity maxima computed in the planes z = l z /4, without any aliasing removal (solid line) and with a dealiasing (open circles), and z = l z /2, without any aliasing removal (dashed line) and with a dealiasing (triangles). The elliptic dealiasing cancels the top one-third modes following[START_REF] Orszag | On the elimination of aliasing in finite-difference schemes by filtering highwavenumber components[END_REF]. The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000.

Figure 2

 2 Figure 2.3: Displacement of the vortex with respect to its initial position, in the plane z = l z /4. Four trajectories are reported, respectively corresponding to η = 90% (black dash-dotted line), η = 50% (black dashed line), η = 10% (black solid line) and η = 5% (grey open circles). The sampling times range from t = 0 up to t = 35. The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000. The horizontal dimensions of the domain are l x = l y = 18.

Figure 2 . 4 :

 24 Figure2.4: Displacement of the vortex with respect to its initial position, in the plane z = l z /4. Three trajectories are reported, respectively corresponding to l x = l y = 18 (triangles), l x = l y = 24 (circles), and l x = l y = 30 (squares). The sampling times range from t = 0 up to t = 40. The control parameters are F h = 0.5, k z = π, U S = 0.1, and Re = 6000. The displacements are those of the potential vorticity centroid obtained with η = 5%.

  Figure 3.1: Left column: three-dimensional contours of the potential vorticity at different times for F h = 0.5, k z = π, U S = 0.2, Re = 6000. Right column: corresponding vertical cross-sections of the shear √ S z = (∂u/∂z) 2 + (∂v/∂z) 2 (color) and of the total density ρ t (black contour lines) in the plane y = 9. The times shown are (a,b) t = 12, (c,d) t = 24, (e,f ) t = 36. In (a,c,e), the isocontours correspond to 20% (yellow) and 60% (red) of the initial maximum value.

Figure 3

 3 Figure 3.2: Vertical profiles of (a) the streamwise velocity component ū and (b) the spanwise velocity component v of the mean flow. The profiles are shown at four different times t = 0 (grey dash-dotted line), t = 10 (black dash-dotted line), t = 20 (grey dashed line) and t = 30 (black dashed line). The control parameters are F h = 0.5, k z = π/2, U S = 0.2, and Re = 6000. The left picture (a) clearly shows that the departures from the initial sinusoidal profile are insignificant, whereas the right picture (b) indicates that a spanwise component grows as time evolves, but remains small relative to the streamwise flow.

  Figure3.3: Evolutions of (a) the total horizontal kinetic energy (black solid line), the mean flow kinetic energy (black dash-dotted line), the vortex horizontal kinetic energy (black dashed line), the vertical kinetic energy (grey solid line), and the potential energy (grey dashed line) and of (b) the total horizontal enstrophy (black solid line), the enstrophy of the mean flow (black dash-dotted line), the vortex horizontal enstrophy (black dashed line), and the vertical enstrophy (grey solid line). The parameters are F h = 0.5, k z = π, U S = 0.2 and Re = 6000.

Figure 3

 3 Figure 3.4: Source terms in (a) the kinetic energy budget (3.11) and in (b) the enstrophy budget (3.21). The parameters are the same as in figure 3.3.

Figures 3 .

 3 Figures 3.4a,b show the source terms of the kinetic energy equation (3.11) and of the enstrophy equation (3.21) of the vortex, for the same run as presented in figure 3.3.As seen in figure3.4a, the kinetic energy of the vortex decreases because of three main effects. There are first a transfer to the potential energy (C * ) and also a transfer to the shear flow (T ). Later, these two terms decrease while the dissipation due to the vertical shear of the horizontal velocity (D * z ) increases and saturates at t 40. The other dissipative terms are very small (D * h + D w ). Figure3.4b shows that the enstrophy of the vortex grows because of the stretching and tilting, first by the shear flow ( S * ), and then by the vortex on itself (S * ) also. In contrast, the baroclinic term (B * ), the transfer to the enstrophy of the vortex (T ), and dissipation (D * ) are all negative. All these terms saturate between t 20 and t 30. The dominant terms are the stretching S * and the dissipation D * while the stretching S * contributes also significantly to the budget but grows later and is at best of the order of S * . Therefore, in order to estimate the maximum enstrophy of the vortex (figure3.3a), we can, as a first approach, approximate the enstrophy equation of the vortex in the form dZ * dt S * + D * .(3.30)

Figure 3 .

 3 Figure 3.5a shows the temporal evolution of the kinetic energy of the vortex E * k for several wavenumbers k z , for F h = 0.5, U S = 0.2 and Re = 6000. The decay is faster as k z increases. When the time is rescaled by τ v (figure 3.5b), the three curves collapse quite well confirming (3.41).Similarly, figure 3.6 displays the evolution of the vortex kinetic energy for various U S , for F h = 0.5, k z = π and Re = 6000. Again, the rescaling of t by τ v (figure 3.6b) leads to a self-similar decay in agreement with (3.41).The effect of the Reynolds number is tested in figure3.7 for F h = 0.5, k z = π and U S = 0.2. It is seen that the evolution of E * k weakly depends on the Reynolds number, especially the initial part t 20. When t is rescaled by τ v (figure 3.7b), the collapse of the three curves remains satisfactory. This is because τ v depends on Re with the power law exponent 1/3. Moreover, the relative variations of Re in figure3.7 are smaller than for the two previous parameters.All these results are consistent with (3.41). Remarkably, this relation is not only valid for large times but also for short times. This suggests that the transfer terms T and C * , although dominant at the beginning, do not have a so important effect.

Figure 3 Figure 3 . 7 :

 337 Figure 3.6: Kinetic energy E * k of the vortex for F h = 0.5, k z = π, Re = 6000, for U S = 0.2 (dash-dotted line), U S = 0.3 (dashed line) and U S = 0.4 (solid line), as a function of (a) the time t and (b) the rescaled time t/τ v .

Figure 3

 3 Figure3.9: Enstrophy of the vortex for F h = 0.5, k z = π, Re = 6000, for U S = 0.2 (dash-dotted line), U S = 0.3 (dashed line) and U S = 0.4 (solid line). The unscaled enstrophy Z * is plotted against t in (a) and the rescaled enstrophy δ 2 Z * is plotted against the rescaled time t/τ v in (b).

Figure 3

 3 Figure3.10: Enstrophy of the vortex for F h = 0.5, k z = π, U S = 0.2, for Re = 4000 (dash-dotted line), Re = 5000 (dashed line) and Re = 6000 (solid line). The unscaled enstrophy Z * is plotted against t in (a) and the rescaled enstrophy δ 2 Z * is plotted against the rescaled time t/τ v in (b).

Figure 4

 4 Figure 4.1: (a) Sketch of the flow configuration studied in this dissertation. The ambient shear flow U (z)e x is represented by the black vectors. The red and yellow contours show two levels of the potential vorticity of the vortex at t = 0. (b) Potential vorticity contours at t = 20 obtained from a direct numerical simulation for the Froude and Reynolds numbers F h = 0.1, Re = 6000 and the shear flow U = U S sin(k z z) with U S = 0.2 and k z = π. See chapter 5 for the definitions of these parameters and the details of the numerical simulations.

Figure 4

 4 Figure 4.2: Variations of (a) β and (b) k z β against the radius r for different values of the vertical wavenumber: k z = 0 (grey solid line), k z = π/6 (grey dashed line), k z = π/2 (black dash-dotted line), k z = π (black dashed line) and k z = 3π/2 (black solid line).

  4.70)-(4.71), (4.74)-(4.75), and (4.77)-(4.78) give the complete expressions of w 0 and b 0 :

Figure 4 . 3 :Figure 4

 434 Figure 4.3: Streamfunctions ψ wc (a) and ψ ws (b) plotted against r at t = 0 (grey solid line), t = 4 (black dash-dotted line), t = 8 (black dashed line) and t = 16 (black solid line), at z = l z /4 where Ω = [1 -exp(-r 2 )] /r 2 regardless of τ . The Froude number is F h = 0.1.

Figure 4

 4 Figure 4.5: Variables (a) A and (b) ϑ plotted against F h , at z = l z /4 where Ω = [1 -exp(-r 2 )] /r 2 regardless of τ .

Figure 4 . 6 :

 46 Figure 4.6: Forcings C w (grey lines) and S w (black lines) for (a) F h = 0.1 and (b) F h = 0.5, at z = l z /4 where Ω = [1 -exp(-r 2 )] /r 2 regardless of τ .

Figure 4 . 7 :

 47 Figure 4.7: Streamfunctions ψ2b (a), ψ2s (b), and ψ2c (c) rescaled by r2 as functions of r. Note that they all rapidly vanish outside the vortex core.

Figure 4

 4 Figure4.9: Angular velocity profile Ω computed numerically from (4.178) (grey solid lines) and predicted by (4.186) (black dashed lines) at different times for γ = 0.1 (a) and γ = 1 (b), as a function of r. The Lamb-Oseen profile (4.190) is also represented by black dash-dotted lines. The times shown are: (a) t = 0, t = 2.65, t = 7.50, t = 13.8, t = 21.2 and (b) t = 0, t = 0.265, t = 0.750, t = 1.38, t = 2.12, from top to bottom.

Figure 4

 4 Figure 4.10: Evolution of the angular velocity at the vortex center Ω c as a function of t1/3 predicted by (4.188) (black solid lines) and computed numerically from (4.178) (grey circles) for γ = 0.1, γ = 1, and γ = 10 (top to bottom).

  .206) where ∆ X = U S t+ X is the total displacement of the vortex in the streamwise direction. Such phenomenological model can be obtained from (4.194)-(4.195) by assuming Re → ∞ and by setting ϑ = ∂C w /∂t = ∂S w /∂t = 0, i.e. by neglecting the waves effect and the three-dimensional correction of the advection by the base flow. The model (4.205)-(4.206) has been introduced for example by

Figure 5 . 1 :

 51 Figure 5.1: Left column: three-dimensional contours of the potential vorticity at different times for F h = 0.1, k z = π, U S = 0.2 and Re = 6000. Right column: corresponding horizontal cross-sections in the plane z = l z /4 where the advection is the most intense. The times shown are (a,b) t = 4, (c,d) t = 13, (e,f ) t = 26. In (a,c,e), the isocontours correspond to 20% (yellow) and 60% (red) of the initial maximum value.

Figure 5

 5 Figure 5.2: Same as figure 5.1 except that F h = 0.5, k z = 3π/2. The times shown are (a,b) t = 4, (c,d) t = 13, (e,f ) t = 26, (g,h) t = 36.

Figure 5

 5 Figure 5.3: Vertical cross-sections of the shear √ S z (color) and of the total density ρ t (black contour lines) in the plane y = 9, for F h = 0.1, k z = π (left column) and F h = 0.5, k z = 3π/2 (right column), for U S = 0.2 and Re = 6000. The times shown are (a,b) t = 4, (c,d) t = 13, (e,f ) t = 26, (g,h) t = 36.

Figure 5

 5 Figure 5.4: Minimum Richardson number as a function of time for U S = 0.2, Re = 6000, and F h = 0.1, k z = π (black solid line) and F h = 0.5, k z = 3π/2 (grey solid line) from the DNS. The asymptotic prediction (5.56) has been also plotted for both cases F h = 0.1, k z = π (black dashed line) and F h = 0.5, k z = 3π/2 (grey dashed line). The horizontal black dash-dotted line shows the critical value Ri c = 0.25.

Figure 5 . 5 :Figure 5

 555 Figure 5.5: Vertical cross-sections of the buoyancy b at t = 26 in the planes y = 9 (a) and x = 9 (b) for F h = 0.5, k z = 3π/2, U S = 0.2, and Re = 6000. The black contours represent the lines where Ri = 0.25.

Figure 5

 5 Figure 5.7: Time evolutions of the global kinetic energy Ēk + E * kh + E * kz (black solid line), the mean flow kinetic energy Ēk (black dash-dotted line), the vortex horizontal kinetic energy E * kh (black dashed line), the vortex vertical kinetic energy E * kz (grey dashed line) and the potential energy E p (grey solid line) for (a) F h = 0.1, k z = π and (b) F h = 0.5, k z = 3π/2 for U S = 0.2, Re = 6000.

Figure 5 .Figure 5

 55 Figure5.7 presents the evolutions of the energies integrated over the whole computational domain for the two simulations for k z = π, F h = 0.1, and k z = 3π/2, F h = 0.5

Figure 5 . 9 :

 59 Figure 5.9: Horizontal cross-sections in the plane z = l z /2 at t = 0.4 of the vertical velocity (a,b) and buoyancy (c,d) computed asymptotically (a,c) and numerically (b,d).The contour interval is 0.002 for w and 0.001 for b. The contour at the level zero is not drawn. The negative iso-contours are colored in grey. The parameters are F h = 0.5, k z = π, U S = 0.2, and Re = 6000.

Figure 5 Figure 5

 55 Figure5.10: Same as figure 5.9 but for t = 1.6. The contour interval is 0.02 for w and 0.05 for b.

Figure 5

 5 Figure5.12: Temporal evolutions of (a) w m and (b) b m computed asymptotically (black lines) and numerically (grey lines) for F h = 0.5, U S = 0.2, and Re = 6000, for k z = π/2 (dash-dotted lines), k z = π (dashed lines) and k z = 3π/2 (solid lines).

Figure 5

 5 Figure5.13: Temporal evolutions of (a) the rescaled maximum vertical velocity w m /U S and (b) the rescaled maximum buoyancy b m /U S computed analytically (black lines) and numerically (grey lines) for F h = 0.5, k z = π, and Re = 6000, for U S = 0.02 (dash-dotted lines), U S = 0.1 (dashed lines), and U S = 0.2 (solid lines).

Figure 5

 5 Figure 5.14: (a,b) Comparison between the vertical vorticity centroids ∆x ζ c and ∆y ζ c (grey dashed lines with crosses) and the asymptotic predictions for the displacements ∆x (a) and δy (b) (black solid lines) in the plane z = l z /4. (c) displays the same comparison as in (b) except that the grey dashed line marked with open circles now represents the potential vorticity centroid ∆y Π c . (d) displays a close-up view of the initial evolution of ∆y ζ c (grey crosses), ∆y Π c (grey open circles), and δy (black solid line). The asymptotic prediction in the absence of internal waves, i.e. by setting C w = S w = 0 in (5.22)-(5.23) has been also reported (black dashed line). The shorttime asymptotic prediction for δy has been plotted as well (grey solid line). The parameters are F h = 0.5, k z = 0.3, U S = 0.2, and Re = 6000.

Figure 5 .

 5 Figure 5.15: Horizontal cross sections in the plane z = l z /2 of (a,b,c) the vertical velocity calculated asymptotically and (d,e,f ) in the DNS for F h = 0.5, k z = 0.3, U S = 0.2, and Re = 6000 at t = 2 (a,d), t = 6 (b,e) and t = 18 (c,f ).

Figure 5 .

 5 Figure 5.16: Horizontal cross sections in the plane z = l z /2 of (a,b,c) the buoyancy calculated asymptotically and (d,e,f ) in the DNS for F h = 0.5, k z = 0.3, U S = 0.2, and Re = 6000 at t = 2 (a,d), t = 6 (b,e) and t = 18 (c,f ).

Figure 5 .

 5 Figure 5.17: Evolution of the maximum values of (a) the vertical velocity w m and (b) the buoyancy b m in the DNS (grey dashed lines with open circles) and predicted by the asymptotics (black solid lines), for F h = 0.5, k z = 0.3, U S = 0.2, and Re = 6000.

Figure 5 .

 5 Figure 5.19: Evolution of the maximum buoyancy b m in the DNS (grey dashed lines with grey open circles) and from the asymptotic predictions (black solid lines) for Re = 6000, for (a) F h = 0.1, k z = 1.5, U S = 0.2, (b) F h = 0.1, k z = 0.3, U S = 0.2, (c) F h = 0.5, k z = 0.1, U S = 0.2, and (d) F h = 0.5, k z = 0.3, U S = 0.1.

Figure 5 Figure 5

 55 Figure5.20: Evolution of (a) ∂u/∂x, (b) ∂v/∂x, and (c) ∂b/∂x at r = 0, z = l z /2 from the DNS (grey solid lines) and predicted from the asymptotics (black solid lines), for F h = 0.1, k z = π, U S = 0.2, Re = 6000, i.e. γ 0.27.

Figure 5

 5 Figure 5.22: Evolution of ∂u/∂z (a,b), ∂v/∂z(c,d), and ∂b/∂z (e,f ) at r = 0, z = l z /2, for F h = 0.1, U S = 0.2, and Re = 6000, for k z = 2 (dash-dotted lines), k z = π (dashed lines), and k z = 4 (solid lines). The quantities ∂u/∂z, ∂v/∂z, ∂b/∂z are plotted against the physical time t in (a), (c), (e), respectively, whereas (b), (d), (f ) show (∂u/∂z)/M u against t/t * sat , and (∂v/∂z)/M v , (∂b/∂z)/M b against t/t sat , respectively.

Figure 5

 5 Figure 5.23: Evolution of ∂u/∂z (a,b), ∂v/∂z(c,d), and ∂b/∂z (e,f ) at r = 0, z = l z /2, for F h = 0.1, k z = 2, and Re = 6000, for U S = 0.2 (dash-dotted lines), U S = 0.3 (dashed lines), and U S = 0.4 (solid lines). The quantities ∂u/∂z, ∂v/∂z, ∂b/∂z are plotted against the physical time t in (a), (c), (e), respectively, whereas (b), (d), (f ) show (∂u/∂z)/M u against t/t * sat , and (∂v/∂z)/M v , (∂b/∂z)/M b against t/t sat , respectively.

Figure 5

 5 Figure5.24: Evolution of ∂u/∂z (a,b), ∂v/∂z (c,d), and ∂b/∂z (e,f ) at r = 0, z = l z /2, for F h = 0.1, k z = 2, and U S = 0.2, for Re = 4000 (dash-dotted lines), Re = 6000 (dashed lines), and Re = 10000 (solid lines). The quantities ∂u/∂z, ∂v/∂z, ∂b/∂z are plotted against the physical time t in (a), (c), (e), respectively, whereas (b), (d), (f ) show (∂u/∂z)/M u against t/t * sat , and (∂v/∂z)/M v , (∂b/∂z)/M b against t/t sat , respectively.

Figure 5

 5 Figure 5.25: Evolution of ∂u/∂z (a,b), ∂v/∂z(c,d), and ∂b/∂z (e,f ) at r = 0, z = l z /2, for k z = 2, U S = 0.2, Re = 6000, for F h = 0.1 (dash-dotted lines), F h = 0.25 (dashed lines), and F h = 0.5 (solid lines). The quantities ∂u/∂z, ∂v/∂z, ∂b/∂z are plotted against the physical time t in (a), (c), (e), respectively, whereas (b), (d), (f ) show (∂u/∂z)/M u against t/t * sat , and (∂v/∂z)/M v , (∂b/∂z)/M b against t/t sat , respectively.

Figure 5

 5 Figure 5.26: Contours of the minimum Richardson number as a function of its initial value Ri 0 and the buoyancy Reynolds number Re b . The grey domain corresponds to Ri 0 < 0.25 i.e. where the shear flow is unstable by itself. The black solid line shows the threshold min (Ri) = 0.25.

  .60) meaning that it is a function of only three parameters: k z F h U S , Re b , and t. This greatly simplifies its study. Thereby, figure5.26 shows the minimum of the Richardson number over time as a function of the initial Richardson number Ri 0 = 1/(k z F h U S ) 2 and Re b . Only the domain Ri 0 > 0.25 is of interest since below the ambient shear is unstable by itself. The threshold min (Ri) = 0.25 is indicated by the black solid line. For Re b ∼ 100, this line is tangent and already very close to the boundary Ri 0 = 0.25. When Re b is increased, it quickly crosses the limit Ri 0 = 0.25 meaning that min (Ri) 0.25 only if already initially, Ri 0 0.25.The limit for Re b 1 can be studied by considering the expressions (5.47), (5.48) and (5.49) of ∂u/∂z, ∂v/∂z and ∂b/∂z when Re → ∞ and t → ∞. This yields min

Figure 5

 5 Figure 5.27: Summary of the runs performed at F h = 0.5, in the plane (Ri 0 , Re b = ReF 2 h ). Empty symbols correspond to cases where the Kelvin-Helmholtz does not develop, while filled symbols correspond to unstable runs. The symbols are associated with constant values of the ambient shear flow amplitude: downward triangles correspond to U S = 0.1, open circles to U S = 0.2, diamonds to U S = 0.3, and upward triangles to U S = 0.4.

  Re b , and where Re b = ReF 2

F

  p = -cos( θ)f pc (r) + t sin( θ)f ps (rt → 0, the right-hand side of (4.106) readsF w = -F p + t 2 cos( θ) f 2

  Figure C.1: Streamfunctions ψν 2b (a), ψν 2c (b), and ψν 2s (c) rescaled by r2 as functions of r.

  Route de l'Orme aux Merisiers RD 128 / 91190 Saint Titre : Décorrélation verticale d'un tourbillon soumis à un champ de cisaillement dans un fluide fortement stratifié Mots clés : décorrélation, tourbillon, cisaillement, stratif Résumé : Cette thèse étudie la décorrélation vertical tourbillon initialement rectiligne par un écoulement sinusoïdalement dans un fluide stratifié. Il a été proposé qu'un tel mécanisme devrait déclencher des instabilités de cisaillement et ainsi contribuer à la production de petites échelles en turbulence fortement stratifiée, présente l'atmosphère et les océans dans une gamme d'échelles où la force de Coriolis est négligeable. La première partie étudie les évolutions des énergies et enstrophie totales du tourbillon calculées au moyen de Simulations Numériques Directes (DNS) est différente de celle des écoulements libres non cisaillés du fait de la présence de l'écoulement de saturation entre les termes d'étirement et de dissipation dans le bilan d'enstrophie globale implique que l'enstrophie maximale du tourbillon est nombre de Reynolds à la puissance 2/3 condition de saturation ne rend pas compte de l'effet observé de la stratification. Afin de l'expliquer, la dynamique locale a théoriquement dans une deuxième partie ou grande longueur d'onde. Cette dernière asymptotique fournit les équations d'évolution Title : Vertical decorrelation of a vortex by a strongly stratified shear flow Keywords : decorrelation, vortex, shear, stratification Abstract : This thesis studies the vertical decorrelation of an initially straight vortex by an ambient sinusoidal shear flow in a stratified fluid. It has been conjectured that such process should trigger shear instabilities and, as such, contribute to the generation of small scales in strongly stratified turbulence, which is encountered in the atmosphere and the oceans in an intermediate range of scales where Coriolis effects are negligible. The first part analyses the evolutions of the total energies and enstrophy of the vortex in Direct Numerical Simulations (DNS). The dynamics differs from unsheared flows: because of the presence of the ambient shear flow, the balance between stretching and dissipation terms in the global enstrophy budget maximum enstrophy of the vortex is Reynolds number to the power 2/3 simplified balance does not account for the observed effect of the stratification. Thus, the local dynamics has been investigated theoretically in a second part for short time or long wavelength. The latter asymptotic Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France Décorrélation verticale d'un tourbillon soumis à un champ de cisaillement dans un fluide fortement stratifié tourbillon, cisaillement, stratification Cette thèse étudie la décorrélation verticale d'un par un écoulement cisaillé sinusoïdalement dans un fluide stratifié. Il a été proposé déclencher des instabilités de cisaillement et ainsi contribuer à la production de petites ce fortement stratifiée, présente dans l'atmosphère et les océans dans une gamme d'échelles où partie étudie les évolutions des énergies et enstrophie totales du tourbillon calculées au moyen de umériques Directes (DNS). La dynamique des écoulements libres non cisaillés : t de la présence de l'écoulement ambiant, la condition de saturation entre les termes d'étirement et de dissipation dans le bilan d'enstrophie globale implique que male du tourbillon est proportionnelle au à la puissance 2/3. Néanmoins, cette on ne rend pas compte de l'effet a dynamique locale a été étudiée théoriquement dans une deuxième partie pour temps courts Cette dernière analyse d'évolution de la vitesse angulaire du tourbillon et des On montre principalement que la vitesse angulaire décroît à cause d'effets dynamique et visqueux dynamique est dûe à un resserrement des lignes isopycnes dans le coeur du tourbillon, qui impliq la vorticité verticale afin que la vorticité potentielle se conserve. Dans une troisième partie, des DNS révèlent que l'instabilité de cisaillement apparaît stratification est modérée. L vertical de vitesse horizontale et du gradient vertical de flottabilité sont finement asymptotique lorsque la stratification est forte. de Richardson asymptotique est critique 1/4 nécessaire au déclenchement de l'instabilité de cisaillement. La saturation du cisaillement vertical est dûe au déclin accéléré du tourbillon dans les régions où le cisaillement ambiant est maximal. Ces résultats suggèrent que difficilement déclenchée par des processus de décorrélation dans les écoulements fortement stratifiés, contredisant ainsi les conjectures précédemment formulées dans la littérature.

  Décorrélation verticale d'un tourbillon soumis à un champ de cisaillement dans un fluide fortement stratifié et des déformations de son axe. montre principalement que la vitesse angulaire décroît à d'effets dynamique et visqueux. La décroissance dûe à un resserrement des lignes isopycnes dans le coeur du tourbillon, qui implique une atténuation de la vorticité verticale afin que la vorticité potentielle se Dans une troisième partie, des DNS révèlent que ité de cisaillement apparaît seulement lorsque la Les évolutions du cisaillement vertical de vitesse horizontale et du gradient vertical de finement prédites par l'analyse lorsque la stratification est forte. Le nombre n asymptotique est toujours supérieur au seuil nécessaire au déclenchement de l'instabilité de La saturation du cisaillement vertical est dûe au déclin accéléré du tourbillon dans les régions où le cisaillement ambiant est maximal.

  These equations can be resolved successively. We first imposeψ w0 = -[ψ ps + λ 0 rΩ + (λ 1 -ϑ)r] (B.12)in order to satisfy the same initial condition (4.130) as for the solution computed numerically. Note that ψ ps is defined in (4.120). Then, (B.8) gives ψ w1 = -ψ pc + (λ 1 -ϑ)rΩ +

	Inserting this expansion in (4.106) yields the sequence of equations:
	∇ 2 r ψ w1 + Ω∇ 2 r ψ w0 -	1 r ∂ζ 0 ∂r	ψ w0 = -f pc ,	(B.8)
	2∇ 2 r ψ w2 -Ω∇ 2 r ψ w1 +	1 r ∂ζ 0 ∂r	ψ w1 = -f ps ,	(B.9)
	3∇ 2 r ψ w3 + Ω∇ 2 r ψ w2 -	1 r ∂ζ 0 ∂r	ψ w2 = -	f 2 h 4F 2	,	(B.10)
	4∇ 2 r ψ w4 -Ω∇ 2 r ψ w3 +	1 r ∂ζ 0 ∂r	ψ w3 = -	f 3 12F 2 h	.	(B.11)
					A 2	r ,	(B.13)
	r dS w dt	+ sin( θ) rΩS w (t) +	r dC w dt	d 2 U dz 2
	5 )	d 2 U dz 2 .				(B.7)
			129

cos( θ) rΩC w (t) --sin( θ)ψ w0 (r) + t cos( θ)ψ w1 (r) + t 2 sin( θ)ψ w2 (r) + t 3 cos( θ)ψ w3 (r) + t 4 sin( θ)ψ w4 (r) + O(t

However, another possible choice for the scaling of the slow time is discussed at the end in the light of the final equations.

As for the slow time scale, another possible choice for the scaling of the Reynolds number is discussed at the end of the chapter.

These equations show that the vortex axis evolves because of four different effects in addition to the advection U te x by the shear flow: 4 With the alternative scalings proposed previously, τ = √ εt, 1/ReF 2 h = √ ε/ Re b , these solutions show that δy = O(1) for τ = O(1) as assumed herein, while δx = O( √ ε). However, the fast time would be still present in the wave terms.

Chapter 4

Evolution of a vortex in a strongly stratified shear flow. Part 1.

Asymptotic analyses.

Article in preparation for submission to the Journal of Fluid Mechanics.

49

Abstract

We conduct direct numerical simulations of an initially vertical Lamb-Oseen vortex in an ambient shear flow varying sinusoidally along the vertical in a stratified fluid. The Froude number F h and the Reynolds number Re, based on the circulation Γ and radius a 0 of the vortex, have been varied in the ranges: 0.1 F h 0.5 and 3000 Re 10000. The shear flow amplitude ÛS and vertical wavenumber kz lie in the ranges: 0.02 2πa 0 ÛS /|Γ| 0.4 and 0.1 kz a 0 3π/2. The results are analysed in the light of the asymptotic analyses performed in chapter 4.

The vortex is mostly advected in the direction of the shear flow but also in the perpendicular direction owing to the self-induction. The decay of potential vorticity is strongly enhanced in the regions of high shear. The shear instability is triggered only when the Froude number is moderate F h = 0.5 and for sufficiently high Reynolds number and vertical wavenumber kz .

The initial non-hydrostatic dynamics is well captured by the short-time asymptotic analysis of chapter 4. The long-wavelength analysis predicts also very well the deformations of the vortex axis. The evolutions of the vertical shear of the horizontal velocity of the vortex and of the vertical gradient of the buoyancy at the location of maximum shear are in good agreement with the asymptotic predictions. These quantities saturate owing to both dynamic and viscous effects. The asymptotic Richardson number is shown to depend only on the initial Richardson number Ri 0 and the buoyancy Reynolds number Re b = ReF 2 h . Its minimum value never goes below the critical value 1/4. This proves that the vertical shear cannot grow indefinitely in strongly stratified flows, contradicting the conjecture of [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF].

The asymptotic analysis is not able to predict the development of the shear instability for F h = 0.5 because of the important effect of internal waves, partially neglected in the long-wavelength analysis.

Appendix A

Study of a linear shear flow

This appendix generalizes the short time asymptotic analysis of chapter 4 ( §4.3), to the case of a linear shear flow. The starting point of the present discussion is the equation (4.31), since it is valid for any vortex profile u θ (r) and any shear flow profile U (z).

We consider the linear shear flow given by:

As in §4.3, we look for a solution of (4.31) in the form (4.32). Then, β is found to satisfy: 1 r

A solution to this equation that is not singular at r = 0 and that vanishes at infinity is:

Note that this solution has been derived without specifying the profile u θ (r). As in the case of the sinusoidal shear flow, we have w = w 2 t 2 and b = b 3 t 3 at leading order for t 1. However, in the present case, the expressions of w 2 and b 3 are:

Thus, a transient non-hydrostatic regime also exists at the start-up of the motion when t F h . A striking difference with the sinusoidal shear flow ( §4.3) is that w and b are now proportional to 1/l z whatever its value. As for the second-order horizontal velocities, they read:

These expressions can be retrieved by assuming that the vortex is translated by ∆x(z, t) = U (z)t and ∆y(z, t) = 0. The second-order horizontal velocity field therefore corresponds to a displacement of the vortex as a whole unlike (4.36)-(4.37) for the sinusoidal shear flow.

127

Appendix C

Determination of ψ ν 12

In this appendix, we compute the streamfunction ψ ν 12 forced by the viscous term in (4.154).

By writing it in the form

the viscous part of (4.154) becomes

(C.4)

Like for the determination of the streamfunction ψ * 12 , we assume that the angular velocity follows the law (4.163). Under this hypothesis, ψ ν 2b , ψ ν 2c and ψ ν 2s are rescaled as follows: 133