guided me into exciting and fascinating research topics that I would not have otherwise thought of. Without their constant encouragement and suggestions, I doubt I would ever reach this point. Their passion and enthusiasm for research is infectious and will continue to be an enduring inspiration to me. I am grateful to my thesis referees and examiners, Profs. Joanna M. Nield, Nicolas Taberlet, Mederic Argentina and Sébastien Charnoz, for their enlightening remarks and for their careful reading of the thesis. My thanks also go to Profs. Yuan Wang and Yang Zhang of Department of Fluid Machinery and Engineering at Xi'an Jiaotong University, for their continuous support and encouragement while my studying at ESPCI. I am forever indebted to my family for their most precious gift to me -love. I can never thank them enough for the love, understanding, patience, and encouragement that they unselfishly gave me when I needed them most.

Besides the ones I have mentioned above, there are still many people who accompanied me all along this adventure. I would prefer not to list people's name here, since every list has an end. My sincerely thanks go to all of you. I would like to thank you, for having faith in me, even when I am lacking it myself. I would like to thank you, for encouraging me to try and never to be afraid of failures.

I would like to thank you, for letting me keep in mind that steady, rather than fast pace, wins real games. I would like to thank you, for letting me know that you are with me, regardless of light and dark. I would like to thank you, for letting me know that even though I lose everything, I still have you.

You know who you are . . .

Abstract

This thesis is devoted to the experimental and theoretical investigations of four instabilities associated with the emergence of regular patterns over erodible/flexible substrates, and related to hydrodynamics over a modulated relief.

First, the instability of a flexible sheet clamped at both ends and submitted to a permanent wind is investigated. The flat sheet solution is unstable towards propagative waves, for strong enough wind. We experimentally study the selection of frequency and wavenumber as a function of the wind velocity. These quantities obey simple scaling laws derived from a linear stability analysis of the problem. This phenomenon may be applied for energy harvesting.

Second, an explanation is proposed for the giant ripples observed by spacecraft Rosetta at the surface of the comet 67P. We show that the outgassing flow across a porous surface granular layer and the strong pressure gradient associated with the day-night alternance are responsible for thermal superficial winds. We show that these unexpected patterns are analogous to ripples emerging on granular beds submitted to viscous shear flows. Linear stability analysis of the problem quantitatively predicts the emergence of bedforms at the observed wavelength and their propagation. This description provides a reliable tool to predict the erosion and accretion processes controlling the evolution of small solar system bodies.

Third, we propose a model for rhythmic, dune-like patterns observed on Sputnik Planum of Pluto. Their emergence and evolution are related to the differential condensation/sublimation of nitrogen ice. We show that the temperature and pressure in Pluto's atmosphere are almost homogeneous and steady, and that heat flux from the atmosphere due to convection and turbulent mixing is responsible for the emergence of these sublimation patterns, in contrast to the penitentes instability due to solar radiation.

Last, we report an analytical model for the aeolian ripple instability by considering the resonant grain trajectories over a modulated sand bed, taking the collective effect in the transport layer into account. The model is tested against existing numerical simulations that match experimental observations. viii

French version

Cette thèse porte sur l'étude expérimentale et théorique de quatre instabilités associées à l'émergence de motifs réguliers sur des substrats érodables ou fortement déformables, instabilités liées à l'hydrodynamique sur un relief modulé.

La première partie porte sur l'étude de l'instabilité d'une plaque élastique fixée aux deux bouts et soumise à un écoulement fluide permanent. La solution plane est instable vis-à-vis d'ondes propagatives, lorsque l'écoulement est suffisamment fort. La sélection de fréquence et de longueur d'onde est caractérisée expérimentalement en fonction de la vitesse de l'écoulement. Ces quantités suivent remarquablement les lois d'échelle obtenues par l'analyse de stabilité linéaire du problème. Le principe de l'expérience pourrait être

Table of contents

List of figures xi

List of tables xv

Chapter 1

General introduction

When a fluid interacts with a movable or deformable structure with an internal or surrounding flow, it exerts forces on the structure, hence causing displacement or deformation. The deformation/displacement, however, changes the boundary conditions of the flow at the same time. Problems involving this kind of coupling have attracted a lot of attention in both fundamental science and engineering [START_REF] Padoussis | Fluid-Structure Interactions: Slender Structures And Axial Flow[END_REF]. Besides the traditional applications such as the design of aircraft [START_REF] Anderson | A History of Aerodynamics: And Its Impact on Flying Machines[END_REF], automobiles [START_REF] Katz | Aerodynamics of race cars[END_REF] and bridges [START_REF] Miyata | Historical view of long-span bridge aerodynamics[END_REF], it has also been considered for medical treatments recently, such as in the analysis of aneurysms in large arteries [START_REF] Gerbeau | Fluid-structure interactions in blood flows on geometries based on medical imaging[END_REF] and artificial heart valves [START_REF] Lim | Cavitation phenomenon in mechanical heart valves: The role of squeeze flow velocity and contact area on cavitation initiation between two impinging rods[END_REF]. In addition, the interaction between a fluid and a structure is also relevant in the context of geological fluid mechanics [START_REF] Huppert | The intrusion of fluid mechanics into geology[END_REF]. Amongst other famous examples, let us cite the dynamics of meanders [8] or the spontaneous formation of ripples and dunes [START_REF] Charru | Sand ripples and dunes[END_REF] in the nature. In this thesis, we focus on four subjects involving the interactions between the flow and a flexible or an erodible substrate. They are respectively related to waves on a flexible sheet, to granular patterns formed on an erodible bed and to sublimation patterns.

Flag flapping instability induced by wind

Waves on a flexible sheet are usually termed as the flag flapping instability (Fig. 1.1a). The archetype setup is the case of a cantilevered flexible sheet lying in an axial flow, attached on the up-stream side and freely flapping at the down-stream end. The interactive motion between wind and flags has drawn people's interests and curiosities since long. is the wind that moves!". Listening the back-and-forth argument, Hui-Neng interrupted and told them that they were both wrong, and that instead "It is the mind that actually moves." Although being short from the physical point view , it also diverted the framing of a fascinating scientific question that might be posed as how the passive motions of the flag interacts the unseen wind dynamics. So far, it is already well known that the flag instability results from the competition between the destabilising effect of the pressure applied by the surrounding flowing fluid, which, by virtue of Bernoulli's principle, is lowered above crests and increased in troughs, pushing the sheet away from its flat equilibrium, and the stabilising effect of the bending rigidity of the solid, which tends to restore the sheet flat (Fig. 1.2a). The occurrence of flag flapping is also related to the inertia of the flag, and the system is driven by the vortex shedding at the trailing end (Fig. 1.2b). Compared with other flow-induced instabilities, such as airfoil flutters involving only a limited number of degrees of freedom (usually 1 or 2), large deformations of the flag itself due to the flexibility always induce a large number of degrees of freedom in the motion. The strong coupling between the fluid flow and the elasticity is must therefore be included in the modeling, which leads to significant difficulties to solve the problem. The characteristics of the dynamical system are expected to be influenced by both the flow field and the properties of the flag. Recently, the flapping of flags has been broadly studied due to its wide applications in the paper industry [START_REF] Watanabe | Piezoelectric coupling in energyharvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Watanabe | An experimental study of paper flutter[END_REF] and airfoil flutter [START_REF] Bisplinghoff | Aeroelasticity[END_REF], as well as the biological situations including snoring [START_REF] Huang | Flutter of cantilevered plates in axial flow[END_REF] and the motion of swimming or flying animals [START_REF] Lighthill | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Huber | Siwimming in flasta[END_REF][START_REF] Liao | Fish exploiting vortices decreasing muscle activity[END_REF][START_REF] Muller | Fish'n flag[END_REF][START_REF] Ramanananarivo | Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance[END_REF][START_REF] Shyy | An introduction to flapping wing aerodynamics[END_REF] -see the recent review by Shelley &

Zhang [START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF] and references therein.

Granular patterns on an erodible bed

The emergence of granular patterns have been studied under various situations [START_REF] Aranson | Granular Patterns[END_REF]. A seminal example is that of ripples and dunes, which form spontaneously from an erodible bed sheared by a fluid flow in a wide variety of environments (Fig. 1.3): in water channels, rivers and coastal areas [START_REF] Best | The fluid dynamics of river dunes: A review and some future research directions[END_REF], in deserts [START_REF] Bagnold | The Physics of Blown Sand and Desert Dunes[END_REF] and snow fields [START_REF] Petrich | Snow dunes: A controlling factor of melt pond distribution on arctic sea ice[END_REF] on Earth and under methane or CO 2 atmospheres on other planets [START_REF] Bourke | Extraterrestrial dunes: An introduction to the special issue on planetary dune systems[END_REF][START_REF] Lucas | Growth mechanisms and dune orientation on titan[END_REF][START_REF] Charnay | Methane storms as a driver of titan's dune orientation[END_REF], in hydraulic engineering and industrial pipe flows [START_REF] Stevenson | The transport of particles at low loading in near-horizontal pipes by intermittent flow[END_REF][START_REF] Schaflinger | An experimental study of viscous resuspension in a pressure-driven plane channel flow[END_REF]. Their size can range from the centimeter scale for subaqueous ripples to one hectometer for large river mega-dunes, from one decimeter for the smallest aeolian dunes to one kilometer for the largest ones. Mature, finite-height bedforms are typically asymmetric, with an avalanche slip face on their lee side. More generally, bedforms exhibit different shapes depending on the symmetries of the fluid forcing, or the boundary conditions [START_REF] Andreotti | Giant aeolian dune size determined by the averaged depth of the atmospheric boundary layer[END_REF]. Since the pioneering work by Bagnold [START_REF] Bagnold | The Physics of Blown Sand and Desert Dunes[END_REF], a number of studies have been stimulated, and significant progresses have been achieved in the understanding of the formation of ripples and dunes [START_REF] Charru | Sand ripples and dunes[END_REF][START_REF] Engelund | Sediment ripples and dunes[END_REF], of coastal forms [START_REF] Blondeaux | Mechanics of coastal forms[END_REF], of fluvial sedimentary patterns [START_REF] Seminara | Fluvial sedimentary patterns[END_REF], and of snow bedforms [START_REF] Filhol | Snow bedforms: A review, new data, and a formation model[END_REF]. The dynamics of these patterns results from the interaction between the fluid flow and the bed topography through sediment transport. In the following, we recall a brief description of particle transport, which will be used in the discussions afterwards.

Sediment transport mode

The particle at the bed surface are set into motion when the hydrodynamical drag is strong enough to overcome a certain threshold related to the bed disorder. The ability of the fluid to put the grains of the bed into motion is quantified by the Shields number defined as

Θ = τ (ρ p -ρ)gd , ( 1.1)
with g the gravity acceleration, d the grain size, ρ and ρ p the fluid and particle densities;

and τ = ρu 2 * the shear stress exerted by the fluid on the bed, where u * is the shear velocity. The threshold Shields number Θ t is closely related to the threshold velocity u t , which is discussed in details in Section 5.2.

Observations have shown that the particles are transported in different modes as schematically shown in Fig. 1.4, associated with the forces acting on them. In the aeolian case, the ratio of the particle density to the fluid density is large so that the grains are mainly transported in saltation, in a succession of jumps. When the impact of saltating grains on the bed is strong enough, they release a splash-like shower of ejected grains that make small hops, and this secondary transport mode is called reptation. In the subaqueous case, the grains and the fluid have comparable densities. The transport is mainly a turbulent suspension when the velocity of turbulent fluctuations is larger than the settling velocity.

When gravity is large enough to confine sediment transport in a layer at the surface of the bed, one refers to bedload: the grains are either hopping in saltation or roll and slide at the bed surface, with long-contacts between the grains, termed as traction.

Aeolian

Saturated flux

Our primary interest is the evolution of the bed surface due to the sediment transport, which is related to the particle flux q, in m 2 /s, counting the volume of the grains (packed at the bed volume fraction) passing a vertical surface of unit width, and per unit time. To understand, let's consider a simplest case of an infinite flat sediment bed submitted to a steady flow. After a transient time, a dynamical equilibrium between flow and particle transport is reached, which means that there are as many particles deposited as particles eroded from the bed. The erosion and deposition fluxes balance each other; i.e. there is no net erosion nor deposition of particles. The dynamical equilibrium is characterised by the so-called saturated flux q sat which determines the resulting sediment transport flux, and q sat is spatially homogeneous. It is obvious that the stronger the flow, the more it can carry grains: q sat is therefore an increasing function of the increasing wind strength, which vanishes below a threshold shear velocity u * = u t . As shown in Fig. 1.5a, the behavior of q sat for bedload has been measured in hydraulic channels, and the data are traditionally fitted by an empirical formula by Meyer-Peter & Müller (1948) [START_REF] Meyer-Peter | Formulas for bed load transport[END_REF]:

q sat ∝ ρ p -ρ ρ p gd 3 (Θ -Θ t ) 3/2 . (1.2)
This scaling law is well understood from Bagnold's original idea [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF]. The saturated flux can be decomposed as the product of the number N of transported grains per unit area by the mean grain horizontal velocity u p : q sat = πd 3 /(6φ b )N u p , with φ b the bed volume fraction. The total shear stress ρ f u 2 * can be decomposed as the sum of the fluid-borne shear stress τ f and the grain-borne shear stress τ p . Now τ p is proportional to the moving grain density N and to the drag force acting on a grain moving at the average velocity u p due to a flow at the velocity u. If the grains are in a steady motion, the drag force balances a resistive force due granular friction, collisions with the bed, etc. These different dissipative mechanisms can be modeled as an overall effective friction force characterized by a friction coefficient. Assuming that the transported grains do not disturb the flow, the flow velocity u around grains must be proportional to the shear velocity, so that u p scales linearly with u * (or Θ, Eq. 1.1). Furthermore, saturation is reached when the fluid-borne shear stress equals the transport threshold at the surface of the static bed. As a consequence, the number of transported particles per unit area is proportional to the excess shear stress, or Θ -Θ th . Taking the product of N and u p , it then gives a transport law as in Eq. 1.2. We will derive the transport law based on Bagnold's idea for the cometary case in Chapter 5.

Saturation length

The sediment flux is no longer uniform when the bed topography is modulated by the patterns (bedforms). There is ample experimental evidence that that transport does not adapt instantaneously to a spatial change of the shear stress but shows a relaxation in space and time [START_REF] Anderson | Simulation of aeolian saltation[END_REF]. Considering that the transport over bedforms is never far from its saturated state, the problem can therefore be described by a linear relaxation towards the saturated state, and one obtains: where L sat and T sat are called the saturation length and times [START_REF] Andreotti | Selection of dune shapes and velocities. part 1: Dynamics of sand, wind and barchans[END_REF][START_REF] Sauermann | Continuum saltation model for sand dunes[END_REF][START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF][START_REF] Claudin | Transport relaxation time and length scales in turbulent suspensions[END_REF][START_REF] Narteau | Setting the length and time scales of a cellular automaton dune model from the analysis of superimposed bedforms[END_REF][START_REF] Parker | Sediment inertia as cause of river anti-dunes[END_REF][START_REF] Andreotti | A two species model of aeolian sand transport[END_REF][START_REF] Anderson | Wind modification and bed response during saltation of sand in air[END_REF][START_REF] Valance | Ripple formation over a sand bed submitted to a laminar shear flow[END_REF]. Comparing to the bedform growth time (∼ 10 2 s for subaqueous ripples and ∼ 10 5 s for aeolian dunes), T sat is much smaller for the ripples and dunes, which justifies the simplifying assumption that the fluid flow can be computed as if the bed were fixed. To understand L sat , we consider a homogeneous flow on a granular bed, which extends only in the half space x > 0. Upstream of this bed (x < 0) is not erodible but has the same hydrodynamic roughness. The flux q at the upwind boundary (x = 0) is zero; it grows and then relaxes exponentially towards q sat over relaxation length L sat , which thus characterise the length scale over which sediment transport relaxes towards equilibrium. When the flux is locally larger than the saturated flux, it decreases over space so that grains are deposited (see Fig. 1.5b). If the flux is locally smaller than the saturated flux, it grows over space, causing erosion of the bed. L sat has been theoretically and experimentally studied for saltation [START_REF] Andreotti | A two species model of aeolian sand transport[END_REF][START_REF] Andreotti | Measurements of the aeolian sand transport saturation length[END_REF][START_REF] Hersen | Relevant lengthscale of barchan dunes[END_REF][START_REF] Durán | On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws[END_REF][START_REF] Iversen | The effect of wind speed and bed slope on sand transport[END_REF]]

T sat ∂ q ∂ t + L sat ∂ q ∂ x = q sat -q, ( 1.3)
L sat ∝ ρ p ρ d, (1.4)
and suspension [START_REF] Van Rijn | Sedimentation of dredged channels by currents and waves[END_REF][START_REF] Ashida | On the calculation method of the concentration of suspended sediment under non-equilibrium condition[END_REF][START_REF] Claudin | Transport relaxation time and length scales in turbulent suspensions[END_REF][START_REF] Jobson | Vertical transfer in open channel flow[END_REF]]

L sat ∝ u * H V fall , (1.5)
with H the fluid depth and V fall the falling velocity.

Contrarily to the aeolian transport, L sat has never been directly measured for bedload, mainly because it is much smaller, of the order of few grain sizes. Indirect estimates by Fourrière et al. give very small values scaling with the grain diameter, typically 10-20d [START_REF] Fourrière | Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by non-linear pattern coarsening[END_REF]. We will start over the discussion of bed load saturation length in Chapter 5.

Stability analysis of a flat granular bed

The stability analysis gives the time and length scales at which bedforms emerge from a flat bed. Studies have shown that aeolian dunes and subaqueous ripples form by the very same linear instability, which has been modeled and quantitatively tested against laboratory measurements [START_REF] Charru | Sand ripples and dunes[END_REF]. The destabilizing effect results from the phase advance of the wind velocity just above the surface with respect to the elevation profile. As in Fig.

1.6a, it shows a linear response of the basal shear stress τ, and we note, especially, the phase advance with respect to the bottom. This means that the maximum shear stress locates upstream of the crest. The stabilizing mechanism comes from the space lag between sediment transport and wind velocity, which is characterized by the saturation length L sat [START_REF] Charru | Sand ripples and dunes[END_REF][START_REF] Sauermann | Continuum saltation model for sand dunes[END_REF][START_REF] Andreotti | Measurements of the aeolian sand transport saturation length[END_REF] (Fig. 1.6b). Grains are eroded (deposited) when the flux increases (decreases) in space. Instability takes place when the crest is in the deposition zone, i.e. when the maximum of the sediment flux is upwind of the crest. The dispersion relation for the linear stability analysis of a flat granular bed has been studied by Charru et al [START_REF] Charru | Sand ripples and dunes[END_REF]. The growth rate σ and propagation velocity c of a bed modulation of the wavenumber k = 2π/λ, where λ is the wavelength, is given by:

σ = k 2 ( -) -kL sat 1 + (kL sat ) 2 , (1.6) c = k + ( -)kL sat 1 + (kL sat ) 2 . ( 1.7)
In these expressions, ≡ τ∂ τ q sat quantifies the sediment transport, and encodes the fact that the threshold for transport is sensitive to the bed slope with = 1 µ τ t /τ, where µ = tan(29 • ) ≃ 0.55 is the tangent of the avalanche angle. and are the components of the basal shear stress in phase and in quadrature with the bottom, respectively, whose determination are discussed in Section 6.4.

Regarding aeolian ripples, it has been long known that they do not share the same 'linear instability' formation mechanism as dunes and subaqueous ripples [START_REF] Wilson | Aeolian bedforms: Their development and origins[END_REF]. Until recently, their emergence was ascribed to a geometrical effect responsible for the modulation of sediment transport. We will revisit this problem in Part IV.

Crest

Maximum wind velocity

Sublimation patterns on an ice bed

Different patterns on a sublimating surface, such as penitentes and ice waves, have been observed on Earth [START_REF] Nichols | Nieves penitentes near boston, massachusetts[END_REF][START_REF] Frezzotti | Snow megadunes in antarctica: Sedimentary structure and genesis[END_REF], as in Fig. 1.7. Some of them also form (or are expected to form) on other planetary bodies [START_REF] Cutts | Evolution of martian polar landscapes: interplay of long-term variation in perennial ice cover and dust storm intensity[END_REF][START_REF] Howard | The role of aeolian processes in forming surface features of the martian polar layered deposits[END_REF][START_REF] Howard | How rough is the surface of europa at lander scale?[END_REF]. These patterns are always found in specific environments, where the partial pressure of vapour (corresponding to the ice substance) in the atmosphere is low, sublimation therefore significantly contributes to ablation and takes a part in the development of patterns on the ice surface.

It is already known that penitentes form in sublimation conditions by differential ablation due to self-illumination, vapor diffusion and heat conduction [START_REF] Betterton | Theory of structure formation in snowfields motivated by penitentes, suncups, and dirt cones[END_REF][START_REF] Claudin | Physical processes causing the formation of penitentes[END_REF]. Regarding ice waves, the emergence and evolution can be described using a similar hydrodynamics as granular (snow or sand) dunes and ripples. They differ by the presence of particle transport in one case and the sublimation at the interface in the other. Recent studies have shows that complex interactions between sublimation-related mass transfer and turbulent flow in the lower-atmosphere lead to the development of stunning spiral-shaped topographic ice waves at the surface of the North Polar cap of Mars, and it is the periodic spatial variations in sublimation-related ablation rates that are responsible for the development of these topographic waves [START_REF] Herny | Sedimentation waves on the martian north polar cap: Analogy with megadunes in antarctica[END_REF][START_REF] Smith | The spiral troughs of mars as cyclic steps[END_REF]. So far, the genesis of these patterns is still not completely known, and further work is needed to understand the role of sublimation in the development of such wavy patterns, the relevant dynamical mechanisms controlling the mass balance, the dynamics at different scales and so on.

Another group of similar phenomena are dissolution patterns, such as scallops and icicle ripples, which always develops in limestone caves and in caves in ice, as well as other precipitation and dissolution interfaces [START_REF] Ogawa | Surface instability of icicles[END_REF][START_REF] Thomas | Size of scallops and ripples formed by flowing water[END_REF][START_REF] Meakin | Geological pattern formation by growth and dissolution in aqueous systems[END_REF][START_REF] Blumberg | Experimental and theoretical studies of dissolution roughness[END_REF]. These patterns result from the interaction of a soluble surface and an adjacent turbulent flow, and this is beyond the discussion in this thesis.

Fluid flow over the rippled patterns

All the patterns mentioned above are resulted from the interaction between substrates and fluid flows over the interfaces. A good understanding of the flow field near the patterns plays a key role in explaining the emergence and evolution of these patterns. As a general description, we briefly introduce here a two-dimensional incompressible flow over the interface, with x, z and ξ denote the flow direction, vertical direction and the interface profile, respectively. Following the standard separation between average quantities and fluctuating ones (denoted with a prime), the mean velocity field u i is governed by Reynolds

∂ i u i = 0, (1.8) ρ(∂ t u i + u j ∂ j u i ) = ∂ j τ i j -∂ i p, ( 1.9)
where τ i j contains the deviatoric part of the Reynolds stress tensor -ρu ′ i u ′ j . The trace of the Reynolds stress tensor is included inside p. We use here a Prandtl-like first order turbulence closure in which the distance to the bed zξ determines the mixing length ℓ, and the mixing frequency is given by the strain rate modulus |γ| = 1 2 γi j γi j , where we have introduced the strain rate tensor γi j = ∂ i u j + ∂ j u i .

In the general case, we can write the stress tensor components as the sum of the viscous and turbulent contributions:

τ xz = ρ ℓ 2 |γ| + ν γxz ,
(1.10)

τ x x = ρ ℓ 2 |γ| + ν γx x , (1.11
)

τ zz = ρ ℓ 2 |γ| + ν γzz , (1.12)
where ν is the gas kinematic viscosity. In these expressions, the strain tensor components are given by

γxz = γzx = ∂ z u x + ∂ x u z , γx x = 2∂ x u x and γzz = 2∂ z u z = -γ x x , (1.13)
and the strain modulus by:

|γ| 2 = 2(∂ x u x ) 2 + 2(∂ z u z ) 2 + (∂ z u x + ∂ x u z ) 2 = 4(∂ x u x ) 2 + (∂ z u x + ∂ x u z ) 2 . (1.14)
For a turbulent boundary layer, the mixing length cannot exceed the distance to the wall, and it is generally estimated by ℓ = κ(zξ), with κ = 0.4 the von Kármán constant.

This expression describes well the fully turbulent part of the boundary layer, yet excluding the viscous sublayer and buffer layer close to the wall. To apply the mixing length model in the entire turbulent boundary layer, various empirical expression have been proposed, and we employ here the van Driest formula by introducing a damping function [START_REF] Van Drist | On turbulent flow near a wall[END_REF],

ℓ = κ(z -ξ) 1 -exp - (τ xz /ρ) 1/2 (z -ξ) νR t , (1.15)
where R t is the transitional Reynolds number. Following Hanratty [START_REF] Charru | Sand ripples and dunes[END_REF][START_REF] Abrams | Relaxation effects observed for turbulent flow over a wavy surface[END_REF][START_REF] Frederick | Velocity measurements for a turbulent nonseparated flow over solid waves[END_REF], we consider that the transitional Reynolds number depends on a dimensionless number which depends on, but lags behind the pressure gradient:

a ν u * ∂ x = ν ρu 3 * ∂ x (τ x x -p) - (1.16)
where a is the multiplicative factor in front of the space lag and u * is the shear velocity.

is positive when the pressure decreases in space, which is stabilizing. The transitional Reynolds number is therefore an increasing function of . We also introduce b as the relative variation of R t due to the pressure gradient:

b = 1 R 0 t dR t d > 0. (1.17)
where R 0 t = 25 is the transitional Reynolds number for the homogeneous case. The generic hydrodynamical description above will be used simplify or further developed in the coming chapters for different situations.

Outline of the thesis

The present thesis is elaborated in four parts:

1. In Part I, we study a situation related to but different from the flag flapping: a flexible or compliant material clamped at both ends, which develops traveling waves when submitted to a flow. We first measure the frequency and the wavenumber of the waves varying the wind velocity and for different materials (paper and plastic sheets), and then we carry out a linear stability analysis based on a simple hydrodynamic assumption and on the Euler-Bernoulli beam theory for the sheet to study the temporal and spatial scalings, which are then compared with the experimental data. 3. In Part III, we propose an explanation for the kilometer scale patterns on the surface of Pluto, revealed by the New Horizons mission of NASA. Firstly, we study Pluto's atmosphere in Chapter 8, and then we model the ice waves in a general case by 1.5 Outline of the thesis linking the emergence and evolution to the instability of an interface between a sublimating ice bed and its turbulent vapor flow in Chapter 9. We test the model against the observed patterns in the end. [START_REF] Miyata | Historical view of long-span bridge aerodynamics[END_REF]. In Part IV, we report an analytical model based on the sand grain trajectories to study the emergence and evolution of the aeolian ripples. The model is then compared to the numerical data; and we also present evidence from wind tunnel experiments.

Introduction

When it comes to the flag flapping instability, the earliest explanation dates back to Rayleigh [START_REF] Rayleigh | On the instability of jets[END_REF], who found that a flag of infinite span and infinite length is always unstable. When considering a flag of finite dimensions, this problem becomes more difficult and depends on its aspect ratio, defined as the ratio of the flag span to the flag length.

The slender body approach and the airfoil theory have been respectively employed to implement the calculations for asymptotically small [START_REF] Lighthill | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Lemaitre | Instability of a long ribbon hanging in axial air flow[END_REF] and large [START_REF] Kornecki | On the aeroelastic instability of twodimensional panels in uniform incompressible flow[END_REF][START_REF] Huang | Flutter of cantilevered plates in axial flow[END_REF][START_REF] Watanabe | An experimental study of paper flutter[END_REF][START_REF] Guo | Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow[END_REF] aspect ratios. A recent unified model by Eloy et al. [START_REF] Eloy | Flutter of a rectangular plate[END_REF][START_REF] Eloy | Aeroelastic instability of cantilevered flexible plates in uniform flow[END_REF] has considered the intermediate case. Experimental studies have been carried out in wind tunnels [START_REF] Taneda | Waving motions of flags[END_REF][START_REF] Huang | Flutter of cantilevered plates in axial flow[END_REF][START_REF] Watanabe | An experimental study of paper flutter[END_REF], in water flumes [START_REF] Shelley | Heavy flags undergo spontaneous oscillations in flowing water[END_REF] and even in flowing soap films [START_REF] Zhang | Flexible laments in a owing soap lm as a model for one-dimensional ags in a two-dimensional wind[END_REF]. The first experiments of Taneda were performed in a vertical wind tunnel with flags of different materials (silk, flannel, canvas, muslin) and shapes (triangles, rectangles), and it is reported that the flags do not flap in slow flows due to the stabilizing effects of both viscosity and gravity [START_REF] Taneda | Waving motions of flags[END_REF]. Later, Datta & Gottenberg conducted experiments with long ribbons hanging vertically in downward flows, and the critical flow velocity for the onset of flapping was studied as a function of the length, width and thickness of the ribbons [START_REF] Datta | Instability of an elastic strip hanging in an airstream[END_REF]. This small aspect ratio regime has recently been experimentally revisited by Lemaitre [START_REF] Lemaitre | Instability of a long ribbon hanging in axial air flow[END_REF]. Experiments for larger and intermediate aspect ratios have also been reported by Huang [START_REF] Huang | Flutter of cantilevered plates in axial flow[END_REF], Watanabe et al [START_REF] Watanabe | An experimental study of paper flutter[END_REF], Yamaguchi et al [START_REF] Yamaguchi | Flutter limits and behavior of a flexible thin sheet in high-speed flow -ii: Experimental results and predicted behaviors for low mass ratios[END_REF] and Eloy et al [START_REF] Eloy | Aeroelastic instability of cantilevered flexible plates in uniform flow[END_REF]. In these experiments, the critical velocity, at which the flag flaps spontaneously with a large amplitude and a well-defined frequency, is systematically found higher than theoretical predictions [START_REF] Argentina | Fluid-flow-induced flutter of a flag[END_REF] and the origin of this discrepancy has recently been related to inherent planarity defects [START_REF] Eloy | The origin of hysteresis in the flag instability[END_REF]. In a more recent experiment, Kim et al.

have investigated the occurrence of the flapping of an inverted flag, with a free leading edge and a fixed trailing edge [START_REF] Kim | Flapping dynamics of an inverted flag[END_REF].

Various numerical approaches have also been used to tackle the different aspects of this problem. The experiments of Zhang et al [START_REF] Zhang | Flexible laments in a owing soap lm as a model for one-dimensional ags in a two-dimensional wind[END_REF]. were reproduced by direct numerical simulation using immersed boundary method (IBM) [START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF] and arbitrary Lagrangian Eulerian (ALE) [START_REF] Sawada | Fluid-structure interaction analysis of a two-dimensional flag-in-wind problem by the ale finite element method[END_REF][START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF]. Both numerical methods have reproduced well qualitatively the stretchedstraight state and bistability-switching. Alben et al employed a flexible body vortex sheet model to compute the fluid forcing and the flow field around a flapping 2D flag [START_REF] Alben | Flapping states of a flag in an inviscid fluid: Bistability and the transition to chaos[END_REF], and they reported a chaotic state, characterised by undefined amplitude and frequency, may appear when the incoming flow velocity is much larger than the critical velocity. Michelin et al also observed this chaotic state by using a unsteady point vortex model [START_REF] Michelin | Vortex shedding model of a flapping flag[END_REF].

In this chapter, we report experiments of wind-generated waves on a flexible sheet clamped at both ends, thus avoiding the flapping phenomenon and subsequent vortex shedding from the trailing edge. These elastic surface waves are induced by the same pressure-related instability mechanism as in the case of the flag flapping [START_REF] Kornecki | On the aeroelastic instability of twodimensional panels in uniform incompressible flow[END_REF][START_REF] Hansen | An experimental study of turbulent flows over compliant surfaces[END_REF][START_REF] Hansen | An experimental study of flow-generated waves on a flexible surface[END_REF][START_REF] Kim | Space-time characteristics of a compliant wall in a turbulent channel flow[END_REF], but their dynamics and length/time scale selection are different, and were not investigated before. The content is organized as follows. The experimental setup and the data processing techniques are described in the following. Then, the linear instability analysis is carried out. Finally, experimental results are compared with the theoretical predictions.

Experimental study

Experimental setup

As schematically illustrated in Fig. 2.1, the experiments are conducted in a wind tunnel with a square cross section. We denote as x, y and z the longitudinal, transverse vertical and horizontal axes, respectively. The wind flow is induced along x by imposing the pressure at the inlet. The wind velocity V is monitored with an anemometer (Testo 405-V1, Appendix A) at a fixed position at the exit of the tunnel. A flexible sheet of width ≃ 4 cm and of length L is placed at the centre of the tunnel and clamped at both ends on fixed masts. We denote as L 0 the distance between the two masts and define ∆L = L -L 0 (different values of ∆L have been used, see below). The reference coordinate x = 0 is chosen at the inlet mast. The air flow is uniformly injected on both faces of the sheet, in order to avoid the formation of vortices. A transverse orientation of the sheet along y and z have both been tested, showing qualitatively similar behaviors (quantitative differences exist, due to gravity, see Section 2.4), and in what follows all data correspond to a vertical orientation. The experiments are conducted either with paper or plastic (bi-oriented polypropylene) sheets of different thicknesses. The values of the relevant physical characteristics of these materials

Experimental data

The wind flow generates waves on the sheet, that propagate downstream. The highfrequency movies allow us to follow in detail the kinematics of these waves. In the case of the paper sheet, they are mostly transverse to the wind, i.e. the sheet is not twisted.

However, due to a smaller flexibility, the plastic sheets show three-dimensional motions for the stronger winds. For a given wind velocity, the experimental data are obtained as follows. The borders of the sheet are detected from the image sequences by maximizing the correlation with an analyzing wavelet (Fig. 2.2a). The average profile ζ(x, t) representing the two-dimensional shape of the sheet is determined with a sub-millimetric resolution. The typical time variation of ζ at a fixed value of x is displayed in Fig. 2.2b, and shows harmonic oscillations: for a given x, ζ(x, t) is well represented by the function Acos(ωtφ), where A is an amplitude; ω is an angular frequency and φ is a phase. All these three quantities a priori depend on x (and V ). However, as shown in Fig. 2.2c, ω turns out to be constant all along the sheet and the phase is linearly related to space φ = kx (Fig. 2.2d), corresponding to a constant wave propagation velocity c = ω/k. k is the wavenumber of the waves, and λ = 2π/k is the wavelength. The behavior of A(x) shows several regimes (Fig. 2.2e). A vanishes at both ends of the sheet, as it should because of the clamping. In between, it first increases rapidly, then slowly decreases in a more noisy way over most of the tunnel, and finally quickly drops at the very end. The first part can be associated with the spatial development of the instability. As shown in section 2.4, the typical value of the amplitude in the second regime is dictated by the geometrical constraint that relates A to k and ∆L/L 0 . We have conducted experiments similar to that corresponding to Fig. 2.2, systematically varying the air flow velocity and the total length of the sheet L, for all three types of sheet. As in Fig. 2.3, the dependence of the angular frequency and wave number on wind velocity is plotted, and some simple monotonical variations are observed as suggested by the phenomenological fitting, for which a rigorous explanation is needed. We now mostly focus on the frequency and the wavenumber of the waves, for which a comparison with an analytical theory is possible.

Theoretical modelling

The purpose of this section is to provide a brief but self-contained summary of the theoretical framework within which we analyze our experimental data. The theoretical description of the flow over a flexible sheet has been treated in a general way several decades ago, as e.g.

summarized by Paidoussis [START_REF] Padoussis | Fluid-Structure Interactions: Slender Structures And Axial Flow[END_REF], see chapter 10 and references therein. Here, we restrict this analysis to a two-dimensional linear perturbation theory. Furthermore, we hypothesize that the air flow can be decomposed into a turbulent inner boundary layer and an outer laminar flow which can be described as an incompressible perfect fluid. As we only need the pressure field, which is almost constant across the inner layer according to the boundary layer theory, we will simply describe the outer layer. Under these simplifying assumptions, we are able to derive analytical scaling laws for the frequency and the wavenumber of the most unstable mode in the asymptotic limit of either very flexible or very rigid sheets, and which were not available in the literature.

Governing equations

We consider a flexible sheet of infinite span and length submitted to an air flow along the x-axis. For later rescalings, we denote by V the characteristic velocity of the wind, i.e.

the average air velocity at a given and fixed altitude z w (in the experiment z w ≃ 7.5 cm).

Assuming that the motion of the sheet is independent of the coordinate y, we denote ζ(x, t)

as its deflection with respect to the reference line z = 0. In the limit of small deflections with respect to a flat reference state, the sheet obeys the linearized Euler-Bernoulli beam equation:

m ∂ 2 ζ ∂ t 2 + D ∂ 4 ζ ∂ x 4 + δp = 0, (2.1)
with δp the air pressure jump across the sheet. m and D are respectively the mass per unit surface and the bending rigidity (Table 2.1).

Neglecting viscous stress in the outer layer and assuming incompressibility (recall that velocities are on the order of a few m/s, i.e. corresponding to very low Mach numbers), mass and momentum conservations for the flow field are therefore expressed by Euler equations:

∇ • u = 0, (2.2) ∂ u ∂ t + (u • ∇) u = - 1 ρ ∇p, (2.3)
where u and p are the velocity and pressure fields, and ρ is the air density. Finally, the fluid velocity on the sheet should be equal to the sheet velocity, in order to ensure the impermeability of the sheet:

u(z = ζ) • n = dζ dt , (2.4)
where n is the unit vector normal to the sheet. Equations (2.2-2.4) must be consistently linearized in the limit of small sheet deflections, and together with (2.1) they form a closed set that we analyze in the next sub-section.

Linearized problem

The base state the perturbations

We consider that the sheet is long enough to ignore the influence of boundaries and the base state corresponds to a flat sheet, ζ (0) = 0. And it reads, for flow:

u (0) x = V ; u (0) z = 0; p (0) = P 0 , (2.5)
where P 0 is the reference pressure. This state satisfies the conservation equations and the boundary conditions above. Perturbing this base state, the velocity, pressure, and sheet position are written as (1) x ; u z = u (1) z ; p = P 0 + p (1) ; ζ = ζ (1) .

u x = V + u
(2.6)

Considering the boundary conditions, one obtains for flow u (1) x = 0; u (1) z = 0; p (1) = 0 for z → ∞;

(2.7) And Eq. 2.4 gives the normal velocity continuity

∂ ζ (1) ∂ t + V ∂ ζ (1) ∂ x = u (1) z | z=0 . (2.8)

Linearized equations and normal modes

Introducing Eq. 2.6 into the governing equations (2.2 and 2.3), one obtains the linearized equations with respect to the base state (2.5) by neglecting the products of the perturbations, ∂ u (1) x

∂ x + ∂ u (1) z ∂ z = 0, (2.9)
∂ u (1) x ∂ t + V ∂ u (1) x

∂ x = - 1 ρ ∂ p (1) ∂ x , (2.10) ∂ u (1) z ∂ t + V ∂ u (1) z ∂ x = - 1 ρ ∂ p (1) ∂ z (2.11)
we obtain a system with constant coefficients independent of x and t with Eq. 2.8, 2.9, 2.10 and 2.11. The dependence of the solutions on x and t is therefore exponential, and the perturbations can be sought in the form of normal modes, characterized by the complex frequency Ω:

ζ (1) = Ae ikx-iΩt , (2.12)
u (1) x = V kAe -iΩt+ikx-κz , (2.13)

u (1) z = V kAe -iΩt+ikx-κz , (2.14)
p (1) = ρV 2 kAe -iΩt+ikx-κz , (2.15)

Ω will be later decomposed into real and imaginary parts as Ω = ω + iσ, where σ is the temporal growth rate of the perturbation. κ is the spatial decay rate along the z-axis. A sets the amplitude of the sheet deflection, and, as it should in the framework of a linear analysis, it will factorize out of all results when solving the governing equations in the asymptotic limit kA ≪ 1. V sets the dimensionful reference for the velocities, and ρV 2 does that for the pressure. The unknowns are thus the dimensionless quantities , and as well as κ, and they are to be determined by the above governing equations.

Plugging (2.12-2.15) into (2.8-2.11), and treating the first order in kA, we first find

κ 2 = k 2 , ( 2.16)
from which we set κ = k (resp. κ = -k) for the fields in the region z > 0 (resp. z < 0), in order for the perturbation to decay away from the sheet. The three other unknowns are found as:

= ± kV -Ω kV , (2.17)
= i kV -Ω kV , (2.18) = ∓ (kV -Ω) 2 (kV ) 2 , ( 2.19)
where the ± sign corresponds to the positive/negative region (in z). For our purpose, the most important quantity is δp the pressure jump across the sheet, which writes:

δp = p(z → 0 + ) -p(z → 0 -) = -2ρ (kV -Ω) 2 k ζ. (2.20)
This allows us to express the dispersion relation from Eq. 2.1 as:

mΩ 2 -Dk 4 + 2ρ (kV -Ω) 2 k = 0, (2.21)
whose properties are studied in the next sub-section. Note again that this equation is a simpler form of the dispersion relation derived in [START_REF] Padoussis | Fluid-Structure Interactions: Slender Structures And Axial Flow[END_REF] (chapter 10) where the distance to the bottom wall has been sent to infinity and the viscous drag, the spring stiffness of the elastic foundation and the plate tension have been set to zero.

Dispersion relation

The above equation can be made dimensionless by setting Ω = mΩ ρV for the frequency, k = mk ρ for the wavenumber, and D = ρ 2 D m 3 V 2 for the bending rigidity. With these rescaled variables, Eq. 2.21 writes

Ω2 -Dk 4 + 2 k ( k -Ω) 2 = 0, (2.22)
and D is its only parameter. Because it is quadratic in Ω, (2.22) can easily be solved as: Paper waves in the wind Above kc , Eq. 2.22 has two real roots, and the complex frequency Ω of any perturbation is real (the expression below the square root is positive), corresponding to a wave which propagates without growth nor decay (σ = 0). On the other hand, any perturbation with a wavenumber below kc has an angular frequency and a growth rate given by: We display σ (positive branch) and ω as functions of k in Fig. 2.5. Both km and kc depend on the rigidity of the sheet, and they are larger for smaller D. In the unstable range of k, ω is independent of D. It increases linearly with k, and eventually saturates to the value 2 when k reaches values on the order of unity. Interestingly, this saturated regime corresponds to waves with vanishing phase (ω/k) and group (dω/dk) velocities. Beyond kc , σ = 0 and ω enters another regime (not shown in Fig. 2.5b) where it asymptotically varies like the square of the wavenumber.

Ω = 2 k k + 2 ± -2 k3 + 2 Dk5 + Dk6 ( k + 2) 2 . ( 2
ω = 2 k k + 2 , (2.25) σ = ± 2 k3 -2 Dk5 -Dk6 ( k + 2) 2 . ( 2
for k < kc , ω is independent of D.
In this temporal stability analysis, both phase and group velocities are found positive, which may suggest a convective instability, as often the case for instabilities generating propagative waves. However, performing the spatial stability analysis of these equations, one finds unstable modes with both positive and negative group velocities, suggesting an absolute instability. Previous theoretical analyses of this issue [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF][START_REF] De Langre | Absolute unstable waves in inviscid hydroelastic systems[END_REF] (and references therein) have shown that the instability is absolute below a threshold in D around 10 2 , i.e. at small enough D or large enough V . This therefore justifies the temporal analysis performed here.

Asymptotic analysis and scaling laws

Scaling laws for the characteristics of the most unstable mode ( km , ωm , σm ) as well as for the cut-off wavenumber kc can be analytically derived in the limits of asymptotically small and large D. kc and km are calculated from (2.24) and (2.27), respectively. ωm and σm are obtained by introducing km into Eqs. 2.25,2.26. Expanding these equations in the limit Paper waves in the wind 10 -14 10 -14 10 -10 10 -6 10 -10 -2 10 2 10 10 6 10 10 10 10 14 10 -14 10 -14 10 -10 10 -6 10 -10 -2 10 2 10 10 6 10 10 10 10 14 10 -14 10 -14 10 -10 10 -6 10 -10 -2 10 2 10 10 6 10 10 10 10 14 10 -14 10 -14 10 -10 10 -6 10 -10 -2 10 2 10 10 6 10 10 10 10 14 D → 0 and D → ∞, one obtains for the wavenumbers:

kc ∼ 2 D 1/3 , km ∼ 1 2 D 1/3 for D → 0, (2.28) kc ∼ 1 D 1/2 , km ∼ 2 5 D 1/2 for D → ∞. (2.29)
Similarly, the angular frequencies and growth rates scale as:

ωm ∼ 2, σm ∼ 27 16 D 1/6 for D → 0, (2.30) ωm ∼ 3 5 D 1/2 , σm ∼ 2 5 1/2 3 5 D 3/4 for D → ∞. (2.31)
The variations of kc , km , ωm and σm with D are displayed in Fig. 2.6, showing a very good agreement between the numerical solution of the equations and this asymptotic analysis.

Comparisons with experiments 2.4.1 Selection of angular frequency and wavenumber

Considering the experimental parameters (see Tab. 2.1 and typical values in Section 2.2), the dimensionless rigidity lies in the range 10 -3 -10 -2 . For the analysis of the experimental data, we shall then make use of the scaling laws (2.28) and (2.30) obtained in the limit of small D. Introducing back physical dimensions in these expressions, we obtain

ω m ∼ 2ρ m V, (2.32) k m ∼ ρ 2D 1/3 V 2/3 . (2.33)
The selected angular frequency purely results from the balance between dynamic pressure and inertia. The selected wavenumber results from the balance between dynamic pressure and elasticity. It is interesting to compare the phase velocity ω m /k m with that of the elastic waves in the absence of wind flow. In the latter case, Eq. 2.1 tells us that the dispersion relation is simply ω = D/mk 2 , which corresponds to a velocity, evaluated at the most unstable wavenumber, D/mk m . This scales as k m with V 2/3 , whereas ω m /k m is here predicted to be proportional to V 1/3 . Our main goal is the experimental check of these scaling laws of ω m and k m with the wind velocity V . The influence of gravity g is not accounted for in the theory. Computing the dimensionless ratio mgh 3 /D, where h is the sheet thickness, which compares the gravity-induced stress and the stress due to elastic bending, we can see that gravity is clearly sub-dominant: this number is typically in the range 10 -9 -10 -8 . However, gravity does break the up-down symmetry by slightly twisting the sheet, and, importantly, it sets a velocity scale v c that breaks the predicted scale-free power laws. Balancing the dynamic pressure ρv 2 c with the weight of the sheet mg per unit surface, we get the characteristic velocity:

v c = mg/ρ.
(2.34)

Assuming that the observed waves correspond to the most unstable mode, we therefore expect a data collapse when plotting data in the following way:

1 2 m ρ g ω ∼ V v c , (2.35) 3 2D mg k ∼ V v c 2/3 . (2.36)
We emphasize that v c is not an adjustable parameter whose value would depend on the experimental setting, but part of the theoretical analysis.

The experimental measurements of the angular frequency and the wavenumber, rescaled as proposed above, are shown in Fig. 2.7. The data collapse is effectively pretty good, especially if one keeps in mind that the dimensionful wavenumber typically vary by a factor of 6 and ω by a factor of 10 at a given wind velocity from the paper to the thin plastic sheet, for which D changes by two orders of magnitude (Tab. 2.1). The expected 2/3-power of k with the wind velocity is nicely consistent with the data, although the limited accessible range of V gives a low sensitivity on the value of the exponent. The adjusted multiplicative factor in front of (V /v c ) 2/3 is furthermore only 20% below the prediction. The collapse for ω is less impressive and one observes that the expected proportionality relationship (2.35) only holds for large velocities: extrapolating the data, ω would vanish for V ≃ 1.5v c . Such a multiplicative factor of order one shows that the dimensional analysis of perturbative effects is correct. Similarly, the slope of this linear law is around 2, which is the correct order of magnitude, but quantitatively larger than the prediction.

Finite amplitude effects

Although the results displayed in Fig. 2.7 show a good agreement of the selection of angular frequency and wavenumber in the experiments with the prediction of the linear stability analysis of the problem, we have also found some evidence for finite amplitude effects.

Focusing on the paper material, we have systematically varied the sheet length. Data corresponding to different values of ∆L are displayed in Fig. 2.8, showing ω and k in the same rescaled way as in Fig. 2.7. The scaling law obeyed by the angular frequency is found to be independent of the sheet length, whereas that of the wavenumber shows small but systematic variations with ∆L. This shows the presence of non-linearities that are not described here, the linear regime corresponding to the limit of vanishing ∆L. In Fig. 2.9 Schematic for the geometrical relation between the arclength and the length between two lamps. fact, wavenumber and amplitude of the waves can be related to each other as follows. As in Fig. 2.9, we take a sinusoidal shape ζ = Asin(kx) for the sheet over its entire length between the two clamps, and the arc length is then

L 0 0 ( 1 + ζ ′2 )dx.
Therefore, the geometrical constraint that the extra-length ∆L accommodates these undulations without any longitudinal extension can be written as:

∆L = L -L 0 = L 0 0 ( 1 + ζ ′2 -1)dx.
(2.37)

In the regime of small perturbation where kA ≪ 1, and assuming that L 0 is much larger than the wavelength, the integral 1 L 0 sin(kx)dx vanishes and the above relation can be simplified into: ∆L

L 0 = k 2 A 2 4 . (2.38)
This suggests to take data such as those displayed in Fig. 2.2d, and to produce rescaled amplitude profiles of the form

L 0 ∆L kA(x) 2
vs x L 0 . This is done in Fig. 2.10, with data corresponding to different ∆L and different wind velocities. Although scattered, the data collapse and the order of magnitude indicate that this geometrical constraints capture the finite amplitude selection of the waves in these experiments.

Results and Discussions

Combining experiments and a linear analysis in the study of propagative waves on a flexible sheet submitted to a permanent wind, we have shown that the selection of the frequency and of the wavenumber obey the derived scaling laws. The experiments have been performed with different materials, varying the their rigidity by two orders of magnitude. We have shown that these scaling laws result from the balance between dynamic pressure and inertia or elasticity. However, we have here performed the simplest theoretical analysis, based on an unbounded homogeneous sheet. As a consequence, the theory can only work in the limit kL 0 ≫ 1. The boundary conditions actually break the invariance along the x-axis. In principle, one should therefore find the temporal modes whose spatial shape Finally, this work can be of applicative interest for energy harvesting. It is generally based on a fluid flow, or on surface waves inducing a relative motion between articulated parts, which is then converted into an electrical current. However, these moving parts are subject to mechanical wear and may be noisy. It would therefore be interesting to use instead deformable systems without rotating parts, like those investigated here: using the permanent relative flow between the device and the surrounding fluid to produce energy [100-104, 99, 105]. As in Fig. 2.11, it shows an energy harvesting mechanism by using a piezoelectric flag waving in a uniform axial flow. The hydrodynamic flow deforming the surface on which electrical charges are deposited will lead to a motion of charges, and therefore to a current. Such a soft system therefore enables to transform mechanical energy into electrical energy.

Unexpected bedforms on the comet

Rosetta misson has triggered a number of research reporting new discoveries on 67P, including aspects on the nucleus [START_REF] Pätzold | A homogeneous nucleus for comet 67p/churyumov-gerasimenko from its gravity field[END_REF][START_REF] Filacchione | Exposed water ice on the nucleus of comet 67p/churyumovgerasimenko[END_REF][START_REF] Auster | The nonmagnetic nucleus of comet 67p/churyumov-gerasimenko[END_REF][START_REF] Massironi | Two independent and primitive envelopes of the bilobate nucleus of comet 67p[END_REF][START_REF] Sierks | On the nucleus structure and activity of comet 67p/churyumovgerasimenko[END_REF], the coma [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67p/churyumovgerasimenko[END_REF][START_REF] Haessig | Time variability and heterogeneity in the coma of 67p/churyumovgerasimenko[END_REF][START_REF] Rotundi | Dust measurements in the coma of comet 67p/churyumovgerasimenko inbound to the sun[END_REF][START_REF] Nilsson | Birth of a comet magnetosphere: A spring of water ions[END_REF], the surface features [START_REF] Mottola | The structure of the regolith on 67p/churyumov-gerasimenko from rolis descent imaging[END_REF][START_REF] Thomas | The morphological diversity of comet 67p/churyumov-gerasimenko[END_REF] and the interior properties [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF][START_REF] Vincent | Large heterogeneities in comet 67p as revealed by active pits from sinkhole collapse[END_REF][START_REF] Rubin | Dust measurements in the coma of comet 67p/churyumovgerasimenko inbound to the sun[END_REF][START_REF] Altwegg | 67p/churyumov-gerasimenko, a jupiter family comet with a high d/h ratio[END_REF].

What interests us is the bedforms there. The OSIRIS imaging instrument on board the ESA's Rosetta spacecraft has revealed some bedforms on the neck of the comet 67P/Churyumov-Gerasimenko [START_REF] Thomas | The morphological diversity of comet 67p/churyumov-gerasimenko[END_REF][START_REF] El-Maarry | Regional surface morphology of comet 67p/churyumovgerasimenko from rosetta/osiris images[END_REF][START_REF] Thomas | Redistribution of particles across the nucleus of comet 67p/churyumov-gerasimenko[END_REF][START_REF] De Sanctis | The diurnal cycle of water ice on comet 67p/churyumovgerasimenko[END_REF] and on both lobes (Fig. 3.2). Several features suggest that these rhythmic patterns belong to the family of ripples and dunes [START_REF] Charru | Sand ripples and dunes[END_REF]. The pattern presents the characteristic asymmetry of ripples and dunes, with a small steep side and a five times longer gentle slope, which appears darker in Fig. 3.2a. Moreover, two photographs of the same location -one well before perihelion and the other well after it -show the evidence for the bedform propagation (Yellow marks in Fig. 3.2a).

However, the existence of sedimentary bedforms on a comet, comes as a surprise -it requires sediment transport along the surface, i.e. erosion and deposition of particles. When heated by the sun, the ice at the surface of comets sublimates into gas, that is, outgassing.

As gravity is extremely small due to the kilometer scale of the 67P's dimensions [START_REF] Sierks | On the nucleus structure and activity of comet 67p/churyumovgerasimenko[END_REF][START_REF] Pätzold | A homogeneous nucleus for comet 67p/churyumov-gerasimenko from its gravity field[END_REF],

the escape velocity is much smaller than the typical thermal velocity. Outgassing therefore constitutes an extremely rarefied atmosphere, called the coma as in planetary case, around the nucleus (Fig. terrestrial analogue? Making a naive analogy with sediment transport processes on larger planetary bodies, the analog of aeolian dunes [START_REF] Charru | Sand ripples and dunes[END_REF][START_REF] Claudin | A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples[END_REF][START_REF] Elbelrhiti | Field evidence for surface-wave-induced instability of sand dunes[END_REF], whose wavelength scales as λ ∼ 10ρ p d/ρ, would have an emergent wavelength of 10 8 m due to the extremely large density ratio on the comet, i.e. much larger than the comet itself. Similarly, using the comet's values, the analogue for aeolian ripples [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF], whose wavelength scales as λ ∼ ρ/ρ p u 2 t /g, would produce a pattern of wavelength 10 4 m, also larger than the comet size.

Outline of the part

Our goal is to understand the emergence and evolution of the bedforms on 67P, and to constrain the modelling of dynamical processes in the superficial layer of the comet nucleus. To explain the nature of these patterns, a good understanding of the sediment transport on 67P's surface is necessary, which is based on the discussion of the transport threshold and mode. Yet, a description of the atmosphere should come first, which offers the hydrodynamical parameters, like vapor density and viscosity, for the determination of transport threshold and mode. This part is therefore constructed as follows. We first model the dynamics of 67P's atmosphere in Chapter 4. Then, the transport threshold and mode are determined in Chapter 5. Finally, we discuss the nature of the bedforms in Chapter 6.

Chapter 4

Thermo-hydrodynamics of comet 67P/Churyumov-Gerasimenko's atmosphere In this chapter, we propose a model for the dynamics of 67P's atmosphere of our point view, to serve the basics for the subsequent discussions on the sediment transport and on the nature of the bedforms. Despite this two-lobe shape, we work below in spherical coordinates, simplifying the geometry of the comet to a sphere of effective radius R c . We denote by r the radial coordinate that originates at the centre of the nucleus, by θ the ortho-radial (azimuthal) angle, and by ϕ the polar angle. We shall also make use of the distance z to the comet's surface, counted positive downwards. Furthermore, we neglect the effect of the comet's obliquity. More detailed and recent modelling can be found in [START_REF] Bieler | Abundant molecular oxygen in the coma of comet 67p/churyumovgerasimenko[END_REF][START_REF] Haessig | Time variability and heterogeneity in the coma of 67p/churyumovgerasimenko[END_REF][START_REF] Rotundi | Dust measurements in the coma of comet 67p/churyumovgerasimenko inbound to the sun[END_REF][START_REF] Nilsson | Birth of a comet magnetosphere: A spring of water ions[END_REF].

In the following, we first provide derivations to estimate the gravity acceleration, which enters the hydrodynamical equations of the coma, in the region of the neck, where the bedforms that we have primarily studied are located. Then, the thermal process is studied.

Finally, the hydrodynamics of the coma is discussed, and the vapour flow along the surface is computed, which could be responsible for the sediment transport (Chapter 5).

Gravity

The gravity field on 67P has been studied by [START_REF] Sierks | On the nucleus structure and activity of comet 67p/churyumovgerasimenko[END_REF]. The thick neck relating the two lobes is of radius R n ≃ 1 km. The large lobe has dimensions of 4.1 × 3.2 × 1.3 (in km). It can be approximated as a sphere of effective radius R l = (4.1 × 3.3 × 1.8) 1/3 /2 ≃ 1.5 km, leading to a gravity acceleration at the surface g l = 4π

3 ρ c R l ≃ 1.9 10 -4 m/s 2 , where = 6.67 10 -11 m 3 kg -1 s -2 is the gravitational constant and ρ c ≃ 470 kg/m 3 an estimate of the comet's bulk mass density. Similarly, the small lobe is 2.6×2.3×1.8 (in km), which gives an effective radius R s ≃ 1.1 km, and a gravity acceleration at the surface g s ≃ 1.5 10 -4 m/s 2 .

In the region of the neck, the gravity acceleration is given by

g n = (g l sin θ l + g s sin θ s ) 2 + (g l cos θ l -g s cos θ s ) 2 1/2 , ( 4.1)
where we have defined the two angles tan θ l = R n /R l and tan θ s = R n /R s . This expression gives g n ≃ 2.2 10 -4 m/s 2 . This value leads to an escape velocity on the order of g n R n ≃ 0.5 m/s, which is three orders of magnitude smaller than the thermal velocity V th ≃ 500 m/s, as we have declared in Section 3.2.

Thermal process of the comet's nucleus

We describe in this section the thermal and ice sublimation processes, which are responsible for the vapour outgassing near the comet surface and also couple into the hydrodynamics of the coma.

Thermal diffusion

Inside the nucleus, we write the heat conductive flux as ⃗ J = -k c ⃗ ∇T , where T is the temperature field and k c is the thermal conductivity. Denoting by C the bulk heat capacity of the comet and ρ c its bulk mass density, the heat conservation equation reads:

ρ c C∂ t T = k c ∇ 2 T. (4.2)
All three parameters k c , C and ρ c are assumed to be homogeneous. Equivalently, a temperature diffusion equation can be written with a thermal diffusivity κ c = k c /(ρ c C). The material constituting the bulk of the comet is a mixture of dust and ice, with a rather large porosity on the order of 75% [START_REF] Sierks | On the nucleus structure and activity of comet 67p/churyumovgerasimenko[END_REF]. Its effective thermal inertia I = k c ρ c C has been estimated to be in the range 10-50 Jm -2 K -1 s -1/2 [START_REF] Gulkis | Subsurface properties and early activity of comet 67p/churyumovgerasimenko[END_REF]. Taking ρ c ≃ 470 kg/m 3 and

C ≃ (1 -) × 10 3 J/kg/K, we obtain k c ≃ 10 -2 W/m/K and κ c ≃ 10 -7 m 2 /s.
The time-evolution of the temperature of the comet's surface T s can be decomposed in Fourier modes. Diffusion being linear, we can do the reasoning one particular mode of angular frequency ω, written in complex notations as Ts (ω). Assuming that the flux vanishes at infinity (deep inside the bulk of the comet), the temperature field takes the form:

T = T s exp (iωt + qz) (4.3)
Introducing Eq. 4.3 into Eq. 9.1, one obtains:

q 2 = iω κ , (4.4)
and the negative root is chosen since the temperature decays to zero when z → ∞.

Therefore, the solution of the diffusion equation in the Fourier space for the temperature field takes the form:

T (z, ω) = Ts (ω) exp -(1 -i)z |ω| 2κ c for ω ≤ 0, (4.5) T (z, ω) = Ts (ω) exp -(1 + i)z |ω| 2κ c for ω > 0. (4.6)
The penetration length δ is defined:

δ = 2κ c |ω| . ( 4.7)
The rotation period of the comet is Γ d = 12.4 hours, or, equivalently, ω d = 2π/Γ d = 1.4 10 -4 s -1 . This gives a diurnal penetrating length δ d ≃ 4 cm, which means that a few tens of cm below the surface, the day-night alternation has no influence on the temperature field. Regarding the seasonal variations, the orbital period is Γ y = 6.44 years, corresponding to a penetrating length δ y ≃ 3 m. Conversely, one can compute the time scale corresponding to the size of the comet δ h = R c , which gives Γ h ≃ 10 6 years. This is the time scale required to get a homogeneous temperature T a across the whole body. It is much smaller than the age of the comet, which is that of the solar system, i.e. about 4.5 10 9 years.

Ice sublimation

We hypothesis that the vapour outgassing comes from the sublimation of ice just below the surface of the comet. To sublimate ice at a rate corresponding to a vapour mass flux q m (in kg per second and per unit surface), a power per unit surface q m is absorbed.

≃ 3 10 6 J/kg is the latent heat of water ice sublimation. The corresponding power balances write:

(1 -Ω)ψ = σεT 4 s + J s + q m , (4.8)
where Ω = 0.05 is the estimated albedo, σ = 5.67 10 -8 W/m 2 /K 4 is the Stefan constant and ϵ ≃ 0.9 is the estimated emissivity [START_REF] Gulkis | Subsurface properties and early activity of comet 67p/churyumovgerasimenko[END_REF].

ψ = sin ϕ ψ E (η E /η) 2 φ
J s = -k c ∂ z T | z=0 . (4.9)
Pluging Eq. 4.3 into Eq. 4.9, one easily gets J s from its Fourier transform:

Ĵs (ω) = (1 -i)k c |ω| 2κ c Ts (ω)
for ω ≤ 0, (4.10)

Ĵs (ω) = (1 + i)k c |ω| 2κ c Ts (ω) for ω > 0. (4.11)
The integration of Eq. 4.8, coupled to those describing the vapor flow in the atmosphere as well as in the porous surface layer, is used to predict the time variations of the vapor flux q m at both daily and yearly scales.

Hydrodynamics of the comet's atmosphere

The vapour flow in the comet's atmosphere is described by the conservation of mass, momentum and energy [START_REF] Kundu | Fluid Mechanics[END_REF]:

∂ ρ ∂ t + ⃗ ∇ • (ρ⃗ u) = 0, (4.12
)

∂ ρ⃗ u ∂ t + ⃗ ∇ • (ρ⃗ u⃗ u) = ρ ⃗ g -⃗ ∇p + ⃗ ∇ • ⃗ ⃗ τ, (4.13) ∂ ∂ t ρ ε + 1 2 u 2 + ⃗ ∇ • ρ w + 1 2 u 2 ⃗ u = ρ ⃗ g • ⃗ u + ⃗ ∇ • ( ⃗ ⃗ τ • ⃗ u) -⃗ ∇ • ⃗ J, ( 4.14)
with the mass density ρ, the velocity ⃗ u, the pressure p, the stress tensor ⃗ ⃗ τ, the specific energy ε, the specific enthalpy w = ε + p/ρ, the heat flux ⃗ J and the gravity acceleration ⃗ g.

Taking the density weighted time averaging to get so-called Favre averaged Navier Stokes (FANS) equations, the averaged stress tensor can be expressed as the sum of viscous and turbulent contributions:

τ i j = ρνγ i j + ρν t γi j - 1 3 Kδ i j , (4.15)
where we have introduced the shear rate γi j = ∂ j u i + ∂ i u j -2 3 ∂ k u k δ i j , and where K = χ 2 |γ|, with χ ≃ 2.5 a phenomenological constant and |γ| the modulus of the shear rate tensor. In the ideal gas approximation, the molecular viscosity ν can be related to the mean free path

ℓ = m 2πd 2 w ρ , ( 4.16)
and to the thermal velocity

V th = 8k B T πm , (4.17)
defined as the mean magnitude of the velocity of the molecules, by

ν = 1 3 V th ℓ. ( 4
J i = -ρ γ γ -1 ν Pr + ν t Pr t ∂ i p ρ , (4.19)
with γ = 4/3 the adiabatic expansion coefficient and ν t the turbulent viscosity of water vapour. Pr and Pr t are the Prandtl and turbulent Prandtl numbers, both typically on the order of unity for gases. The averaged energy density has also an internal and a turbulent contribution:

e = ρε = 1 γ -1 p + 1 2 ν t ρK. (4.20)
Finally, the additional term u j τ i j complements the enthalpy contribution ρwu i .

Outer layer flow

Eqs. (4.12-4.14) can be solved averaging over the polar angle, and assuming steady state. We describe the atmosphere as a two-layer flow: an outer layer where viscosity and turbulent fluctuations can be neglected (perfect flow) and an inner turbulent layer of thickness δ i ≪ R c matching with the surface conditions. Here we take R c ≃ 1.95 km, corresponding to an equivalent surface S c ≃ 47.7 km 2 [START_REF] Gulkis | Subsurface properties and early activity of comet 67p/churyumovgerasimenko[END_REF]. We separately note U r and U θ the velocity components in outer layer, and u r and u θ those in the inner layer. The pressure and density in the inner layer are inherited from the outer layer. This hydrodynamical description of the comet's atmosphere loses it validity when the mean free path of the vapour becomes on the order of the comet size itself. Neglecting all dissipative terms in (4.12-4.14), the steady equations for the outer layer are, for mass conservation:

1 r 2 ∂ ∂ r r 2 ρU r + 1 r ∂ ∂ θ (ρU θ ) = 0; (4.21)
for momentum conservation in the radial direction:

1 r 2 ∂ ∂ r r 2 ρU 2 r + 1 r ∂ ∂ θ (ρU r U θ ) - 1 r ρU 2 θ + ∂ p ∂ r = 0; (4.22)
for momentum conservation in the ortho-radial direction:

1 r 3 ∂ ∂ r r 3 ρU r U θ + 1 r ∂ ∂ θ ρU 2 θ + 1 r ∂ p ∂ θ = 0; (4.23)
and for the energy conservation:

1 r 2 ∂ ∂ r r 2 1 2 ρ U 2 r + U 2 θ + γ γ -1 p U r + 1 r ∂ ∂ θ 1 2 ρ U 2 r + U 2 θ + γ γ -1 p U θ = 0. (4.24)
The asymptotic analysis of these equations gives:

U r ∝ r 0 , U θ ∝ r 2(1-γ) , ρ ∝ r -2 and p ∝ r -2γ . ( 4.25)
One concludes that orthoradial terms are subdominant in the outer layer, and the equations at the leading order reduce to

1 r 2 ∂ ∂ r r 2 ρU r = 0, (4.26) 1 r 2 ∂ ∂ r r 2 ρU 2 r + ∂ p ∂ r = 0, (4.27) 1 r 2 ∂ ∂ r r 2 1 2 ρU 2 r + γ γ -1 p U r = 0. (4.28)
These equations can be integrated as:

U r = U 0 G(r), (4.29) ρ = ρ 0 R c r 2 1 G(r)
, (4.30)

p = p 0 + γ -1 2γ ρ 0 U 2 0 [1 -G(r)] R c r 2 1 G(r) , (4.31) = p 0 G ∞ -G(r) G ∞ -1 R c r 2 1 G(r) (4.32)
where the function G satisfies G(R c ) = 1, so that ρ 0 and p 0 are the vapour density and the pressure at the surface of the comet r = R c and U 0 the velocity at top of the surface layer.

We have introduced

G ∞ = 1 + 2γ γ -1 p 0 ρ 0 U 2 0 . ( 4
G ′ - γ -1 γ + 1 G ∞ G ′ G + (G ∞ -G) 4 r = 0, (4.34)
This first order differential equation solves into:

G 1 2 (γ-1) G ∞ -G G ∞ -1 = R c r 2(γ-1) , (4.35)
expressed inexplicitly, which could be easily handled in a numerical way. The outer flow is then entirely determined by the three parameters ρ 0 , U 0 and p 0 .

Turbulent boundary layer

We need to compute the vapor wind flow close to the surface, which may entrain the surface grains into motion. This flow is controlled by the momentum balance in the boundary layer approximation, in which the horizontal diffusion of momentum is negligible:

1 r 3 ∂ ∂ r r 3 (ρu r u θ -τ rθ ) + 1 r ∂ ∂ θ ρu 2 θ + 1 r ∂ p ∂ θ = 0. (4.36)
To compute an approximate solution, we write the velocity profile in the inner layer under the form:

u θ (r) = u * κ ln 1 + r -R c z 0 , (4.37)
parametrised by the shear velocity u * defined from the basal shear stress τ 0 rθ ≡ ρ 0 |u * |u * , and the phenomenological von Kármán constant κ ≃ 0.4 has been calibrated on turbulent flows along a wall. For the sake of simplicity, we will use here the logarithmic law of the wall, but more complicated profiles could be easily accommodated. z 0 is the aerodynamic roughness and here we take z 0 = 0.11ν/u * corresponding to the smooth aerodynamic regime. We introduce the notation

Λ ≡ ln 1 + δ i z 0 , (4.38)
where δ i is the thickness of the boundary layer. The density and the pressure in the inner layer are inherited from the outer layer, so that p ≃ p 0 and ρ ≃ ρ 0 . The radial velocity at the top of the of the boundary layer is U 0 . Integrating (4.36) between R c and R c + δ i , for δ i ≪ R c , we obtain:

ρ 0 |u * |u * + ρ 0 Λ κ U 0 u * + d dθ 2 -2Λ + Λ 2 δ i κ 2 R c ρ 0 u 2 * + δ i R c dp 0 dθ = 0, (4.39)
where we have used the fact that the velocity u θ vanishes at the comet's surface, and that the shear stress vanishes at the top of the inner turbulent boundary layer, when one reaches the outer perfect flow.

The radial component of the velocity in the inner layer u r is deduced from u θ by the mass conservation equation:

1 r 2 ∂ ∂ r r 2 ρu r + 1 r ∂ ∂ θ (ρu θ ) = 0. (4.40)
By integration across the boundary layer, we similarly obtain:

U 0 = u 0 - 1 ρ 0 d dθ (Λ -1)δ i κR c ρ 0 u * . (4.41)
Using this expression for U 0 in (4.36), we deduce:

ρ 0 |u * |u * + ρ 0 Λ κ u 0 u * + d dθ 2 -2Λ + Λ 2 δ i κ 2 R c ρ 0 u 2 * - Λ κ u * d dθ (Λ -1)δ i κR ρ 0 u * = - δ i R c dp 0 dθ , (4.42)
The boundary layer thickness corresponds to the crossover altitude at which one makes the transition from the inner to the outer layer, i.e. where the inertial terms are comparable to the pressure gradient:

|u * | + Λ κ u 0 ρ 0 u * ≈ d dθ 2 -2Λ + Λ 2 δ i κ 2 R c ρ 0 u 2 * - Λ κ u * d dθ (Λ -1)δ i κR c ρ 0 u * , (4.43)
so that (4.42) simplifies into:

- δ i 2R c dp 0 dθ = |u * | + Λ κ u 0 ρ 0 u * . (4.44)
Eq. 4.43 is further simplified under the assumption that variations of all quantities along θ are slow, essentially equivalent to sinusoidal variations, i.e. with d dθ ≈ 1 2π . We then obtain:

(Λ -2)δ i 2πκ 2 R c ≈ |u * | u * + Λ κ u 0 u * . ( 4

Porous sub-surface layer

The vapor production rate from outgassing, defined as the product of the vapor density ρ 0 by the outward vapor velocity u 0 , has been measured at different heliocentric distances [START_REF] Gulkis | Subsurface properties and early activity of comet 67p/churyumovgerasimenko[END_REF][START_REF] Hanner | The dust coma of periodic comet churyumov-gerasimenko[END_REF][START_REF] Osip | Comets: groundbased observations of spacecraft mission candidates[END_REF][START_REF] Schleicher | Compositional and physical results for rosetta's new target comet 67p/churyumov-gerasimenko from narrowband photometry and imaging[END_REF][START_REF] Ootsubo | Akari near-infrared spectroscopic survey for cometary volatiles[END_REF][START_REF] Bertaux | The water production rate of rosetta target comet 67p/churyumov-gerasimenko near perihelion in 1996,?2002 and 2009 from lyman α observations with swan/soho[END_REF] (Fig. 4.5). Common models assume that ice sublimation takes place at the surface and produces a radial flow at the thermal velocity [START_REF] Skorov | Simulation of gas flow in a cometary knudsen layer[END_REF]. This would result in a density ρ 0 an order of magnitude smaller than that necessary to induce a fluid drag force large enough to overcome the threshold for grain motion (as discussed in Chapter 5). We suggest that most of the vapor is emitted from sub-surface ice and must travel through the porous surface granular layer (Fig. 4.1). Sublimation makes the ice trapped in the pores recede, releasing unglued grains in surface that can be eroded. We assume that the ice level remains at a constant distance from the surface, comparable to the grain size d.

We describe in the following the close sub-surface as a thin porous granular layer of thickness h. The picture is that of a chaotic billiard, where a water molecule, emitted at depth z = h where the ice is, experiences collisions with the grains of the packing but not with the other molecules. The mean free path of the molecules is then a fraction of grain size d. The probability for a molecule to cross this layer rather than going back to z = h and being adsorbed by the ice again is p c ∝ d/h, depending on porosity and grain shape.

We assume that the water molecules emitted from ice have a half Maxwell-Boltzmann velocity distribution:

P i (⃗ v) = m 2πk B T i 3/2 exp - m|⃗ v| 2 2k B T i Θ(⃗ v • ⃗ e r ). (4

.46)

T i is the temperature of the ice at z = h. Θ is the Heaviside function and ⃗ e r is the unit vector pointing upwards. The vapor mass flux of molecules emitted by the ice surface is then

F ρ sat +∞ -∞ d v x +∞ -∞ d v y +∞ 0 d v r v r P i (⃗ v) = 1 4 F ρ sat V i th , (4.47)
Fig. 4.1 Schematics of the porous granular layer at the comet's surface. The water molecules are emitted by the ice (dark blue) at the thermal velocity corresponding to the ice temperature. Experiencing collisions with the grains of the packing (blue arrows), the molecules have a probability to cross the layer decreasing as the inverse of its thickness. The mean free path of the molecules in the layer is comparable to the pores between the grains, i.e. a fraction of the grain size. Molecules just above the surface may also enter the porous layer and be absorbed if they reach the ice. This layer is typically 1.5d thick, so that the surface grains, not glued to ice, are potentially free to move if the wind is above the transport threshold.

where F is the ice surface fraction, and where we have introduced the thermal velocity Eq. 4.17). ρ sat is the saturated vapor density, here also evaluated at the temperature of the ice T i .

V i th = V th (T i ) = 8k B T i /(πm) (see
At the comet's surface (z = 0), where the temperature of the vapor is T 0 , we assume furthermore that the vapor flow has an average velocity u 0 ⃗ e r , so that the water molecules have a velocity distribution given by:

P 0 (⃗ v) = m 2πk B T 0 3/2 exp - m|⃗ v -u 0 ⃗ e r | 2 2k B T 0 . (4.48)
The vapor mass flux of molecules entering in the porous layer from the atmosphere, whose density is ρ 0 , is then

q -= ρ 0 +∞ -∞ d v x +∞ -∞ d v y 0 -∞ d v r (-v r )P 0 (⃗ v) = 1 4 f (Υ 0 )ρ 0 V 0 th , (4.49)
where we have introduced the thermal velocity V 0 th = V th (T 0 ), the velocity ratio Υ 0 ≡ u 0 /V 0 th and defined the function:

f (Υ ) = e -4Υ 2 π -2Υ 1 -erf 2Υ π . (4.50)
Υ is similar to a Mach number, as the speed of sound in an ideal gas is γk B T /m = π/6 V th for an adiabatic index γ = 4/3 used here.

Assuming perfect absorption of the water molecules when they come back to ice (a vanishing probability of rebound), the vapor mass flux coming out at the surface q m = ρ 0 u 0 is then the result of the following balance:

q m = Υ 0 ρ 0 V 0 th = p c F 4 ρ sat V i th -q -. (4.51)
In the limit of an unlimited (F = 1) and vanishingly thin (T i = T 0 ) layer, the Hertz-Knudsen sublimation law, with a vapur flux proportional to (ρ sat -ρ 0 )V th is recovered. Similarly, the momentum flux ρ 0 u 2 0 + p 0 reads:

Υ 2 0 + π 8 ρ 0 V 0 th 2 = π 4 1 4 F p c ρ sat V i th 2 + (2 -p c )q -V 0 th . (4.52)
Finally, the energy flux 1 2 ρ 0 u 2 0 + γ γ-1 p 0 u 0 reads:

1 2 Υ 0 Υ 2 0 + π ρ 0 V 0 th 3 = 7π 16 p c 1 4 F ρ sat V i th 3 -q -V 0 th 2 . (4.53)
Introducing the expression for q -(4.49) into Eqs. 4.51 and 4.52, we solve for ρ 0 and V 0 th :

ρ 0 = F p c π[ f (Υ 0 )(p c -2) + 2] + 16Υ 2 0 π[ f (Υ 0 )p c + 4Υ 0 ] 2 ρ sat , (4.54)
V 0 th = π[ f (Υ 0 )p c + 4Υ 0 ] π[ f (Υ 0 )(p c -2) + 2] + 16Υ 2 0 V i th . (4.55)
The final equation for Υ 0 is obtained introducing these expressions into (4.53):

-7π 2 + 32π 2 -112π Υ 2 0 + (32π -448) Υ 4 0 + 7(p c -1) f 2 (Υ 0 ) + (14 -7p c )π 2 + 15p c π 2 Υ 0 + (112π -56p c π)Υ 2 0 + 8p c πΥ 3 0 f (Υ 0 ) = 0. (4.56)
To solve numerically this equation, values must be chosen for the different parameters.

Consistently with the value of the porosity of the comet's ground, we take F = 0.2 for the ice surface fraction. The porous layer thickness is set to h = 1.5 d, which corresponds to a monolayer of grains not attached to the icy bed, and free to move by the wind. The probability for a water molecule to cross the porous layer is set to p c = 0.1d/h ≃ 0.07, in order to adjust the vapor density at the comet's surface. With these numbers, the velocity ratio Υ 0 = u 0 / 8k B T 0 /(πm), which compares the outgassing velocity to the thermal velocity of the vapour at the comet's surface can be computed as the solution of Eq. 4. it is evaluted by the temperature at the ice surface [START_REF] Murphy | Review of the vapour pressures of ice and supercooled water for atmospheric applications[END_REF]. The results will be discussed in the following.

Results and discussions

Temperatures, density and pressure

The calculated results show both seasonal and diurnal time variations of the atmosphere characteristics. At perihelion, the vapor temperature peaks around 200K at the surface of the comet (Fig. 4.2), and the corresponding thermal velocity around 500 m/s, which is much larger than the escape velocity. The vapour mean free path ℓ is about 3 cm at the surface of the comet (Fig. 4.3). As ℓ is significantly smaller than the bedform emergent wavelength λ ≃ 5 m (Fig. 3.2a), hydrodynamics accurately describes the flow above relief.

The vapor density is around 10 times larger than the previous estimations due to the presence of the granular porous layer. As we can see, there exists an asymmetry between sunrise and sunset for temperature, density as well as the pressure (Fig. 4.4), which is simply results from thermal inertia, as some heat is accumulated in the superficial layer during in the morning and released in the afternoon.

Vapour flux

Regarding the vapor flux, observations [START_REF] Gulkis | Subsurface properties and early activity of comet 67p/churyumovgerasimenko[END_REF][START_REF] Hanner | The dust coma of periodic comet churyumov-gerasimenko[END_REF][START_REF] Hanner | Compositional and physical results for rosetta's new target comet 67p/churyumov-gerasimenko from narrowband photometry and imaging[END_REF][START_REF] Bertaux | The water production rate of rosetta target comet 67p/churyumov-gerasimenko near perihelion in 1996, 2002 and 2009 from lyman α observations with swan/soho[END_REF][START_REF] Ootsubo | Akari near-infrared spectroscopic survey for cometary volatiles[END_REF] provide data at a different heliocentric distances η (Fig. 4.5a), which we use to calibrate some parameters of the model. From the vapor mass flux q m coming out at the surface, integrated over the whole comet, the global vapor flux reads:

qm (η) = α 4π π -π dθ π 0 sin ϕ q m (θ , ϕ) dϕ, ( 4.57)
where the factor α accounts for the fraction of the surface where sublimation is effective.

Assuming that all points of the surface receiving the same insolation would produce the same vapor rate, one can solve Eq. 4.8 at the equator only (ϕ = π/2) and compute the vapor rate as

qm (η) = α 4 π -π | sin θ | q m (θ ) dθ , (4.58)
where the angle θ = 0 points in the direction of the sun. This assumption is valid as long as the heat flux term J s in (4.8) is negligible, so that the surface points can be considered as thermally decoupled. This is the case in the illuminated side of the comet

(-π/2 ≤ θ ≤ π/2)
, where most of the vapor flux comes from. This approximation is uncontrolled on the night-side, where J s , due to the thermal inertia of the comet's body, is the source of heat for sublimation, but corresponding to a negligible part of q m (Fig. 4.5b).

The fit of the observational data allows us to set the porous layer thickness to h = 1.5 d.

Larger values lead to a dependence of the vapor flux that decreases to fast with the heliocentric distance η. Also, the fraction of active (sublimating) surface is adjusted to α = 0.1 in order to reproduce the value of the flux at perihelion. , where θ = 0 corresponds to the noon. Therefore it is night for 67P for θ = π/2 ∼ 3π/2, when q m is negligible.

Wind

From the thermal calculation, we see that the pressure drops by ten orders of magnitude from day to night (Fig. 4.4). The comet's atmosphere therefore presents a strong pressure gradient that drives a tangential flow from the warm, high pressure towards the cold, low pressure regions, in a surface boundary layer (Fig. 4.6). The extension of the halo of vapor on the dark side of the comet is a signature of this surface wind (Fig. 3.3). To get u * , we iteratively solve Eqs. 4.38, 4.44 and 4.45 by pluging in the density ρ 0 (θ ) and the pressure p 0 (θ ) profiles from the thermal calculations. The calculated results show that u * reverses direction during the day and is maximal at sunrise and sunset due to the direction change Fig. 4.6 Schematic of the vapour flow at the comet surface driven by strong pressure gradients from illuminated to shadow areas, with the solar radiation labelled in yellow, the outgassing in blue and the resulting winds in red.

of the pressure gradient, with a shear velocity u * on the order of a fraction of the thermal velocity (Fig. 4.2).

A brief summary

We have proposed here a model for the thermo-hydrodynamics of 67P's coma. It couples the kinetics of ice sublimation, driven by the power per unit surface received from the sun, to the hydrodynamical description of the vapour flow in the coma. The calculation is performed in an ideal spherical geometry using the active fraction 10% of the surface as an adjustable parameter to match the measurements of the vapour rate at different heliocentric distances. Using kinetic theory of gas, we predict that for such vapor flow the outgassing velocity is around ten times smaller than that of the spectacular vapor jets streaming from active pits [START_REF] Vincent | Large heterogeneities in comet 67p as revealed by active pits from sinkhole collapse[END_REF][START_REF] Sierks | On the nucleus structure and activity of comet 67p/churyumovgerasimenko[END_REF]. Accordingly, the vapor atmosphere is ten times denser than previous estimates. A vapour flow along the comet surface is resulted from the strong pressure gradient at the day-night alternace. Still, the atmosphere density is 7 orders of magnitude lower than on Earth, can such a flow entrain grains into motion? To answer this question, the shear velocity u * is compared with the threshold velocity in the next chapter.

Chapter 5

Sediment transport

Using the thermo-hydrodynamical parameters of the coma obtained in Chapter 4, sediment transport on 67P is discussed in this chapter, to deduce a few key quantities, namely the sediment transport threshold, the saturated flux and the saturation length for the cometary case. We first estimate the grain size on 67P, which is the key parameter determining the sediment transport threshold. Then, we determine the sediment transport mode. Finally, the transport law is derived, accounting for the peculiar conditions of the comet's atmosphere.

Grain size

There are high resolution pictures taken by Philae, Rosetta's lander module, above the Agilkia landing site (Fig. 5.1b). To estimate the grain size, a series of calibrated photographs of a sand-bed is used to relate the image auto-correlation to the mean grain diameter of the bed. The pictures are taken at resolutions going from 1 to 10 pixels per grain diameter d. The rescaled correlation functions C(δ) collapse on a master curve, when δ is properly divided by d -both expressed in the same units. To determine the mean grain size using a picture whose resolution is known, one computes C(δ), with δ expressed in meters or in pixels. One then determines by a least square method the value of d that should be used as rescaling factor of δ, to collapse the new curve on the calibration curve [START_REF] Claudin | A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples[END_REF]. The best collapse of the correlation functions is obtained for a mean grain diameter d ≃ 9.7 mm on the comet (Fig. 5.1A). Another image available for the estimation was taken by Rosetta just before its touch down in the Ma'at region (Fig. 5.1d), which allowed us to determine the surface granulometry (Fig. 5.1C). The results indicate that the surface is composed of centimeter scale grains, which is consistent with the prediction in the publication [START_REF] Mottola | The structure of the regolith on 67p/churyumov-gerasimenko from rolis descent imaging[END_REF].

Transport threshold

In this section, we first calculate the threshold velocity u t with the effects of the mean free path and the intergrain cohesion taken into account. Then u t is compared to u * to evaluate whether the wind is able to set the surface grains into motion. The dependence of u t on the grain size is discussed at the end.

Threshold velocity u t

We consider a grain of size d at the surface of the comet, on the verge to be entrained into motion. The threshold shear velocity u t is quantitatively determined by the balance between gravity, hydrodynamic drag, a cohesive force at the grain contacts and a resistive force associated with the geometrical effect of the surrounding grains. The later can be modeled by a Coulomb friction of coefficient µ relating the tangential and normal forces.

The grain weight can be expressed as π 6 ρ p gd 3 , where g is the gravity acceleration and ρ p is the mass density of the grains. In most practical cases, the threshold velocity falls in the cross-over between the viscous and turbulent asymptotic regimes. It is thus important to have a model of it valid in both regimes [START_REF] Claudin | A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples[END_REF]. The drag force exerted on a grain reads

F drag = π 8 C d d 2 ρu 2 , ( 5.1)
where u is the velocity of the fluid around the grain and C d is a drag coefficient. In order to account for viscous as well as turbulent regimes, C d can conveniently be written as:

C d = C 1/2 ∞ + s ν ud 1/2 2 , ( 5.2)
where ν is the fluid viscosity. C ∞ and s are phenomenological calibrated constants. For example, we have C ∞ ≃ 1 and s ≃ 5 for natural grains. For the cometary case on 67P, the atmosphere is dilute, i.e. the mean free path ℓ becomes comparable to the grain size, the no-slip condition is no longer correct when calculating the drag force on small particles, and one needs to consider the noncontinuum effects. An empirical correction due to Cunningham [START_REF] Cunningham | On the velocity of steady fall of spherical particles through fluid medium[END_REF] is applied here,

s 2 = 25 1 + 2ℓ d (1.257 + 0.4 exp(-0.55d/ℓ)) . (5.3)
When the grain is at rest at the surface of the bed, we consider that the hydrodynamical stress is exerted on its upper half so that the effective drag force becomes

F drag = β π 8 C d d 2 u 2 , x z Fig. 5.2
The origin of the sediment transport threshold at the scale of a grain.

with β = 1/2. Just at the threshold and neglecting cohesion force A for the moment, this force is balanced by the horizontal bed friction felt by the grain: F t = π 6 µ(ρ p -ρ)gd 3 (Fig. 5.2). Here we take µ = tan(29 • ) ≃ 0.55. We introduce the viscous size

d ν = (ρ p /ρ -1) -1/3 ν 2/3 g -1/3 , (5.4)
and further make the fluid velocity dimensionless as 1/2 ≡ u/ (ρ p /ρ -1)gd. With these notations, the threshold value of the flow velocity at the scale of the grain, denoted as

1/2 t , is solution of (C ∞ t ) 1/2 + s d ν d 3/4 1/4 t - 4µ 3β 1/2 = 0, (5.5)
which resolves immediately into:

t = 1 16C 2 ∞ s 2 d ν d 3/2 + 8 µC ∞ 3β 1/2 1/2 -s d ν d 3/4 4
.

(5.6)

Following [START_REF] Claudin | A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples[END_REF], the corresponding threshold Shields number is the sum of a viscous and a turbulent contribution:

Θ t = 2 d ν d 3/2 1/2 t + κ 2 ln 2 (1 + 1/2ξ) t , (5.7)
where ξ is the hydrodynamic roughness rescaled by the grain diameter. Here we take the experimental value ξ = 1/30.

A similar approach can be used to compute the settling velocity V fall , which also gives the vertical threshold velocity, balancing the drag force and the particle weight. We can proceed as in Eq. 5.5, but with µ/β = 1 and get:

fall = 1 16C 2 ∞ s 2 d ν d 3/2 + 8 C ∞ 3 1/2 1/2 -s d ν d 3/4 4 .
(5.8)

The settling velocity is found always smaller than u * during the fraction of time when sediment transport occurs.

Cohesion

The effect of the cohesion is absent from the discussion above, which is important when the grain size is below a certain value. In the following, we consider the cohesion force A between the grains and include it in the calculation of u t . Investigating this balance highlights the need to apply findings from contact mechanics of rough interfaces [START_REF] Greenwood | The elastic contact of rough spheres[END_REF] to the study of small solar system bodies. The adhesive free energy, resulting from van der Waals interactions, is proportional to the real area of contact between the grains, which is much smaller than the apparent one because of surface roughness. A realistic computation of this cohesion can be achieved under the assumption that contacts between grains are made of elastically deformed nano-scale asperities and that the apparent area of contact follows Hertz law for two spheres in contact. The cohesive force is then found to scale as the maximal load experienced by the grains to the power 1/3 [START_REF] Restagno | Adhesion between weakly rough beads[END_REF]. Considering that this load is typically the weight of a surface grain, this force scales as

A ∝ (ρ p gd/E) 1/3 γd, ( 5.9)
where ρ p is the grain bulk density, E is the grain Young modulus and γ is the surface tension of the grain material. It is therefore much lower than the force γd obtained for ideally smooth grains. Importantly, the gravity force increases as d 3 , while the cohesive force increases as d 4/3 only. This allows to define a cross-over diameter d m at which these two forces are comparable:

d m = γ 3 Eρ 2 p g 2 1/5
.

(5.10)

It gives the typical grain diameter below which cohesive effects become important and are responsible for the increase of the threshold at small d. For silica (quartz) grains on Earth, the cohesive size d m is around 10 µm (Fig. 5.3), and this is why sand grains, with typical diameters on the order of a few hundreds microns, are not affected by cohesion. On the comet, the composition of the regolith dust is not precisely know, but the particle bulk density ρ p has been estimated in the range 1000-3000 kg/m 3 [START_REF] Rotundi | Dust measurements in the coma of comet 67p/churyumovgerasimenko inbound to the sun[END_REF], i.e. close to that of sand on Earth. We make the assumption that the values of E and γ are also similar for the particles on both bodies. According to (5.10), the ratio of the values of d m on Earth and on 67P is then essentially given by the corresponding ratio of the gravities, to the power 2/5. Using the gravity field derived above, we can estimate d m ≃ (9.8/0.00022) 2/5 × 10 µm ≃ 720 µm on the comet. Such a millimeter scale is three orders of magnitude smaller than the capillary length γ/ρ p g ≃ 1 m suggested by the traditional approach, which ignores contact roughness [START_REF] Kok | The physics of wind-blown sand and dust[END_REF]. Accounting for these cohesion effects, the threshold Shields number finally reads:

Θ t = Θ 0 t 1 + 3 2 d m d 5/3 , (5.11)
where Θ 0 t is the expression given by Eq. 5.7 [START_REF] Claudin | A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples[END_REF]. And one thus obtains the threshold velocity with cohesive effect taken into account:

u t = Θ t (ρ p /ρ -1)gd.
(5.12)

Apparently, u t is a function of both the wind properties and the grain size, and the two dependences are respectively discussed in the following.

A comparison of u t and u *

Using the wind parameters (ρ 0 , ℓ and ν), we compute u t at the year scale for a given grain diameter d = 4 mm (discussions Chapter 6) with Eq. 5.12. Fig. 5.4b shows |u * | and u t at perihelion, which indicates that where the wind is strong enough to entrain the grains into motion in the 'afternoon' of one comet day. The ratio of u * and u t shows that there is only a small fraction of the time -typically ≃ 6 10 3 s at perihelion, i.e. ≃ 14% of the comet's day of 12.4 h. The asymmetry between sunrise and sunset winds has an important consequence: the morning thermal wind is not strong enough to entrain grains (Fig. 5.4a). With the wind parameters (ρ 0 = 2.6 10 -6 kg/m 3 , ℓ = 3 cm and ν = 4.6 m 2 /s) at the 'afternoon' peak wind at perihelion (Inset of Fig. 5.4a), u t is computed and plotted as a function of the grain size d in Fig. 5.5, and it and shows a minimum value on the order of 50 m/s for the whole range 10 3 -10 5 µm. The curve shows generally three regimes. In the large d regime, the turbulent drag essentially balances the friction force:

Dependence of u t on d

u t ∼ (ρ p /ρ 0 )gd ∝ d 1/2 .
(5.13)

In the small d regime, the viscous drag balances cohesion:

u t ∼ ρ p gd E 1/6 γℓ ρ p gd 3 1/2 (ρ p /ρ 0 )gd ∝ d -5/6 .
(5.14)

In the intermediate regime for which d m < d < ℓ, the mean free path ℓ of the vapor molecules is larger comparing with the case on Earth, which leads to a reduced drag force for grains smaller than ℓ. In this regime, the viscous drag balances the friction force:

u t ∼ (ρ p /ρ 0 )gℓ ∝ d 0 . (5.15)
This explains that u t presents a plateau extending from the millimeter scale to the meter scale (Fig. 5.5).

In conclusion, we find that, sufficiently close to perihelion, all these grains, and in particular those at the centimeter scale observed by Rosetta near bedforms, can be transported by the afternoon thermal wind (Fig. 5.5).

Transport mode and saturated transport

Transport mode

Given the very large density ratio between grains and vapor, the grains can rebound when colliding with the surface of the comet. At first this would suggest the possibility of a cometary saltation [START_REF] Thomas | Redistribution of particles across the nucleus of comet 67p/churyumov-gerasimenko[END_REF] analogous to aeolian transport, in which the grains move by bouncing or hopping [START_REF] Durán | On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws[END_REF][START_REF] Kok | The physics of wind-blown sand and dust[END_REF] (Fig. 1.4a,b). However, compared to the Earth, there are two important differences that prevent saltation (Fig. 5.6). Rebounding grains would acquire a vertical velocity larger than the escape velocity, on the order of a meter per second. Moreover, the flow is turbulent above a viscous sub-layer, typically 0.7 m thick at perihelion, where turbulent fluctuations are dampened by viscosity. After a rebound, grains with enough energy to reach the turbulent zone would be entrained into suspension, as the settling velocity is much smaller than turbulent velocity fluctuations. These grains are eventually ejected into the coma. Thus, the only mode of sediment transport along the surface is traction, where grains remain in contact with the substratum on which they roll or slide.

Traction Ejection

Viscous boundary layer Fig. 5.6 Schematics featuring the modes of sediment transport on the surface of 67P. In the cometary case, grains rebounding on the bed are eventually ejected in the coma, preventing the existence of saltation. This schematics holds for monolithic (crystalline) grains as well as for agglomerates of smaller particles. Violet background: viscous sub-layer close to the bed, which is typically 10ν/u * ≃ 0.7 m thick in the cometary case at perihelion.

Saturated transport flux q sat

The grains on the comet's bed move in the traction mode, a slow transport, where the energy brought by the flow is dissipated during the collision of moving grains with the static grains of the bed. Sediment transport on the comet is therefore analogous to subaqueous bedload. The behavior of q sat for bedload has been measured and fitted by an empirical formula (Fig. 1.5a). Recalling the idea of Bagnold [START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF] (Section 1.2), we derive here the corresponding sediment flux at saturation q sat , i.e. in the steady and homogeneous case, for 67P. The saturated flux can generally be expressed as the product of the number N of transported grains per unit area and the mean grain horizontal velocity u p :

q sat = 1 φ b π 6 d 3 N u p , (5.16)
where φ b is the bed volume fraction [START_REF] Durán | Numerical simulation of turbulent sediment transport, from bed load to saltation[END_REF]. The basal shear stress τ = ρu 2 * is decomposed into the sum of the grain-born and fluid-born contributions τ p + τ f . The grain-born stress is τ p = N F drag , with F drag the drag force acting on a grain moving at the average velocity u p . The fluid-born stress must be the threshold stress τ t ≡ ρu 2 t at equilibrium transport. We then obtain N = (ττ t )/F drag .

(5.17)

In the terrestrial sub-aqueous bed-load, because the density ratio ρ p /ρ is on the order of a few units, the drag length is equal to a few d. The moving grains then quickly reach a velocity comparable that of the fluid u. In the cometary case, however, this drag length is much larger than the comet size, so that u p remains much smaller than u. This gives an almost constant drag force F drag on the moving grains, equal to that when the grains are static. This situation of constant mechanical forcing then resembles, for the thin transport layer, a granular avalanche, in which dissipation comes from the collisions between the grains and is thus increasing with u p [START_REF] Andreotti | A mean-field model for the rheology and the dynamical phase transitions in the flow of granular matter[END_REF]. In that case, it has been shown that, close enough to the threshold, the grain velocity follows the scaling law

u p ∼ gd, ( 5.18)
with a multiplicative factor around unity [START_REF] Daerr | Dynamical equilibrium of avalanches on a rough plane[END_REF].

Combining Eq. 5.16, 5.17 and 5.18, these expressions give:

q sat ∼ π 6φ b τ -τ t F drag gd. ( 5.19)
For τ on the order of a few τ t , the number of moving grains per unit surface soon reaches N ≃ 1/d 2 , which means that all the grains of this surface transport layer move. Leading to a typical flux on the order of

q sat ≃ g 1/2 d 3/2 , ( 5.20)
which gives q sat at the order of 10 -6 ∼ 10 -5 m 2 /s for the centimeter scale grains on 67P.

Saturation length L sat

The saturation length L sat has been well modeled and measured in saltation and suspension (Eq. 1.4 and 1.5). Regarding the case of bedload transport, it has never been directly measured. As in [START_REF] Fourrière | Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by non-linear pattern coarsening[END_REF], we calibrate here the behavior of the saturation length for bedload, by deducing L sat from controlled experimental measurement of the wavelength of emerging sub-aqueous ripples for various situations, corresponding to the fastest growing mode

A brief summary

In this chapter, we have discussed the sediment transport on 67P. In the first place, the transport threshold is evaluated and compared to the vapour flow along the surface. The comparison has confirmed the existence of sediment transport on 67P, near perihelion.

Then, we derived the saturated flux q sat and saturated length L sat for the cometary case.

We retain these laws for the traction on the comet and study the nature of the bedforms in next chapter.

Chapter 6

The nature of the bedforms

In the previous Chapter 4 and 5, we have modelled the atmosphere, and derived the transport laws on 67P case by making the observed bedforms an analogue to the terrestrial bedload patterns, which simply form by linear instability [START_REF] Charru | Sand ripples and dunes[END_REF]. Based on these results, we employ in this chapter the linear stability analysis to study the emergence and evolution of the bedforms on 67P .

Dispersion relation

We recall here the dispersion relation (Eq. 1.6 and 1.7) for the linear stability analysis of a flat granular bed and make it dimensionless by asking k = kν/u * , and Lsat = L sat u * /ν:

σ = k2 ( -) -kL sat 1 + ( kL sat ) 2 , ( 6.1)
c = k + ( -) kL sat 1 + ( kL sat ) 2 . (6.2) with σ = σ ν u * 2
and c = c ν u * , respectively, the rescaled growth rate and rescaled propagation speed. For , which quantifies the sediment transport, we take for its value the scaling law (5.20) in the cometary case:

≈ g 1/2 d 3/2 ; (6.
3)

The slope effect = 1 µ τ t /τ is estimated by the cometary parameters, with µ = tan(29 • ) ≃ 0.55. and are determined by discussing the basal shear stress on an undulated bed in the following.

The nature of the bedforms

The shear stress exerted by a flow in the x-direction on a fixed granular bed of elevation z = Z(x) can be computed by means of hydrodynamic equations presented in Section 1.4. When the bed is modulated as Z(x) = ζe ikx , these equations can be linearized with respect to the small parameter kζ and solved for non-slip conditions on the bed and vanishing first order corrections at z → ∞. The shear stress takes the generic form τ xz = ρu 2 * 1 + kζe ik x S t where S t is a dimensionless function of the rescaled vertical coordinate kz. and are defined as

S t (0) = + i , (6.4)
and they are functions of kν/u * , which is the inverse of the Reynolds number based on the wavelength, and can be interpreted as a Reynolds number for the perturbation. and in the smooth hydrodynamic regime has been discussed in [START_REF] Charru | Sand ripples and dunes[END_REF]. Depending on kν/u * , three asymptotic regimes can be identified. At large wavenumbers, the disturbed pressure gradient is balanced by inertia and by the viscous stress, and asymptotic expressions of and have been derived in the viscous and inertial asymptotic regime:

+ i = 2 + i 1 2 (kδ i ) -3 (6.5)
+ i = γ L (kδ i ) -1 e iπ/6 (6.6)

where γ L ≃ 1.06 is a constant [START_REF] Benjamin | Shearing flow over a wavy boundary[END_REF][START_REF] Charru | phase diagram' of interfacial instabilities in a two-layer couette flow and mechanism for the long-wave instability[END_REF], and δ i ∼ (ν 2 /ku 2 * ) 1/3 represents the penetration depth of vorticity disturbances. For small wavenumbers, turbulent Reynolds stress balances the disturbed pressure gradient. In the turbulent asymptotic regime, and are deduced as:

+ i = 2 U 2 (δ m ) U 2 (δ i ) 1 + 1 + ln(π/2) + 4γ E + iπ ln(δ i /z 0 ) , (6.7)
where U(z) is the logarithmic velocity profile, γ E ≃ 0.577 is Euler's constant and δ m /z 0 = (π/(2kz 0 )) ln -1/2 (π/(2kz 0 )) [START_REF] Kroy | Minimal model for aeolian sand dunes[END_REF][START_REF] Weng | Air flow and sand transport over sand dunes[END_REF].

With the wind and granular parameters on 67P estimated in the previous chapters, and are calculated over the whole range of the wave number, using Hanratty's model for the turbulent closure (black solid lines in Fig. 6.1). Apparently, the laminar regime is separated from the turbulent regime by a transitional regime where a 'crisis' can be observed, which nicely fits the experimental measurements [START_REF] Charru | Sand ripples and dunes[END_REF]. In this regime, the perturbation partly penetrates the turbulent region so that both viscosity effects and turbulent fluctuations are needed to be taken into account. The flow is accelerated/deaccelerated on the upstream/downstream face of a bump, where the streamwise pressure gradient tends to damp the turbulent fluctuations so that the viscous sublayer thickens, and therefore there exists a periodic transition from the viscous to the turbulent inner layer. However, it

is still an open question for a true understanding of the interplay between a wavy bottom, even though Hanratty's model provides a convenient parameterization.

Most unstable mode

In this section, we select the fastest growing mode, i.e., most unstable mode, to predict the temporal and spatial characteristics of these patterns.

Wave length selection

A most unstable mode km corresponding to the maximum growth rate is deduced from Eq. 6.1 as the solution of d σ d k = 0. The corresponding growth rate is σm = σ( km ), and the propagation speed is cm = c( km ).

We compute the most unstable mode, and get a relation between k and Lsat . Using

L sat /λ = Lsat k/(2π), one then obtains λ L sat = f L sat u * ν , (6.8)
as in Fig. 6.2, along with a phenomenological fit Empirically, we find that it approximatively scales as λ/L sat ≈ (L sat u * /ν) -2/5 . With Eq 5.21 and 6.8, we compute the most unstable wavelength as a function of grain size d. As in Fig. 6.3, λ ranges from 10 to 20 m for a mean grain diameter d between 10 and 40 mm, in good agreement with the observed crest-to-crest distance, as in Table . 6.1.

Bedform growth and propagation

For grain size d between 10 and 40 mm, the traction sediment flux is on the order of 4 10 -5 m 2 /s. The corresponding ripple growth time, deduced from the linear stability analysis is Fig. 6.1 Basal shear stress components in phase (a) and in quadrature (b) with respect to the bed elevation, as functions of the rescaled wave number kν/u * . Three asymptotics are identified: turbulent regime (blue dashed line, Eq. 6.7), inertia regime (green dashed line, Eq. 6.6) and viscous regime (red dashed line, Eq. 6.5); in addition, a transitional region where a 'crisis' can be observed. Accordingly, comparing pictures before and after perihelion (Fig. 3.2a), we observe that some ripples have disappeared at the downwind end of the field and others have nucleated at the upwind entrance. In between, ripples have propagated. To get the propagation length, the photographs are mapped one on the other using fixed elements of relief (cliffs, rocks, holes, etc) that can be recognized on both pictures. The mapping is performed through a projection, assuming in first approximation that the landscape is planar. The ripple crests and avalanche slip faces can be easily detected and can be transposed from one photograph to the other, and it is estimated on the order of 10 m over one revolution.

Mature ripples at a wavelength of 18 m (kν/u * ≃ 0.024) propagate at a velocity:

c ≃ 0.18 u * /ν ≃ 10 -5 m/s. (6.10)
Therefore, the propagation length is around ≃ 10 m for the time during sediment transport occurs ≃ 10 6 s, this is also consistent with the observations (Fig. 3.2a). Fig. 6.3 Relation between the wavelength and the mean grain diameter predicted at perihelion, for afternoon conditions. The most unstable mode of the linear instability selects the emergent wavelength, which depends on the grain diameter through the saturation length L sat (Fig. 5.7). Yellow mark: range of measured crest-to-crest distance and grains size (Table . 6.1).

Conclusions

In this part, we have proposed an explanation for the bedforms observed on 67P. The discussions are mainly three parts, namely thermo-hydrodynamics of atmosphere, sediment transport and the nature of the bedforms.

In modelling the dynamics of the atmosphere, we find that the existence of the ripples on 67P changes the understanding of outgassing on the comet surface. If the vapor emits with a velocity like the spectacular jets, it feeds an extremely thin atmosphere, which is unlikely to set the grains on 67P's surface into motion. We thus suggest the presence of a porous surface granular layer that slows down the radial outgassing velocity and results in a thicker atmosphere. The pressure gradient due to the day-night alternace results in the surface vapor flow. It is verified to be responsible for the sediment transport of centimeter scale grains on the comet surface, by comparing to the transport threshold.

Considering the low escape velocity and a thick viscous sublayer flow on 67P's surface, we have proved that the transport mode is only traction. Making an analogue to the terrestrial subaqueous ripples, the transport laws are derived using Bagnold's idea, and controlled experimental measurements. Finally, the temporal and spatial scales are predicted for the cometary case by selecting the fastest growing mode in the linear instability analysis, and the predicted results are consistent with the observations. These bedforms live in the viscous sublayer of 67P's atmosphere, they are thus 'viscous ripples'. Previous and recent investigations have shown that surface pressure of the atmosphere of Pluto is about p 0 = 1 Pa, roughly 100,000 times less than Earth's atmospheric pressure [START_REF] Gladstone | The atmosphere of pluto as observed by new horizons[END_REF][START_REF] Schenk | New horizons at pluto[END_REF]. The surface temprearture is around T 0 = 40 K, but it quickly grows with altitude due to methane-generated greenhouse effect (Fig. 8.2). Pluto's atmosphere is thus stably stratified. Near the altitude 30 km it reaches 110 K, and then slowly decreases [START_REF] Stern | Pluto' in Encyclopedia of the Solar System[END_REF][START_REF] Gurwell | Sub-arcsecond scale imaging of the pluto/charon binary system at 1.4 mm[END_REF][START_REF] Dias-Oliveira | Pluto's atmosphere from stellar occultations in 2012 and 2013[END_REF]. We could estimate the saturated state of the atmosphere on the surface of Pluto. The saturation vapour pressure p sat is a calibrated function of the solid-vapour interface temperature [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF], evaluated in bars by Excluding the sublimation latent heat in Eq. 8.2 actually means a saturated atmosphere near Pluto's surface. Thus, we can find the pressure and density with Eq. 8.1 and ideal gas law. The results show that the gentle variation in temperature also gives the similar phenomenon in the pressure (Fig. 8.4) and density (Fig. 8.5). To understand the gentle variations in these parameters, Eq. 8.2 is linearised and transformed into the Fourier space:

ln p sat = A 0 + j A j /T j (8.
(1 -Ω N ) ψ(1) = ΣεT 3 (0) T(1) + Ĵ(1) , (8.3)
where the subscript (0) and ( 1) indicate the base state and its first order correction.

Considering

T (0) ∼ [ψ/(Σε)] 1/4 and Ĵ(1) ∼ k s |ω| 2κ s T(1)
, one compares the heat flux and the radiation term, which gives the dimensionless number controlling the day-night temperature difference:

T = k 2 s ρ 2 s C 2 s Σε Γ 2 d ψ 3 1/4 . (8.4)
Plugging into the values, we have P T = 2 10 7 for Pluto and C T = 3 10 3 for comet 67P. This number is large for both comet 67P and Pluto, which means most of the first order correction of solar radiation contributes into the heat flux, rather than modulating the temperature. Moreover, P T is larger than C T over four orders of magnitude, and this is why the day-night variation on Pluto is much smaller than that for comet 67P (Fig. 4.2).

Hydrodynamical description

We consider an atmospheric boundary layer to compute the wind flow close to Pluto's surface. As the situation for comet 67P (Section 4.3.2), we assume the same velocity profile

u θ (r) = u * κ ln 1 + r-R p z 0
and dimensionless number Λ ≡ ln 1 + δ i z 0 as in Eqs. 4.37 and 4.38. Neglecting the horizontal diffusion, the momentum balance in the boundary layer approximation writes in Eq. 4.36. Integrating Eq. 4.36 from the ground (r = R p ) to the at the top of the inner turbulent boundary layer (r = R p + δ i ), where we assume a vanishing leak, we obtain:

ρ 0 |u * |u * + d dθ 2 -2Λ + Λ 2 δ i κ 2 R p ρ 0 u 2 * = - δ i R p dp 0 dθ , (8.5)
where we have used the fact that the velocity u θ vanishes at the Pluto's surface, and that the shear stress and u r vanish at r = R p + δ i . The inertial terms are comparable to the pressure gradient at the crossover altitude of the boundary layer thickness, one thus has:

ρ 0 |u * |u * ≈ d dθ 2 -2Λ + Λ 2 δ i κ 2 R p ρ 0 u 2 * , (8.6)
and Eq. 8.5 simplifies into:

- δ i 2R p dp 0 dθ = ρ 0 |u * |u * . (8.7)
We make here the same assumption as the comet 67P case that variations of all quantities along θ are slow, which allows us to take d dθ ≈ 1 2π , and thus Eq. 8.6 simplifies into

2 -2Λ + Λ 2 δ i 2πκ 2 R p ≈ 1. (8.8)
With the density ρ 0 and pressure p 0 profiles from Section 8.2.1, we finally solve (8.7) and (8.8) to obtain u * (Fig. 8.6). The results suggest the daily peak values of u * are almost invariant at the revolution scale, and direction of u * shows periodic variation, resulted from the fact that the directions of pressure gradient are opposite before and after noon.

The computed δ i is a around 10 kilometers, much larger than the thickness of the viscous sublayer of a few meters.

A brief summary

In this chapter, we have simply considered Pluto's atmosphere. The results show that thermodynamical parameters (T 0 , ρ 0 and p 0 ) show much smaller variation during the day-night alternance, compared to those of comet 67P; so are they at the revolution scale.

Regarding the hydrodynamics, the computed u * also shows gentle variation at both rotation and revolution scales, in spite of the periodic variation in direction.

In the present model, we have neglected the ice sublimation when considering the energy balance at Pluto's surface, which allows us to consider the thermal processes and hydrodynamics separately. The computed temperature and pressure are close to the measured data [START_REF] Gladstone | The atmosphere of pluto as observed by new horizons[END_REF], and this confirms this simplification at a certain extent. Admittedly, more detailed model is needed to consider the ice sublimation, which couples thermal and hydrodynamical calculations, yet beyond the discussion in this thesis.

Governing equations

Thermal processes

We consider the thermal processes separately in the solid and in the vapour, and couple them later on at the interface. In the solid, we write the heat flux as ⃗ j = -k s ⃗ ∇T , where T is the temperature and k s is the thermal conductivity. Introducing the bulk mass density of the solid ρ s and its heat capacity C s , the heat conservation equation reads:

ρ s C s ∂ t T = k s ∇ 2 T. (9.1)
All three parameters k s , C s and ρ s are assumed to be homogeneous in space. Equivalently, a temperature diffusion equation can be written, which involves the thermal diffusivity

k s /(ρ s C s ).
In the vapour, the convective heat flux, ρC v T ⃗ u is not a galilean invariant and does not vanish when there is flow velocity. This flux, from the point of view of heat conservation, is ρC v (T -T 0 )⃗ u; and thus we write the heat flux as the sum of a convective term and a diffusive term:

⃗ j = ρC v (T -T 0 )⃗ u -k v ⃗ ∇T (9.2)
where T, ρ, C v and k v are respectively the temperature, the mass density, the heat capacity and the thermal conductivity of the vapour. Then, the heat conservation equation reads:

ρC v ∂ t T + ∂ x j x + ∂ z j z = 0. (9.
3)

The thermal conductivity k v will be specified in the following.

Hydrodynamics

For the vapour flow, we assume an incompressible turbulent situation described by Reynolds Averaged Navior-Stokes equations (Eqs. 1.8 and 1.9), with the Reynolds stress closed by the Prandtl mixing length (Eq. 1.15 ) and Hanratty's model (Eq. 1.16), as presented in Section 1.4. For the thermal description, it should be noted that k v composed of a turbulent and a viscous term:

k v = ρC v ℓ 2 |γ| β t + ν β ν , ( 9.4)
where β t and β ν are the turbulent and viscous Prandtl numbers. For an ideal gas, β ν = 1.

As for the turbulent Prandtl number, a typical value for gases is β t = 0.7.

Sublimation interface

The vapour sublimation/condensation occurs at the interface, and the flux q (in mass per unit surface and unit time) is governed by the Hertz-Knudsen law, and we recall here:

q = αV th (ρ sat -ρ) , (9.5)
where V th is the thermal velocity defined in Eq. 4.17 as the average velocity of the molecules. α is a dimensionless factor, which encodes the probability that a molecule of the vapour is absorbed or emitted at the surface of the solid. Its typical value is between 10 -2 and 1, and, for simplicity, we neglect its dependence with temperature. The saturation vapour density ρ sat is related to saturated pressure p sat by the ideal gas law and p sat is a calibrated function of the temperature as in Eq. 8.1 [START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF], here evaluated at the solid-vapour interface and denoted as T i . Later on, the vapour density and the pressure at the interface will similarly be denoted as ρ i and p i respectively.

To obtain q, we can expand the saturation density ρ sat , the density ρ and the thermal velocity V th around the reference state (T 0 , ρ sat (T 0 ), p sat (T 0 )) at the interface. As for ρ sat ,

ρ sat T i = ρ sat (T 0 ) + ρ ′ sat (T 0 ) T i -T 0 , (9.6)
where the prime means the derivative with respect to the temperature. Using the ideal gas law, we can express ρ ′ sat = ρ sat T T p ′ sat p sat -1 . Now using the Clausius-Clapeyron equation p ′ sat = ρ sat /T [START_REF] Wark | Generalized Thermodynamic Relationships[END_REF], we can write

ρ ′ sat (T 0 ) = ρ sat (T 0 ) T 0 m k B T 0 -1 , (9.7)
plugging in 9.6, it gives

ρ sat T i = ρ sat (T 0 ) + ρ sat (T 0 ) T 0 m k B T 0 (T i -T 0 ) - ρ sat (T 0 ) T 0 (T i -T 0 ), (9.8)
where is the latent heat of ice sublimation, and m is the mass of a vapor molecule.

Similarly for thermal velocity, one has V th (T i ) = V th (T 0 ) + V ′ th (T 0 ) T i -T 0 , and it further gives

V th T i = V th (T 0 ) 1 + T i -T 0 2T 0 . ( 9.9)
Regarding the density, one needs to do the expansions in both T i at T 0 and p i at p sat (T 0 ):

ρ i = mp i k B T i = m p i -p sat (T 0 ) + p sat (T 0 ) k b [(T i -T 0 ) + T 0 ] , (9.10
)

ρ i = ρ sat (T 0 ) + m k B T 0 p i -p sat (T 0 ) - ρ sat (T 0 ) T 0 (T i -T 0 ). (9.11)
From Eq 9.5, 9.8, 9.9 and 9.11, finally obtain:

q = αV th (T 0 )ρ sat (T 0 ) 1 + T i -T 0 2T 0 m k B T 0 (T i -T 0 ) T 0 - p i -p sat (T 0 ) p sat (T 0 ) . ( 9.12)
As we will consider small variations of T i around T 0 , one can neglect the variation associated with the thermal velocity so that:

q = αV th (T 0 )ρ sat (T 0 ) m k B T 0 (T i -T 0 ) T 0 - p i -p sat (T 0 ) p sat (T 0 ) . ( 9.13)
At the interface, the vertical and horizontal velocities of the vapour are respectively related to vapour flux q and no-slip condition:

u z (ξ + ) = q ρ i , (9.14)
u x (ξ + ) = 0. (9.15)
For the thermal interfacial conditions, the temperature must be continuous, but the heat flux across it is not: to sublimate ice at a rate corresponding to a vapour flux q, an energy flux q is needed. This quantity must be provided by the heat flux difference. These two conditions at the interface elevation z = ξ (later in Eq. 9.25) thus write:

T (ξ -) = T (ξ + ), (9.16
)

j(ξ -) = q + ΣεT 4 i -j(ξ + ) -j ψ , (9.17)
where ΣεT 4 i is the surface radiation, with Σ and ε the Stefan constant and emissivity, and j ψ is a light power flux, as the ice surface is submitted to an insolation. j ψ will be specified in the following.

Self-illumination on a modulated surface

In this subsection, we determine the illumination of a modulated surface z = ξ(x) in a way similar to [START_REF] Claudin | Physical processes causing the formation of penitentes[END_REF][START_REF] Betterton | Theory of structure formation in snowfields motivated by penitentes, suncups, and dirt cones[END_REF]. Due to a finite albedo Ω N , a unit surface re-emits a light power flux Ω N j ψ , proportional to the power received j ψ . Note that j ψ is (1 -Ω N ) times the power emitted by the sun. When the interface is flat, none of the re-emitted photons reach the surface again. However, when the surface is modulated, its illumination is partly due to these photons. Assuming isotropy of the re-emission, and a one-dimensional profile, one obtains:

j ψ (x, ξ) = j 0 ψ + Ω N x b x a | | π j ψ (X , ξ(X )) dX , (9.18)
where j 0 ψ is the power received by a flat solid, due to sun illumination, and dξ is the solid angle through which the element dξ at position X is seen from point x, which reads

(X ) = 1 x -X ξ ′ (x) - ξ(x) -ξ(X ) x -X . ( 9.19)
In this expression, ξ ′ is the derivative of the interface profile. The bounds x a and x b of this integral both depend on x too. They correspond to positions beyond which the interface cannot be seen from position x, due to shadowing (Fig. 9.2). As the reference state considered is homogeneous, at the linear order, eigen-modes of the illumination operator (9.18) are periodic. However, due to the non-local nature of , they are not Fourier modes, as known in the general context of Fredholm equations. In particular, the illumination of a sinusoidal profile ξ(x) = H cos(kx) is not strictly sinusoidal. For such a function, the contribution to the integral term giving the illumination at first order in kH reads:

I(η x ) = η b η a sin η x + cos η x -cos η ′ x η x -η ′ x dη ′ x |η x -η ′ x | , ( 9.20)
where η x = kx. The boundaries of the integral, η a (η x ) and η b (η x ), correspond to rays that are tangent to the surface (Fig. 9.2). For 0 < η x < π, they are solutions of For π < η x < 2π, the bounds are obtained by symmetry. The dependence of these bounds on η x is displayed in Fig. 9.3a. However, the non-harmonic contribution of the modes turn out to be negligible and the integral (9.20) is numerically found to be very close to the function 1 -cos(kx) (Fig. 9.3b). For the linear stability performed here, the light volumetric power at the interface can be approximately written as

sin η x + cos η a -cos η x η a -η x =0 for 0 ≤ η x ≤ π/2 (9.21) sin η a + cos η x -cos η a η x -η a =0 for π/2 ≤ η x ≤ π (9.22) sin η b + cos η b -cos η x η b -η x =0 (9.23)
j ψ (x) = j 0 ψ [1 + ΩkH (1 -cos(kx))] , (9.24)
with Ω = Ω N /π.

Linearised problem

In this section, we perform the linear analysis of the problem. For small enough amplitudes, we can consider a bottom profile (Fig. 9.1) of the form

ξ(x) = H e σt+ikx ( 9.25)
without loss of generality (real parts of expressions are understood). λ = 2π/k is the wavelength of the bottom and H the amplitude of the corrugation. The case of an arbitrary relief can be deduced by a simple superposition of Fourier modes. σ is the growth rate of the elevation. We make the usual assumption that the time 1/σ is much larger than all time scales involved in the hydrodynamical and thermal processes, so that the solid surface can be considered as fixed in the computation of the velocity, stress and temperature fields.

Base state

The actual base state is inhomogeneous and time dependent. There is a partition of the energy at the surface between latent heat, conductive flux in the solid, conductive and convective flux in the atmosphere, illumination by the sun and radiation. Furthermore, gradients lead to horizontal transport of energy. For the sake of simplicity, we consider a flat solid surface (Dashed line in Fig. 9.1) which presents a uniform temperature T 0 , and thus the surface temperature T i = T 0 ; and we assume that the thermal flux is both homogeneous in the solid and in the atmosphere, which are respectively noted j 0 -and j 0 + . Moreover we consider that the sunlight flux above the surface is uniform, equal to j 0 ψ , and induces a constant migration of the surface at a velocity -q 0 /ρ s . The pressure p = p 0 above the surface is also homogeneous. It relates to the vapour density as ρ 0 = ρ sat (T 0 ) + m k B T 0 (p 0 -p sat (T 0 )).

Following (9.13), the sublimation/condensation flux in the base state obeys the equation:

q 0 = -αV th (T 0 )ρ sat (T 0 ) p 0 -p sat (T 0 ) p sat (T 0 ) . ( 9.26)
q 0 can be positive (net ablation) or negative (net accumulation), depending on whether the pressure p 0 is below or above its saturated value p sat (T 0 ). As in Eqs. 9.14 and 9.17, the corresponding vertical velocity of the vapour and the power balance at the interface read:

w 0 ≃ q 0 ρ 0 , (9.27)
j 0 -= q 0 + ΣεT 4 0 -j 0 + -j 0 ψ . (9.28)
w 0 is positive when sublimation occurs.

Regarding the vapor flow, a steady and homogeneous situation along the x-axis is considered. The strain rate reduces to ∂ z u x . The x component of Eq. 1.9 reduces to

ρu z ∂ z u x = ∂ z τ xz .
Because in the base state u z = w 0 and ρ = ρ 0 are independent of z, this balance can be written as

∂ z (τ xz -ρ 0 w 0 u x ) = 0, or equivalently τ xz -ρ 0 w 0 u x ≡ ρ 0 |u * |u * .
Summing up the turbulent and viscous contributions of the shear stress τ xz , we obtain:

ρ 0 ℓ 2 |∂ z u x |∂ z u x + ν∂ z u x -ρ 0 w 0 u x = ρ 0 |u * |u * . (9.29)
Similarly, Eq. 4.9 reduces to ρ 0 C v w 0 Tk v ∂ z T = j 0 + . With the expression in Eq. 9.82 for k v , it reads:

ρ 0 C v w 0 (T -T 0 ) -ρ 0 C v ℓ 2 |γ| β t + ν β ν ∂ z T = j 0 + . (9.30)
In the solid, the temperature field obeys Fourier's law as -k s ∂ z T = j 0 -and it easily solves into:

T = T 0 - j 0 - k s z , (9.31)
where we have used the boundary condition T = T 0 at the interface.

First order fields

With respect to the small parameter kH, we can write all the relevant quantities f above the solid surface under the generic form:

f (x, z) = f (η) + kH e σt+ikx F (η), (9.32)
where f corresponds to the homogeneous (base) state, and F encodes the profile of the first order correction. More explicitly, for all fields in the atmosphere, we define non-dimensional functions U, W , ..., Y of the dimensionless variable η = kz by:

u x = u * + kH e σt+ikx U , (9.33
)

u z = u * + kH e σt+ikx W , (9.34
)

τ xz = τ z x = ρ 0 u 2 * 1 + kH e σt+ikx S t , (9.35)
τ zz -p = -p 0 -ρ 0 u 2 * kH e σt+ikx S n , (9.36)
p = p 0 + ρ 0 u 2 * kH e σt+ikx P, (9.37
)

τ zz = ρ 0 u 2 * kH e σt+ikx S zz , (9.38)
τ x x = ρ 0 u 2 * kH e σt+ikx S x x , (9.39
)

kℓ = Υ + kH e σt+ikx L, (9.40)
T = T 0 + u 2 * T 0 -+ kH e σt+ikx Θ , (9.41
)

j z = j 0 + -ρ 0 C v u 3 * T 0 kH e σt+ikx Y, (9.42)
j ψ = j 0 ψ 1 -ΩkH e σt+ikx . (9.43)
With these perturbed fields, we will develop the linearised governing equations later on.

Base state in a dimensionless form

Starting from the base state, we now make the normalized and identity the dimensionless numbers. The base state flow field reads:

u x = u * , and u z = w 0 = u * , (9.44)
and we have the dimensionless mixing length Υ = kℓ (Eq. 9.40). From Eq. 1.15, it is rescaled as:

Υ (η) = κη 1 -exp(-η/R 0 t ) , (9.45)
where we have introduced the Reynolds number with respect to the wavelength:

= u * kν . ( 9.46)
Plugging Eq. 9.44 and 9.45 in 9.29, the function (η) giving the wind profile can be computed by:

Υ 2 | ′ | ′ + -1 ′ - = 1, (9.47)
or equivalently

′ = -1 + 1 + 4Υ 2 2 (1 + ) 2Υ 2 , ( 9.48)
which must be solved with the boundary condition (0) = 0 corresponding to the no-slip condition of the wind at the solid interface (Eq. 9.15). When Υ 2 2 ≪ 1, Eq. 9.48 simplifies into:

′ = (1 + ), (9.49)
and it further gives ′ (0) = . The rescaled vertical velocity reads:

= w 0 u * = q 0 u * ρ 0 . ( 9.50)
= 0 corresponds to a reference state that is not sublimating.

The base state temperature field, strain rate modulus and thermal conductive coefficient are

T =T 0 - u 2 * T 0 , (9.51)
|γ| = u * k ′ , (9.52
)

k v = ρC v u * k ′ Υ 2 β t + -1 β ν . ( 9.53)
Plugging Eq. 9.51-9.53 in 9.30, the temperature profile obeys:

Υ 2 | ′ | β t + -1 β ν ′ - = j 0 + ρ 0 C v T 0 u 3 * = + 3 , (9.54)
where we have introduced the dimensionless quantities

+ = j 0 + k 3 B T 5 0 C 2 v α 3 m 3 5 V 3 th ρ 0 and = k B T 2 0 C v u * α 2 mV th , (9.55)
we will specify later how heat flux j + is rescaled to get + . It is obvious that sublimation (net ablation) conditions correspond to negative values of + . Eq. 9.54 must be solved with the boundary condition that the temperature is T 0 at the interface, i.e. (0) = 0, which gives:

′ (0) = + β ν 3 -1 = + β ν 2 , (9.56)
with = -1 .

Interfacial equations

From Eq. 9.41, the temperature at the interface z = ξ reads:

T i -T 0 = u 2 * T 0 -(ξ) + kH e σt+ikx Θ(ξ) . (9.57)
Considering (ξ) = (0) + ′ (0)kξ and Θ(ξ + ) = Θ(0) + Θ ′ (0)kξ at the linear order, Eq. 9.57 becomes

T i -T 0 = u 2 * T 0 kH Θ(0) -′ (0) e σt+ikx , (9.58)
where we have used (0) = 0 and neglected a second order term. According to the temperature continuity at the interface (9.16), one has

T (ξ -)-T 0 + j 0 - k s ξ -= T (ξ + )-T 0 + j 0 - k s ξ -. Introducing 9.
T (ξ -) -T 0 + j 0 - k s ξ -= u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx .
In the solid, the temperature field obeys the Laplace equation, and can therefore be expressed as:

T -T 0 + j 0 - k s z = u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx+kz . ( 9.59)
From Fourier's law of heat conduction, the thermal fluxes in the solid are:

j x = -ik s k u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx+kz , (9.60)
j z -j 0 -= -k s k u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx+kz . ( 9.61)
The heat flux at the interface reads (Eq. 9.42):

j z (ξ + ) -j 0 + = -ρ 0 C v u * u 2 * T 0 kH e σt+ikx Y (0), (9.62)
where we have used Y (ξ + ) = Y (0) + Y ′ (0)kξ and neglected the second order terms. On the other hand, Eq. 9.61 gives

j z (ξ -) -j 0 -= -k s k u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx . ( 9.63)
Recalling 9.17 and 9.28, one obtains the conditions at the interface:

ρ 0 C v u * u 2 * T 0 kH e σt+ik x Y (0) = k s k u 2 * k T 0 Θ(0) -′ (0) + j 0 - k s H e σt+ikx + (q -q 0 ) + 4 u 2 * ΣεT 4 0 Θ(0) -′ (0) kH e σt+ikx + j 0 ψ ΩkH e σt+ikx . ( 9.64)
where we have used Eq. 9.43 and linearised the surface radiation term Σε(T 4 i -T 4 0 ). Furthermore, it gives

q -q 0 = ρ 0 C v u 3 * T 0 2 Y (0) - u 2 * k s k 2 T 0 Θ(0) -′ (0) - j 0 - -4 u 2 * 2 ΣεT 4 0 Θ(0) -′ (0) - j 0 ψ Ω kH e σt+ikx . ( 9.65)
The growth rate σ is estimated by the mass conservation as:

q -q 0 ρ s = -σH e σt+ikx . ( 9.66)
From (9.13) and (9.26), we obtain:

q -q 0 ρ s = αmV th (T 0 ) k B T 0 ρ s ρ sat (T 0 ) T 0 T i -T 0 -p i -p 0 (9.67)
where we have taken the fact that ρ 0 ≃ ρ sat (T 0 ). Introducing Eqs. 9.58 and 9.69, Eq. 9.67 becomes:

q -q 0 ρ s = α mV th ρ s k B T 0 ρ 0 u 2 * Θ(0) -Sn(0) -2i - + β ν 2 kH e σt+ikx , (9.68)
where we have used Eq. 9.56, 9.90, and the interfacial pressure:

p i = p 0 + ρ 0 u 2 * kH e σt+ikx P(0). ( 9.69)
Recalling Eq. 9.66, we therefore get the complex growth rate:

σ = f σ -1 S n (0) + 2i -Θ(0) + + β ν 2 , (9.70)
whose real and imaginary parts give the temporal growth rate and the propagation speed respectively, and with f σ = αmV th ρu 3 * ρ s k B T 0 ν . Note that we have the normalized growth rate: σ = σ/ f σ . Another equation relating the flux to the velocity is:

q -q 0 = ρ 0 u * kH e σt+ikx W (0), (9.71)
which gives, recalling Eq. 9.65:

W (0) = Θ(0) -S n (0) -2i - + β ν 2 , (9.72)
where we have defined the dimensionless number:

= α 2 m 2 2 V 2 th k 2 B T 3 0 C v . ( 9.73)
Eliminating qq 0 from Eqs. 9.65 and 9.68, we obtain the equation relating the flux to the temperature:

αmV th k B T 0 Θ(0) -′ (0) -S n (0) -2i = C v u * T 0 2 Y (0) - k s T 0 k ρ 0 2 Θ(0) -′ (0) + j 0 - ρ 0 u 2 * - 4ΣεT 4 0 ρ 0 2 Θ(0) -′ (0) - j 0 ψ Ω ρ 0 u 2 * . ( 9.74)
Solving for Θ(0), it gives the lower boundary condition in a dimensionless form:

Θ(0) = + β ν 2 + 1 1 + + 4 Σ Y (0) + S n (0) + 2i - -+ Ω ψ 2 , (9.75)
where we have used q 0 = ρ 0 u * , and rescaled all the heat flux terms by a reference flux

α 3 m 3 5 V 3 th ρ 0 k 3 B T 5 0 C 2 v : -= j 0 -k 3 B T 5 0 C 2 v α 3 m 3 5 V 3 th ρ 0 , Σ = Σεk 3 B T 9 0 C 2 v α 3 m 3 5 V 3 th ρ 0 , ψ = j 0 ψ k 3 B T 5 0 C 2 v α 3 m 3 5 V 3 th ρ 0 and = , (9.76) with
= k s ρ 0 C v ν , = α 2 m 2 3 V 2 th k 2 B T 4 0 C 2 v and = k 2 B T 3 0 C v α 2 m 2 2 V 2 th .
(9.77) characterizes the heat diffusion in the solid. With these rescaled fluxes, one obtains the normalised form for Eq. 9.28,

-= + Σ -+ -ψ . (9.78)
We could consider a reduced model by putting Y (0) = 0 and S n (0) = 0 in Eq 9.75. Plugging into Eq 9.70, one gets the reduced growth rate:

σ r e = f σ -+ Ω ψ 3 (1 + + 4 Σ ) , (9.79)
and σre = σ r e / f σ . This analytical scaling serves a test on the full solution at the limit of neglecting the heat flux and pressure effect.

Linearised system

In this section, we linearise the governing equations and the boundary conditions to get the mathematical description for the system. For starters, the linearised strain rate modulus, viscous and turbulent stress coefficient and the thermal conductivity are computed at the linear order:

|γ| = u * k ′ + (U ′ + iW )kH e σt+ikx , (9.80)
l 2 |γ| + ν = u * k ′ Υ 2 + + (U ′ + iW )Υ 2 + 2Υ L ′ kH e σt+ikx , (9.81)
k v = ρC v u * k ′ Υ 2 β t + -1 β ν + 1 β t (U ′ + iW )Υ 2 + 2Υ L ′ kH e σt+ikx . (9.82)
These formulas are of use in the following linearisation.

Linearised equations

We first linearise the Hanratty equation (Eq. 1.16), and it gives at the linear order:

( + ia) = i(S x x -S zz -S n ) = - 4U ′ -iS n . (9.83)
Now, we can express the disturbance to the mixing length:

L = κ -1 + exp(-η/R t ) 1 - η R 0 t + η 2 R 0 t (S t /2 -b ) , (9.84)
where b is the relative variation of R t due to the pressure gradient:

b = 1 R 0 t dR t d > 0. (9.85)
Plugging Eq 9.81, 9.33 -9.34, and 1.10 -1.12 into 9.35 -9.39, one obtains, at the linear order, the stress functions as follows:

S t = -1 + 2Υ 2 ′ (U ′ + iW ) + 2Υ ′2 L, ( 9.86)
S x x = 2 -1 + Υ 2 ′ iU = 2iU ′ , (9.87)
S zz = -2 -1 + Υ 2 ′ iU = -2iU ′ , (9.88)
P = S n + S zz = S n - 2iU ′ . (9.89)
where we have used -1 + Υ 2 ′ = 1/ ′ at the zeroth order. Considering Eq. 9.106, it gives from Eq. 9.89:

P(0) = S n (0) + 2i (9.90)
Similarly, we plug Eqs. 9.33 -9.39 into Navier-Stokes equations (Eqs. 1.8 and 1.9), and the linear expansion gives rise to:

W ′ = -iU, (9.91)
S ′ t = i U + ′ W + U ′ + iS n -iS x x + iS zz , (9.92)
S ′ n = -i W -W ′ + iS t . (9.93)
where we have introduced Eq. 9.89 when finding 9.92.

As for the scalar equations, we plug Eqs. 9.82, 9.41 and 9.42 into Eq. 9.2 -9.3, it gives

Y ′ = -i U + i + Υ 2 | ′ | β t + -1 β ν Θ, ( 9.94)
Y = Υ 2 | ′ | β t + -1 β ν Θ ′ - 1 β t ′ Υ 2 (U ′ + iW ) + 2Υ L ′ + W -Θ. ( 9.95)
With U ′ and Θ ′ expressed from Eq. 9.86 and 9.95, we then obtain six closed equations:

U ′ = -iW + S t -2Υ ′2 L -1 + 2Υ 2 ′ , (9.96)
W ′ = -iU, ( 9.97)
S ′ t = i + 4 ′ U + U ′ + ′ W + iS n , (9.98)
S ′ n = -i W -W ′ + iS t , (9.99)
Θ ′ = Y + Θ -W + ′ Υ 2 S t + 2Υ L ′ -1 + Υ 2 ′ β t ( -1 + 2Υ 2 ′ ) / Υ 2 | ′ | β t + -1 β ν , (9.100)
Y ′ = -i U + i + Υ 2 | ′ | β t + -1 β ν Θ, (9.101)
where

L = κ -1 + exp(-η/R t ) 1 - η R 0 t + η 2 R 0 t 2 S t + b 1 + ia -1 4U ′ + iS n . (9.102)

Boundary conditions

Boundary conditions are needed to solve the system. The upper boundary corresponds to the limit η → ∞, in which the vertical fluxes of mass and momentum vanish asymptotically.

This means that the first order corrections to the shear stress and to the vertical velocity must tend to zero: W (∞) = 0 and S t (∞) = 0. In practice, a boundary at finite height D (at η D = kD) is introduced, at which we impose a null vertical velocity and a constant tangential stress ρ 0 u 2 * so that:

W (η D ) = 0, (9.103)
S t (η D ) = 0. ( 9.104)
This corresponds to a physical situation where the fluid is entrained by a moving upper plate, for instance a stress-controlled Couette annular cell. Then, we consider the limit D → +∞, i.e. when the results become independent of D.

On the interface z = ξ, the horizontal velocity component must vanish u x = 0, and we have, from Eq. 9.33, (ξ) + U(ξ)kξ = 0; (9.105) considering (ξ) = (0) + ′ (0)kξ and U(ξ) = U(0) + U ′ (0)kξ, Eq. 9.105 gives:

U(0) = -′ (0) = - (9.106)
where we have used (0) = 0, and (0) = from Eq. 9.49. We recall the conditions for vertical velocity and temperature at the interface here:

W (0) = Θ(0) -S n (0) -2i - + β ν 2 , (9.107
)

Θ(0) = + β ν 2 + 1 1 + + 4 Σ Y (0) + S n (0) + 2i - -+ Ω ψ 2 . (9.108)
Furthermore, the heat flux through the upper boundary remains equal to j 0 + , because the first order correction from the bed profile modulation becomes negligible when it is far from the bed. We therefore have:

Y (η D ) = 0 (9.109)
To sum up, six boundary conditions (Eq. 9.103, 9.104, 9.106, 9.107 -9.109) are specified to solve the linearised governing eqautions (Eq. 9.96-9.101).

Solving strategy

Introducing the vector X = (U, W, S t , S n , Θ, Y ), we finally get at the first order in kH the following compact form of the equation to integrate:

d dη X = X + S (9.110)
with the coefficient matrix and S the source term. The general solution of this equation is the linear superposition of all solutions of the homogeneous system (i.e. with S = 0), and a particular solution.

In practice, we solve the equations using a fourth-order Runge-Kutta scheme with a logarithmic step. Because we have some boundary conditions at the bottom of the domain yet the others are at the top, we firstly integrate twice over half of the domain, with one from the bottom to the center and the other from the top to the center, and we then match the two solutions at the center of the domain. To be specific:

1. We cross upward a first time from the bottom to the top to compute the value of and at the top of the domain;

2. We cross upward from the bottom to the center of the domain to find one solution:

Dispersion relation

In this section, we analyse the physical content of the model by studying the dependence of the dispersion relation on the different parameters. Firstly, the parameters are summarised and/or estimated for the problem of Pluto's sublimation patterns. We take the measured values for the temperature and pressure at Pluto's surface: T 0 ≈ 40 K and p 0 ≈ 1 Pa. Other parameters for the N 2 atmosphere at the surface of Pluto are computed and presented in Table 9.1, and the probability of molecule absorption/emission is taken α ≈ 1 for Hertz-Knudsen law (Eq. 9.5). With these numbers, we evaluate the dimensionless numbers is related to u * , and we have = 3.3 10 -5 , for u * ≃ 1 m/s. Accordingly, the multiplicative factor f σ = 1.4 10 -8 s -1 in front of the growth rate (Eq. 9.70 and 9.79).

In the following, discussions will be divided into two parts according to the two dominant destabilising factors: the penitentes instability, associated with solar radiation and due to the self-illumination of the surface, and a diffusion-convection instability associated with the heat flux coming from the atmosphere towards the surface. We will study the influence of the parameters on dispersion relation in a rather abstract way, regardless the values of the dimensionless numbers on Pluto. Then we will apply the analysis to the patterns observed on Pluto.

Table 9.1 Parameters of N 2 -atmosphere at the surface of Pluto.

Item Value

Mass density ρ = mp 0 /(k B T 0 ) 8.4 10 -5 kg/m 3 Thermal velocity V th = 8k B T 0 /(πm) 1.7 10 2 m/s Mean free path ℓ m = k B T 0 /( 2πa 2 P 0 )

1.5 10 -3 m Kinematic viscosity ν = 1 3 V th ℓ 8.6 10 -2 m 2 /s Specific heat capacity C v = 5 2 k B /m 7.4 10 2 J/K/kg

Instability due to heat diffusion and convection

Heat flux from the atmosphere is characterized by + . It is due to the heat diffusion, and mostly the convection and turbulent mixing governed by hydrodynamics. We set Ω = 0 to study alone the instability characteristics due to heat flux.

As a start, we take the simplest situation by taking a fixed value for the rescaled wind velocity, = 3. ). The maximum growth rate increases when rising + , with a most unstable mode located around ≈ 10 -2 . To understand this instability, we computed the perturbed fields for the most unstable mode. The diffusion coefficient is modulated by the bottom (in first order approximation) and follows it, due to the fact that it is based on the distance to the ground (Fig. 9.5a). The velocity field is shown with stream functions in Fig. 9.5b. The streamlines get closer at the crests, which means a larger velocity there than at the troughs. The resulting temperature field is plotted in Fig. 9.5c.

One sees that the temperature iso-coutours are closer to each others in the trough. So the heat flux is higher at troughs than that at crests, which means more energy at the troughs for sublimation. The temperature field is resulted from the power balance at the surface, which combines temperature, heat flux and pressure, from the heat transport along the streamlines by convection, and from diffusion which spreads the field. So, the wavelength is selected by the balance between diffusion and convection, and scales as the viscous length ν/u * . Also, one sees in Fig. 9.4b, there is a threshold + for the instability. This is due to a stabilising pressure effect, which will be discussed later on, will dominate the system if + is too small. For the current = 3.3 10 -6 , the threshold flux is around + ≈ 10 -15 .

The pressure effect is characterized by , since p ∝ ρu 2 * ∝ 2 . The dispersion relation and instability diagram is therefore computed by varying , keeping + = 3.0 10 -12 fixed as well as = 0, -= 0 and Σ = 0. As shown in Fig. 9.6a, we observed the similar instability as in Fig. 9.4a, while the growth rate is decreasing with increasing . Thus, pressure is a stabilising effect. This comes from the fact that a higher pressure corresponds to a higher density and therefore a higher condensation rate, as described in Hertz-Knudsen law (Eq. 9.5). On crests, pressure is lower because of a higher velocity. On troughs, conversely, pressure is higher (Fig. 9.5b). This favors condensation in the troughs and sublimation on crests, which is stabilising. This is also why there is an upper bound in for instabilities (Fig. 9.6b), where the heat flux is not large enough the overcome the stabilising pressure. For + = 3.0 10 -12 in this calculation, the maximum is about 5 10 -5 .

The influence of the base state sublimation rate is characterized by the rescaled vertical velocity . We fix = 3. To study the influence of the thermal radiation of the solid surface, we compute the dispersion relation by varying Σ , with = 3.3 10 -6 , + = 3.0 10 -12 , -= 0 and = 0.

Σ enters the instability from temperature boundary condition on the surface (Eq. 9.75). Introducing it into the growth rate (Eq. 9.70), it gives:

σ = f σ -1 1 + + 4 Σ [( + 4 Σ )(S n (0) + 2i) -Y (0)] . (9.112)
One can see that Σ appears in front of the stabilising pressure term. When Σ is small, the growth rate is close but lower than that from Σ = 0 (red and black curves in Fig. 9.8).

As Σ is increased, the stabilising pressure becomes dominant and stabilises the flat state (violet, blue and green curves in Fig. 9.8).

Instability due to solar radiation

To study the behavior of the system due to the solar radiation and surface illumination, we set + = 0, and take a constant Ω = 0.25 (defined in Eq. 9.24). Similar to Subsection 9.3.1, we take Σ = 0, -= 0 and = 0 as a start. Dispersion relation and instability diagram are computed varying , as shown in Fig. 9.10. When is small, the self-illumination term is larger than the pressure effect and the flat base state is unstable. As shown by the orange curve in Fig. 9.10a, the growth rate increases linearly with -1 which points to the fact that self-illumination is a scale-free effect. The dispersion relation superimposes well on the analytical scaling (Eq. 9.79). As rises, the system shows selective instability over -1 , and the width of unstable zone and the growth rate are both decreasing with increasing due to the stablising pressure effect (red, blue and green curves in Fig. 9.10a).

As shown in Fig. 9.10b, there are three zones on the stability diagram in --1 plane. For small , the base state presents a scale free instability because the system is driven by the solar radiation. As increases, the instability becomes selective in wave number, which comes from the competition between the wind pressure effect and solar radiation. For even larger , the pressure effect becomes dominant and stabilises the flat state. It is easy to understand: the marginal stability curve goes up if we increase ψ , because it increases the scale-free unstable zone. Also, the maximum for the instability (the maximum values of dashed curves in Fig. 9.10b) rises, when pressure effect becomes dominant. . Regarding the surface radiation, as shown in Fig. 9.12, it behaves the same as in the same way as for the heat diffusion-convection case (Fig. 9.8). It has a stabilising effect.

The influence of the heat flux in the solid ice is shown in Fig. 9.13. It is destabilising as expected from the Mullins-Sekerka instability (as in Fig. 9.9).

Application to Pluto

Knowing the effects of all these parameters, we now study the patterns on Pluto. All the fluxes discussed above are coupled by the power balance, as in Eq. 9.78. With the collected parameters, we can estimate Fig. 9.12 Dispersion relation computed at different Σ . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. ψ = 10 -10 , = 7.6 10 -5 , -= 0 and = 0 for all the curves; and Σ = 0 for blue ones, Σ = 10 -10 for red ones, Σ = 10 -8 for green ones, Σ = 10 -6 for black ones. It should be noted that we use here the solar radiation 1.3 W/m 2 at the measurement spot of New Horizons for ψ . From the study on penitentes instability, we know that the heat flux in the solid is negligible [START_REF] Claudin | Physical processes causing the formation of penitentes[END_REF], so we take here -= 0. What remains is the fraction of the net flux ψ -Σ that goes into sublimation and that goes into heat flux from atmosphere + . We introduce the ratio

Σ = 3.
ε = + ψ -Σ , ( 9.114)
and vary it from 0 to 1 to see its effect on the instabilities. From the high resolution images taken by New Horizons, we have measured the wavelength of the patterns, and we take an average wavelength λ = 400 m for the following discussions.

In showing a local maximum around = 10 -2 . However, it is still more unstable at large wave number (orange curves in Fig. 9.15). One therefore would expect the patterns to develop at small wavelength, which is not consistent with the observed patterns length scale either.

When half of the heat flux arriving at the ground comes from the atmosphere (ε = 0.5),

the instability shows a nice selectivity over wave numbers, with a most unstable mode locating around = 10 -2 . Moreover, it is stable at large wave numbers (blue curves in Fig. 9.15). If one goes on increasing the fraction of heat flux from atmosphere, one gets the same instability properties (green and red curves in Fig. 9.15), except that the unstable plateau disappears for ε = 1 case because = 0. Taking the most unstable mode, we estimate the wind velocity with respect to the observed wave length on Pluto, and it gives u * = 0.13 m/s. This value is consistent with the one involved in in the present calculation, and also consistent with the prediction from the atmosphere model in Chapter 8. In Fig. 9.16, we show the influence of the wind velocity on the instability. We can clear see that larger wind induces larger pressure, stabilising the system.

km

Fig. 9.14 Rhythmic dune-like patterns observed on the south of the Sputnik plain on Pluto, showing that the pattern adapts to the relief and has a direction that is not selected by the orientation with respect to the sun.

Conclusions

In this part, we have developed an model for the rhythmic, dune-like patterns observed on Pluto. In modelling the dynamics of the atmosphere, we find that the thermal parameters and the velocity do not vary much at both rotation and revolution time scales compared to that of comet 67P. We developed a dimensionless number (Eq. 8.4) showing that most of the first order correction of the solar radiation contributes to the heat flux, instead of modulating the temperature, and this effect is even stronger for Pluto. Considering this quasi-steady uniform atmosphere, we have investigated the instability in the parameter space.

The discussions have shown the presence of four destabilising mechanisms: the diffusion/convection of heat, the blowing velocity from the surface, the Mullins-Sekerka instability associated with a thermal flux coming from the ice and the penitentes instability associated with self-illumination of the surface. In turn, from a Bernoulli-like effect, pressure turns out to be stabilising. In the end, the model is tested against the patterns on Pluto, and we have shown that heat flux from the atmosphere due to convection and turbulent mixing are likely to be responsible for the the emergency of these patterns. Therefore,

Aeolian sand ripples instability

As we have mentioned in Section 1.2, granular ripples form in various environments, resulted from the interaction between topography and sediment transport. The formation of subaqueous ripples is now rather well understood. Regarding aeolian ripples, ubiquitous in deserts or at the surface of sand beaches, it has long been recognized that they do not form via the same dynamical mechanism as dunes [START_REF] Wilson | Aeolian bedforms: Their development and origins[END_REF]. Their formation and evolution are not fully explained yet, despite the significant attention that has been given ever since last century. In this chapter, we revisit this question and present progresses we made on the subject.

Introduction

The modeling of aeolian ripple's formation dates back to the work of Bagnold [START_REF] Bagnold | The Physics of Blown sand and Desert Dunes[END_REF], who suggested that the ripple wavelength is equal to the characteristic length of the saltation jump. This idea was however proved problematic because experimental observations showed that grains jump over distances that are much larger than the ripple wavelength, and moreover this model cannot predict the merging of ripples and the resulting "coarsening" of ripple patterns observed in both in the field and in controlled experiments [START_REF] Sharp | Wind ripples[END_REF]. Anderson has proposed an alternative to Bagnold's model, pointing out that reptation plays a key role in aeolian ripple formation [START_REF] Anderson | Simulation of aeolian saltation[END_REF]. The hop length distribution has later been described in terms of a 'splash function' [START_REF] Ungar | Steady state saltation in air[END_REF]. This model was further improved by taking possible rolling of reptating particles down the slopes of sand surfaces into account, and the improved model was able to simulate realistic shapes of aeolian ripples and ripple coarsening [START_REF] Bouchaud | A model for the dynamics of sandpile surfaces[END_REF][START_REF] Terzidis | A model for ripple instabilities in granular media[END_REF][START_REF] Prigozhin | Nonlinear dynamics of aeolian sand ripples[END_REF]. From Anderson's model, Yizhaq et al [START_REF] Yizhaq | Blown by wind: nonlinear dynamics of aeolian sand ripples[END_REF] proposed a continuous model considering ripples as two-dimensional structures on the sand surface, and the variability of saltation flux due to surface undulations. This model, however, does not include sand segregation effect and is not able to capture the disappearance of ripples at high wind velocity. All the above models (and others like Manukyan and Prigozhin [START_REF] Manukyan | Formation of aeolian ripples and sand sorting[END_REF], Terzidis et al [START_REF] Terzidisa | A model for ripple instabilities in granular media[END_REF] and Csahók et al [START_REF] Csahók | Dynamics of aeolian sand ripples[END_REF]) relate the emergence and evolution of aeolian ripples to the reptation transport caused by the impacts on the ground of grains entrained by the wind into saltation in the transport layer. The characteristic height of the transport layer is indicated by Bagnold focal point. The majority of the particle are transported and all the reptation happens below this point. The transport substantially slows down the wind close to the surface.

Theoretical and numerical models [START_REF] Anderson | Simulation of aeolian saltation[END_REF][START_REF] Ungar | Steady state saltation in air[END_REF] report that below the focal point, the air velocity is almost invariant with the friction velocity u * . Therefore, the characteristics of the grain population that dominate transport (trajectories and velocities) is independent of u * , so is the length at which aeolian ripples emerge. However, field and laboratory measurements [START_REF] Sharp | Wind ripples[END_REF][START_REF] Seppälä | Wind tunnel studies of ripple formation[END_REF][START_REF] Andreotti | Aeolian sand ripples: experimental evidence of fully developed states[END_REF] have shown that the wavelength of ripples, which is of the order of ten centimetres on Earth, is much larger than the reptation hop length, of a few grain sizes, and both the wavelength and the propagation speed increases linearly with u * (Fig. 10.1). This strongly suggests that a fundamental ingredient is missing in the existing models of ripples instability mentioned above. Further studies are necessary to explain this discrepancy. Recently, this question was revisited numerically in the group in which this thesis was done [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. The results have recovered the linear increase of the wavelength and propagation speed with u * (Fig. 10.1), and have shown that 'resonant grain trajectories', tuned with the ripple wavelength, are closely related to the dynamical mechanisms controlling the formation of aeolian ripples, and that the product of the ripple wavelength and velocity is a proxy for the sediment flux. In collaboration with Orencio Durán and A. Brad Murray, we propose in the following an analytical model for the aeolian ripple instability, as a successor of the precursive work [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. The content is organised as follows. Firstly, a simplified transport model is described. Then, we explain a model for the bed evolution considering as a collective effect by a crater function. The models are compared to the numerical data, and we also show some preparatory experimental work at the end.

A simplified transport model

We propose here a simplified transport model, in which successive trajectories of a single grain are computed. We assume a time scale separation between sediment transport and bed evolution, i.e. we hypothesize that the bed profile is fixed during the time we compute the grain trajectories. For the base state, we consider a particular grain trajectory over a flat bed for a steady and homogeneous case, starting at position x ↑ , of length ℓ 0 and thus arriving at position x ↓ = x ↑ + ℓ 0 , as in Fig. 10.2a. The grains are assumed to start their trajectory with a random velocity, whose distribution is equivalently described by the hop-length distribution P(ℓ 0 ). Here we idealize the prediction from the numerical simulation [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF], and describe it as

P(ℓ 0 ) = 1 ln (ℓ s /ℓ d ) 1 ℓ 0 , ( 10.1)
where we have explicitly introduced the lower and upper cutoffs, ℓ d and ℓ s respectively. ℓ s is transport length scale, proportional to q sat /ϕ ↑ (see Eq. 10.9 below), i.e. increasing linearly with the wind shear velocity. ℓ d is on the order of the grain size d. We thus have ℓ s ≫ ℓ d , and we can consider the normalisation ln (ℓ s /ℓ d ) as fairly constant, even for varying wind strengths.

Hop length modulation

In the case of a modulated bed Z(x) = Ẑ e ikx , trajectories are modified, with λ = 2π/k the wavelength of and Ẑ the amplitude of the modulated bed. At the first order in the bed elevation profile, we can write their disturbed length ℓ in a generic way as:

ℓ -ℓ 0 = Ẑ e ikx ↓ , (10.2)
where is a complex number -note: as usually understood with complex notations, real part must be taken. This hop length modulation comes from the geometry of the rippled bed, and is also due to the modulation of the wind induced by the bed. We successively analyse below these two independent contributions and compute the corresponding modulus and the phase of as functions of the primary parameter kℓ, and for different secondary parameters.

Geometrical effect

We note α ↓ the arrival angle of the grain on the bed (counted positive) and consider it as a constant, independent of ℓ 0 . On a flat bed (Fig. 10.2a), the equation of the trajectory close to its end is z = tan α ↓ (x ↑ + ℓ 0 -x). Starting on a modulated bed (Fig. 10.2b,c), at altitude Z(x ↑ ), the trajectory of a grain starting with the same initial velocity crosses the bed again at x ↓ , where ℓ is the modified hop length. Neglecting wind modulation, the trajectory remains identical as that on a flat bed and one then gets, at the linear order, the

condition Z(x ↓ ) = Z(x ↑ ) + tan α ↓ (x ↑ + ℓ 0 -x ↓ ), which solves into ℓ -ℓ 0 = 1 tan α ↓ Z(x ↑ ) -Z(x ↓ ) (10.3)
Introducing Z(x ↑ ) = Ẑ e ik x ↑ and Z(x ↓ ) = Ẑ e ikx ↓ into Eq. 10.3 and comparing with Eq. 10.2, one obtains an explicit expression for the geometrical contribution to :

g = 1 tan α ↓ e -ikℓ -1 . ( 10.4)
where we have used the definition x ↓ ≡ x ↑ + ℓ.

Wind effect

We now consider the influence of the modulation of the wind velocity. The corresponding contribution w is computed numerically, from the integration of the particle trajectories [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. The wind velocity ⃗ u enters the equation of the grain motion through the drag force.

For the sake of simplicity, we choose a wind velocity field of the form:

u x (x, z) = u * κ µ(z -Z) (10.5) u z (x, z) = Z ′ (x) u * κ µ(z -Z) (10 .6)
In the molecular dynamics simulations, the profile µ(z) is determined from the case of steady and homogeneous transport. It takes into account the negative feedback of transport on the wind. For comparison, we have also considered the analytical formula µ(z) = ln(1 + z/ξ 0 ), with ξ 0 the hydrodynamical roughness. For further test of robustness of the results, we have also used two other wind models where the effect of transport is ignored: the perfect flow, and a Reynolds averaged description of the flow on a relief, in the spirit of [START_REF] Jackson | Turbulent wind flow over a low hill[END_REF] [START_REF] Jackson | Turbulent wind flow over a low hill[END_REF], as done for dunes [START_REF] Fourrière | Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by non-linear pattern coarsening[END_REF][START_REF] Kroy | Minimal model for aeolian sand dunes[END_REF].

The modulus and the phase of the function g and w are displayed in Fig. 10.3 as a function of ℓ/λ, for the profile µ(z) computed from the simulations [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. It turns out that, to a first approximation, the two contributions g and w are similar in magnitude and behave roughly in the same way.

Flux modulation

We introduce the distribution ψ(ℓ, x ↓ ), such that ψ(ℓ, x ↓ )dℓd x ↓ is the volume of the grains

ψ(ℓ, x ↓ )dℓd x ↓ = P(ℓ)dℓ ϕ ↑ (x ↑ )d x ↑ , (10.7)
Both ψ and ϕ ↑ have a zeroth order contribution corresponding to the flat bed, and a first order correction for the modulated bed. We denote the decomposition by ψ = ψ 0 + ψ 1 and ϕ ↑ = ϕ 0 ↑ + ϕ 1 ↑ . Considering the homogeneous and steady base state, we have ψ 0 (ℓ, x ↓ ) = ψ 0 (ℓ) and ϕ 0 ↑ (x ↑ ) = ϕ 0 ↑ , and Eq. 10.7 then simply gives

ψ 0 (ℓ) = P(ℓ)ϕ 0 ↑ . (10.8)
The saturated flux is given by q sat = ℓ s ℓ d ℓ ψ 0 (ℓ) dℓ. Using the expression for the hop length distribution (10.1) and Eq. 10.8, it gives

q sat = ϕ 0 ↑ ℓ s -ℓ d ln (ℓ s /ℓ d ) ≃ ϕ 0 ↑ ℓ s ln (ℓ s /ℓ d ) . ( 10.9)
where we have used ℓ s ≫ ℓ d .

For the first order contribution ψ 1 , several contributions must be added. The first term is due to the disturbance ϕ 1 ↑ evaluated at the starting point x ↑ = x ↓ -ℓ. The second term comes from the change of variables (ℓ 0 , x ↑ ) ↔ (ℓ, x ↓ ) by Eq. 10.2 and x ↓ ≡ x ↑ + ℓ, whose Jacobian reads:

∂ ℓ 0 ∂ ℓ ∂ x ↑ ∂ ℓ ∂ ℓ 0 ∂ x ↓ ∂ x ↑ ∂ x ↓ = 1 -k Ẑ e ik x ↓ ′ -1 -ik Ẑ e ik x ↓ 1 = 1 -k Ẑ e ikx ↓ ′ + i , ( 10.10)
where ′ = ∂ /∂ (kℓ). If we consider only the geometrical contribution by plugging Eq. 10.4 into 10.2 , the Jacobi reduces to 1 + ik tan α ↓ Ẑ e ikx ↓ , which corresponds to the slope effect and is to be discussed at the end of this section. Last, there is a term involving the derivative of P(ℓ) when accounting for the disturbance of the hop length. Therefore, the expression of ψ 1 finally reads: (10.11) From the decompostion of ψ and ϕ, one could take the generic form for the fisrt order contributions:

ψ 1 (ℓ, x ↓ ) = P(ℓ)ϕ 1 ↑ (x ↓ -ℓ) -P(ℓ)ϕ 0 ↑ ′ + i k Ẑ e ikx ↓ -ϕ 0 ↑ P ′ (ℓ) k k Ẑ e ikx ↓ .
ϕ 1 ↑ (x ↑ ) = ϕ 0 ↑ F ↑ k Ẑ e ikx ↑ (10.12) ψ 1 (ℓ, x ↓ ) = ϕ 0 ↑ P(ℓ) k Ẑ e ikx ↓ ( 10
= F ↑ e -ikℓ - ′ + i - P ′ (ℓ) kP(ℓ) , ( 10.14)
The problem closes under the assumption that sediment transport is in a steady state, i.e.

when deposition and erosion rates exactly balance:

ϕ ↓ (x) ≡ ψ(ℓ, x)dℓ = ϕ ↑ (x) , (10.15)
with Eq. 10.11 and 10.12 plugged into 10.15, it gives

F ↑ = - P(ℓ) ( ′ + i ) dℓ + P ′ (ℓ) k dℓ 1 -P(ℓ)e -ikℓ dℓ . ( 10.16)
The analytical expression for the factor corresponding to the geometrical contribution (10.4) can be derived. Using the normalisation condition P(ℓ)dℓ = 1 and the fact that the integrals P(ℓ)e -ikℓ dℓ as well as 1 k P ′ (ℓ) dℓ are dominated by the behaviour of their integrand at small ℓ when ℓ d → 0 because of P(ℓ) ∼ ℓ -1 . Eq. 10.16 thus reduces to .17) ives: Plugging Eq 10.17 and 10.4 into 10.14, one obtains

F ↑ ∼ i/ tan α ↓ . ( 10
tan α ↓ ∼ i(1 + e -ikℓ ) + e -ikℓ -1 kℓ , ( 10.18)
where we have used P ′ (ℓ)/P(ℓ) ∼ ℓ -1 . As can be seen from the linear structure of Eqs. (10.14) and (10.16) in , the contribution due to the wind modulation adds up.

For comparison, it is interesting to compute the flux modulation rate corresponding to Anderson's picture [START_REF] Anderson | A theoretical model for aeolian impact ripples[END_REF][START_REF] Anderson | Eolian ripples as examples of self-organization in geomorphological systems[END_REF]. In Anderson's model, the geometrical effect taken into account is the change of the cross section due to the bed slope. In our formalism, this is encoded in the Jacobian, which relates the lengths and coordinates. With the sole contribution of Eq. (10.10), and for = g , one simply gets = i/ tan α ↓ . In other words, there is a prediction for a modulation of the impacting flux, but this modulation is independent of kℓ: There is no resonance when ℓ is a multiple of λ. [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF], with u * /u t = 4, and with increasing wave numbers k labelled from red to violet. Analytical prediction (Eq. 10.18) corresponding to the geometrical contribution, with an adjustable overall prefactor on the order of 2.

A simplified model for bed evolution

The simplified model described in Section 10.2 is able to reproduce the shape of the flux modulation rate, but, as far as ripples are concerned, shows two shortcomings:

1. It models trajectories starting from the bed and which therefore cross the interfacial boundary layer, but ignores mid-air collisions occurring in this layer; 2. It ignores erosion/deposition processes and therefore, by construction, does not lead to any evolution of the bed.

In order to fix these problems, we present in this section an approach based on a crater function. This function ∆(xx ↓ , ℓ) quantifies the change of the bed profile at position x, after a grain making a hop of length ℓ has collided with the static grains at position x ↓ . ψ(ℓ, x ↓ ) expresses the probability that a grain arrives at x ↓ after a hop of length ℓ.

Following these definitions, the equation governing the bed evolution reads: .19) This approach assumes a clear separation of moving grains that rebound after they hit the static bed, and those which do not. The latter may be called reptons, and are effectively encoded in the crater function. The former are the saltons, but it should be emphasised again that they are not distinguished from the reptons by a criterion on the hop length or the hop height. The hop length distribution of the saltons is an input of the description -and is not computed self-consistently by means of a given splash function. We further assume that it is well approximated by P(ℓ), which makes no distinction between saltons and reptons.

∂ t Z = ψ(ℓ, x ↓ )∆(x -x ↓ , ℓ)d x ↓ dℓ. ( 10
The reason is that the computation of the modulation rate of the distribution ψ with P (see Section 10.2), although leading to a stationary bed, gives a quantitative fit of the curve (kℓ) measured in the full simulations (Fig. 10.4), where the ripples develop. Another important point is that the processes at work in the interfacial layer can be effectively encoded in P, taking it as a power law, as measured in the simulations. Finally, we assume for ∆ a self-similar form .20) This form is consistent with the presence of scale-free processes in the interfacial layer.

∆(x -x ↓ , ℓ) = 1 ℓ x -x ↓ ℓ . ( 10
Importantly, there is no explicit dependence on the grain diameter d. is a dimensionless function -following Eq. 10.19, the function ∆ has the dimension of the inverse of a length -which we now wish to determine in order to reproduce the properties of the ripples, their growth rate and propagation velocity in particular. Expression (10.20) is for a flat bed. In the case of a modulated bed as considered below, following the zeroth order (10.20), the first order correction to the crater function can be written as

∆ 1 (x -x ↓ , ℓ) = - 1 ℓ x -x ↓ ℓ ik Ẑ e ikx ↓ .
(10.21)

Note that we consider the modulation of the bed at the point x ↓ of arrival. The minus sign in this expression comes from the fact that the crater is increased when the bed slope at x ↓ is negative.

Taking Eq. 10.19 for a flat steady bed, and recalling that ψ 0 = P(ℓ)ϕ 0 ↑ in this case, we obtain 0 = ϕ 0 ↑ P(ℓ) This expresses the mass conservation.

For an undulated bed, two contributions must be considered. First we take ∆ at zeroth order, while ψ is taken at first order in Eq. 10.19. This contribution is associated with a complex growth rate Ω d -later identified as the destabilising term of the total growth rate Ω. σ = ℜ(Ω) is the actual growth rate of the bed undulations. c = -ℑ(Ω)/k is the propagation velocity of the pattern. Following the definition of ψ 1 (Eq. 10.13) with the modulation rate , we obtain: with ˆ (r) the Fourier transform of (s).

The second (later identified as the stabilising term) contribution to the growth rate comes from ∆ taken at first order (Eq. 10.21) and ψ taken at zeroth order (Eq. 10. where we have changed the variables and have used the definition of Fourier transform, similarly as in Eq. 10.25 and 10.26, and ˆ (r) the Fourier transform of (s).

For the sake of a simple analytical form, we take for the crater function the derivative of the Gaussian (s) = (ss)e -[a(s-s)] 2 .

(10.28)

The coefficient a is a positive real number, which encodes for the spatial extension of the crater (width ∝ 1/a). s is the shift between the centre of the crater function and the point of impact. This functional form satisfies the condition (10.23). One can then analytically compute ˆ , which reads:

ˆ (r) = ∞ -∞
e -i rs-[a(s-s)] 2 (ss)ds = -i π 2a 3 r e -irs-r 2 /(4a 2 ) . (10.29)

The remaining integral over r can then be computed for the function given by Eq. 10.18, in the limit ℓ d → 0 and ℓ s → ∞. The integration can be done analytically and involves special functions. In the limit of small s, it gives at first order:

ln (ℓ s /ℓ d ) tan α ↓ k ϕ 0 ↑ Ω d = -i π 2a 3
∞ 0 e -irs-r 2 /(4a where erf is the error function and 2 F 2 is the generalized hypergeometric function. Note that the analytical expression for an arbitrary s exists, but is not compact enough to be easily displayed. For a ≃ 1, this expression gives a positive real part and a negative imaginary part, which are both of order one (Fig. 10.5). Furthermore, we have checked that, for ℓ d → 0 and small s, this asymptotic behaviour independent of k is reached as soon as kℓ s ≳ 2a. where we have used ˆ (r) = -i π 2b 3 r e -irs-r 2 /(4b 2 ) , and Eq. 10.9 for q sat . In the limit of ℓ d → 0 and s → 0, this expression reduces to With the analytical model, we would like to recover the same structure for the dispersion relation as the one reported in the numerical work [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF], in which the growth rate is derived as a parabolic function σ = c 1 ϕ 0 ↑ kc 2 q sat k 2 , (10.34) where c 1 and c 2 are multiplicative constants of order 1. It contains a destabilising term, linear in k, and a stabilising term, quadratic in k (Fig. 10.7). In the model, we also consider the growth rate in two parts, a destabilising term Ω d and a stabilising term Ω s . Ω d is expressed in Eq. 10.25. It is the product of k, and an integral that is independent of k. In the limit of ℓ d → 0 and ℓ s → ∞, the integral is finite, as calculated in Eq. 10.30 and shown in Fig. 10.5. It is indicated that Ω d is positive and proportional to k. As for Ω s , it is given in Eq. 10.27, and it is found negative and proportional to k 2 in the limit of small ℓ d and small s (Eq. 10.32 and 10.33). So, we successfully recover the two contributions in the dispersion relation.

Ω s q sat k 2 = -

Discussions

In this chapter, we consider the aeolian ripple instability in an analytical way. Sediment transport is modelled computing the successive trajectories of the single grain and the bed evolution is characterised by a crater function. We have compared the model with numerical data that well recover the experimental facts, and good agreement has been seen in the modulation of saltation flux rate and in the dispersion relation.

Experimentally, one has observed the modulation of saltation flux. As in Fig. 10.8, one can clearly see the saltation rain above the rippled bed is modulated and the grains preferentially hit the bed upwind of the ripple crests. To make a systematic and quantitatively test, we built a wind tunnel. As sketched in Fig. 10.9, it is composed of a fully transparent Plexiglas rectangular tank of length L = 1.5 m, width W = 0.2 m, and depth h = 50 mm, fitted to the bottom of a horizontal channel of rectangular cross-section, with a honeycomb over a sinusoidal rippled bed. For this, we will prepare the sinusoidal rippled sand beds with different the wavelength and amplitude, and a fast camera will be employed to record the motion of the sand grains from the side view to study the hop length distribution, and its modulation.

Chapter 11

Conclusions and perspectives

In this thesis, we have focused on patterns observed on a flexible or an erodible substrate induced by hydrodynamical instabilities. Four different problems have been studied whose major conclusions and perspectives are summarized here:

1. Travelling waves on a flexible sheet induced by an axial wind are modelled experimentally and theoretically. We conclude that the frequency is selected by the balance between the system inertia and dynamic pressure, while the wavelength results from the balance between sheet elasticity and dynamic pressure; and the amplitude distribution comes mainly from the geometrical constraints. In the theoretical frame, we have assumed an unbounded homogeneous sheet, which thus constrain the predicted scalings to the limit kL 0 ≫ 1. In addition, we have hypothesized an inviscid description for the fluid flow, and it is thus impossible to study the viscous effect on the system, such as shear stress over the sheet surface, which could be of importance in some cases. As a consequence, future work will have to seek for modes which satisfy the boundary conditions at both end, breaking the invariance along the wind direction, and to include the viscous effect in the model. Regarding the application, this work can be of interest for energy harvesting mechanism by using a piezoelectric flag waving in a uniform axial flow [START_REF] Xia | Fluid-solid-electric lock-in of energy-harvesting piezoelectric flags[END_REF][START_REF] Pineirua | Influence and optimization of the electrodes position in a piezoelectric energy harvesting flag[END_REF], as in Fig. 2.11. 3. An explanation is proposed for the rhythmic, dune-like patterns observed on Pluto.

These patterns are hypothesized to be sublimation dunes created by differential condensation/sublimation of nitrogen ice given their location on Pluto. Firstly, we study thermo-hydrodynamics of the atmosphere, and the results show gentle variations in the parameters for both rotation and revolution scales. Considering the characteristics of the atmosphere and the linear patterns in the nature, we model the sublimation patterns in a general way, by linking their emergence and evolution to the thermal and hydrodynamic instability of an interface between a sublimating ice bed and its turbulent vapor flow. The discussion has focused on two major destabilising factors:

the heat flux from the atmosphere due to hydrodynamics and the solar radiation. In the end, the model has been tested against the patterns on Pluto. Given the observed wave length, we conclude that heat flux from the atmosphere due to convection and turbulent mixing is responsible for the the emergency of these patterns. Admittedly, more detailed atmosphere model is necessary to support this conclusion in the future work.

4. Aeolian ripple instability is revisited in an analytical way. In the model, successive trajectories of a single grain are computed to model the sediment transport, and a crater function is employed to describe the collective effect of the bed evolution.

Good agreement has been seen when comparing the theoretical prediction with the numerical data [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. The growth rate is considered containing a destabilising contribution and a stabilising contribution, and the results have shown that the analytical model recovers well the same structure for dispersion relation as reported in [START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF]. For the further work, it is necessary to perform experimental tests for the model, and it will be also interesting to study the nonlinear characteristics of the model.

List of symbols

A Amplitude of paper waves. S zz Dimensionless function related to the first order correction of τ zz .

T i Temperature at the sublimation interface.

T 0 Base state temperature.

U r Radial velocity components in outer layer.

U θ Azimuthal velocity components in outer layer.

U Dimensionless function related to the first order correction of x velocity.

V fall Settling velocity.

V th Thermal velocity.

V Wind velocity. Dimensionless function related to the first order correction of psi.

Phase componet of the basal shear stress with the bottom.

Quadrature componet of the basal shear stress with the bottom.

Gravitational constant.

Dimensionless number relating transitional Reynolds number to pressure lag.

Latent heat of solid ice sublimation.

Dimensionless function related to the first order correction of pressure field.

the susceptibility with respect to the shear stress.

Reynolds number with respect to wavenumber.

the susceptibility with respect to the slope.

Dimensionless function related to base state temperature field.

Dimensionless function related to the first order correction of x direction velocity.

Dimensionless function related to the first order correction of z direction velocity.

µ Effective friction coefficient. w Specific enthalpy.

General introduction 1 . 1 of figures 1 . 1 xivof figures 2 . 8 33 2. 9 35 3. 1 43 4. 1 54 4. 2 57 4. 3 57 4. 4 58 4. 5 59 4. 6 60 5. 1 62 5. 2 64 5. 3 65 5. 4 67 5. 5 68 5. 6 of figures xv 6 . 1 104 9. 3 120 9. 7 120 9. 8 121 9. 9 123 xviof figures 9 .of tables 2 . 1

 111128339351431542573574585596601622643654675686611043120712081219123921 Flag flapping instability induced by wind . . . . . . . . . . . . . . . . . . . . . 1.2 Granular patterns on an erodible bed . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Sublimation patterns on an ice bed . . . . . . . . . . . . . . . . . . . . . . . . . List Flapping national flags and Huineng's story . . . . . . . . . . . . . . . . . . . . 1.2 Sketch of Bernoulli's principle and flow visualization of the wake structure downstream of a flexible filament showing the vortex shedding . . . . . . . . 1.3 Granular patterns under various environments . . . . . . . . . . . . . . . . . . 1.4 Schematics featuring the modes of sediment transport in the aeolian and subaqueous cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Experimental measurements of the saturated flux in the case of bed load transport under water and of L sat in the case of suspension . . . . . . . . . . 1.6 Schematic of the ripple instability mechanism . . . . . . . . . . . . . . . . . . 1.7 Patterns on the ice surfaces due to sublimation . . . . . . . . . . . . . . . . . . 2.1 Sketch of the experimental setup for the paper waves . . . . . . . . . . . . . . 2.2 Typical experimental measurements and data processing . . . . . . . . . . . . 2.3 Measured angular frequency (a) and wavenumber (b) a fucntion of wind velocity V for paper sheet with ∆L = 4 cm. . . . . . . . . . . . . . . . . . . . . 2.4 Dimensionless growth rate σ (positive branches only) and the dimensionless angular frequency ω as a function as the dimensionless wavenumber, with the bending rigidity D = 10 -3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Dimensionless growth rate σ (positive branches only) as a function as the dimensionless wavenumber, for k < kc , Eq. 2.26; and dimensionless angular frequency ω as a function as of k, Eq. 2.25 . . . . . . . . . . . . . . . . . . . . 2.6 Cut-off wavenumber kc (a), maximum wavenumber km (b), maximum angular frequency ωm (c) and maximum temporal growth rate σm (d) as functions of the bending rigidity D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Experimental values of the angular frequency and wavenumber as functions of the rescaled wind velocity for the different sheet material only but with the same sheet lengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List Experimental values of the angular frequency and wavenumber as functions of the rescaled wind velocity for the paper material only but with different sheet lengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schematic for the geometrical relation between the arclength and the length between two lamps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.10 Rescaled longitudinal profiles of the wave amplitude . . . . . . . . . . . . . . 35 2.11 An energy harvesting setup by a piezoelectric flag in a uniform axial flow . Image of 67P by Rosetta's OSIRIS narrow-angle camera, orbit of 67P and the relative paths of the Rosetta probe and it, Philae touch-down site . . . . 40 3.2 View of the comet's neck (Hapi) region by OSIRIS narrow-angle camera bwfore and after perihelion, and Photograph of ripplemarks on the two lobes 42 3.3 Picture of the comet and its close coma . . . . . . . . . . . . . . . . . . . . . . . Schematics of the porous granular layer at the comet's surface . . . . . . . . Time evolution of the vapour temperature (left axis) and corresponding thermal velocity (right axis) just above the comet's surface . . . . . . . . . . Time evolution of the vapour density (left axis) and corresponding mean free path (right axis) just above the comet's surface . . . . . . . . . . . . . . . Time evolution of the vapour pressure just above the comet's surface, and a schematic for the surface flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . Outgassing flux qm as a function of the comet's heliocentric distance η . . . Schematic of the vapour flow at the comet surface driven by strong pressure gradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Estimation of the grain size on the surface of 67P . . . . . . . . . . . . . . . . The origin of the sediment transport threshold at the scale of a grain . . . . Dependence of the threshold shear velocity u t on the grain diameter d under water on the earth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time evolution of the velocity ratio u * /u t , calculated along the comet's orbit around the sun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dependence of the threshold shear velocity u t with the grain diameter d at perihelion, for afternoon conditions . . . . . . . . . . . . . . . . . . . . . . . . . Schematics featuring the modes of sediment transport in the cometary cases 70 5.7 Saturation length L sat in units of d as a function of the flow velocity at a grain size above the surface u(d) rescaled by the grain settling velocity V fall 72 List Basal shear stress components in phase (a) and in quadrature (b) with respect to the bed elevation, as functions of the rescaled wave number kν/u * 6.2 Relation betwwen λ/L sat and the rescaled saturated length L sat u * /ν predicted from the most unstable mode . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Relation between the wavelength and the mean grain diameter predicted at perihelion, for afternoon conditions . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Dispersion relation: dimensionless growth rate (a) and dimensionless propagation speed (b) as functions of the rescaled wavenumber kν/u * , computed at perihelion for d = 4 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 A global view of Pluto and New Horizons spacecraft, along with their trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Phase diagrams of four geologically relevant substances: nitrogen, ethane, carbon dioxide and water. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Dune-like rhythmic patterns on the surface of Pluto revealed by NASA's New Horizons spacecraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Departure shot of Pluto by New Horizons, showing Pluto's atmosphere backlit by the Sun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Pressure and temperature in Pluto's lower atmosphere . . . . . . . . . . . . . 8.3 Time evolution of the surface temperature of Pluto . . . . . . . . . . . . . . . 8.4 Time evolution of the atmosphere pressure near Pluto's surface . . . . . . . . 8.5 Time evolution of the atmosphere density near Pluto's surface . . . . . . . . 8.6 Time evolution of the friction velocity u * at Pluto's surface . . . . . . . . . . . 9.1 Schematic of a unstable sublimation surface . . . . . . . . . . . . . . . . . . . . 9.2 A given point x receives light from a portion of the surface . . . . . . . . . . . Bounds of the integral and integral I(η) giving the illumination profile for a sinusoidal surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 9.4 (a) Dispersion relations for different + . (b) Stability diagram in + --1 plane118 9.5 Perturbed fields showing the instability due to heat diffusion and convection 119 9.6 (a) Dispersion relations for different . (b) Stability diagram in --1 plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dispersion relations for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dispersion relation computed at different Σ . . . . . . . . . . . . . . . . . . . Dispersion relation computed at different -. . . . . . . . . . . . . . . . . . . 122 9.10 (a) Dispersion relations for different . (b) Stability diagram in --1 plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List 11 Dispersion relations for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.12 Dispersion relation computed at different Σ . . . . . . . . . . . . . . . . . . . 9.13 Dispersion relation computed at different -. . . . . . . . . . . . . . . . . . . 9.14 Rhythmic dune-like patterns observed on the south of the Sputnik plain on Pluto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.15 Dispersion relation computed at different flux distributions for Pluto . . . . 9.16 Dispersion relation computed at different for Pluto . . . . . . . . . . . . . 10.1 Selection of the ripple wavelength and propagation speed . . . . . . . . . . . 10.2 A particular grain trajectory over a flat (a) and a rippled (b and c) bed . . . 10.3 Modulus (a) and phase (b) of the factor involved in Eq. 10.2 as functions of ℓ/λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.4 Modulus (a) and phase (b) of flux modulation rate . . . . . . . . . . . . . 10.5 Complex growth rate Ω d as a function of the rescaled wave number k ℓ . . . 10.6 Growth rate Ω s as a function of the rescaled wave number k ℓ/b . . . . . . . 10.7 Dispersion relation measured from numerical simulations . . . . . . . . . . . 10.8 Experimental evidence of a modulated saltation flux . . . . . . . . . . . . . . 10.9 Sketch of the experimental setup for aeolian sand ripples . . . . . . . . . . . 10.10Aeolian sand ripples emerging from a flat bed in the wind tunnel . . . . . . A.1 A photo of Testo 405-V1 Metric Thermal-Anemometer . . . . . . . . . . . . . A.2 A photo of fast camera Phantom Miro M340 . . . . . . . . . . . . . . . . . . . List Values of the mass per unit surface m and the bending rigidity D of the three different sheets used in the experiments. . . . . . . . . . . . . . . . . . . . . . . 3.1 Parameters of 67P and its current orbit. . . . . . . . . . . . . . . . . . . . . . . 6.1 Ripple crest-to-crest distance measured on different regions of 67P . . . . . 7.1 Parameters of Pluto and its orbit. . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Parameters of nitrogen (N 2 ) and nitrogen ice (N s 2 ) . . . . . . . . . . . . . . . . 9.1 Parameters of N 2 -atmosphere at the surface of Pluto. . . . . . . . . . . . . . . A.1 Technical data of Testo 405-V1 Metric Thermal-Anemometer . . . . . . . . . A.2 Technical data of fast camera Phantom Miro M340 . . . . . . . . . . . . . . .

Figure 1 .

 1 1b describes a well-known and possibly apocryphal story about the Chinese Buddhist master Hui-Neng (AD 638-713), arguing with two other monks when they watched a temple flag flapping in the breeze. One claimed "The flag flaps!" and the other replying "No, it

Fig. 1 .

 1 Fig. 1.1 (a) National flags flapping; (b) An artistic rendering of the story of Hui-Neng and two monks in Tang dynasty. (Photo sources: web)

1. 1 Flag flapping instability induced by wind 3 Fig. 1 .

 131 Fig. 1.2 (a) Sketch of Bernoulli's principle, with the bold back line and the blue lines respectievly, showing the sheet and the streamlines, and indicating the pressure difference at troughs and crests with the black arrows. (b) Flow visualization of the wake structure downstream of a flexible filament made of silk thread and immersed in a running soap film [10].

Fig. 1 . 3

 13 Fig. 1.3 Granular patterns under various environments. (a) An aeolian dune and the ripples on its surface; (b) and (c) show a snow dune and snow ripples; (d) subaqueous sand ripples in the stream in the Altai Republic, Russia; (e) ripples observed on a Martian dune. (All photos' source: web.)

Fig. 1 . 4

 14 Fig. 1.4 Schematics featuring the modes of sediment transport in the aeolian (a) and subaqueous cases (b).

Fig. 1 .

 1 Fig. 1.5 (a)Measurements of the saturated flux in the case of bed load transport under water. (•) Luque & van Beek (1976)[START_REF] Luque | Erosion and transport of bedload sediment[END_REF], (•) data collected by Julien (1998)[START_REF] Julien | Erosion and sedimentation[END_REF]. Red solid line: q sat ∝ (Θ -Θ t ) 3/2 . (b) Experimental measurements of L sat in the case of suspension, with the distance x rescaled by the deposition length L dep = U d/V fall , where U the flow velocity and V fall the grain falling velocity. Markers and lines: green circles (van Rijn 1986[START_REF] Van Rijn | Sedimentation of dredged channels by currents and waves[END_REF]) and gray circles(Ashida & Okabe 1982 [41]) for net erosion, blue squares(Ashida & Okabe 1982 [41]) for net deposition, solid and dotted lines in red for exponential fits.

Fig. 1 . 6

 16 Fig. 1.6 Schematic of the ripple instability mechanism. (a)Rescaled basal shear stress τ/τ 0 on a smooth sinusoidal bottom (black line). Measurements for 2ζ 0 /λ = 0.0125 (red squares) and 0.05 (blue squares) [58]; solid lines: best fit with three harmonics. τ 0 is the basal shear stress for the reference state: flat bed. ξ 0 and λ are the amplitude and wave length of the bedform. (b) The wind velocity close to the surface (red arrow) is modulated by the topography. The maximum of wind, the crest of the bedform and the maximum the sediment flux are labeled by red dotted line, black dotted line and orange dotted line.

Fig. 1 .

 1 Fig. 1.7 (a) Field of Penitentes on the Upper Rio Blanco, Central Andes of Argentina. The blades are between 1.5 and 2m in height, image from enWiki. (b) Ice waves, image from Internet. (c) Ice waves on the ice surface observed in the cryosphere in Antarctica, image from website of National Snow and Ice Data Center of USA. (d) Spiral-shaped topographic ice waves at the surface of the North Polar cap of Mars, from the sublimation and condensation of CO 2 , with superimposed katabatic wind streamlines from Massé et al [65]. (e) A satellite image of the ice waves (or termed as ice megadunes) with 25 km resolution, a top right inset showing the map of Antarctica and location of the ice waves, credit: NASA Earth Observatory.

2 .

 2 The work in Part II is inspired by the unexpected bedforms on Comet 67P/Churyumov -Gerasimenko revealed by the Rosetta Mission of European Space Agency. In Chapter 4, we model the thermo-hydrodynamics of 67P's atmosphere. Based on the specific atmospheric properties, the sediment transport are analyzed in Chapter 5, and finally we discuss the nature of the bedforms and compare the predicted time and length scales with the observations in Chapter 6.

5 .

 5 Finally, the general conclusions and perspectives are given in Chapter 11.

Fig. 2 . 1

 21 Fig. 2.1 Sketch of the experimental setup. V is the air flow velocity, measured at a fixed position at the outlet of the tunnel. The distance between the two masts is L 0 = 1.5 m. The square section of the tunnel is 0.15 m × 0.15 m. The sheet is 4 cm wide and its total length is denoted as L = L 0 + ∆L. ∆L has been varied from 4 to 12 cm.

Fig. 2 . 2

 22 Fig. 2.2 Typical experimental measurements and data processing. (a) Single frame image of the experimental record, with the red lines showing the detected borders of the sheet. (b) Temporal variation of the deflection ζ of the sheet at a given location (here, x = 0.235 m). Solid line: Fourier fit of the form ζ(x, t) = Acos(ωtφ). (c) Spatial variation of the angular frequency ω along x. Solid line: fit by a constant. (d) Spatial variation of the phase φ along x. Solid line: linear fit of the form φ = kx. (e) Spatial variation of amplitude A along x. Solid line: linear fit to estimate the slope ≃ 10 -2 . All these data correspond to the paper sheet, at V = 8.6 m/s and with ∆L = 4 cm. Data for the plastic sheets or for other control parameters look very similar. In panels (b-d), 8 time series corresponding to 2000 frames at 1000 Hz (2 s in total) have been averaged.

Fig. 2 . 3

 23 Fig. 2.3 Measured angular frequency (a) and wavenumber (b) a function of wind velocity V for paper sheet with ∆L = 4 cm. Solid lines are the phenomenological fittings: (a) ω ∝ V and (b) k ∝ V 2/3 .

. 26 )

 26 All wavenumbers between 0 and k c are unstable. In between, σ shows a maximum, corresponding to the most unstable wavenumber km . It is the solution of d σ/d k = 0, which gives: 2 Dk 4 m + 9 Dk 3 m + 10 Dk 2 m -km -6 = 0. (2.27)

Fig. 2 . 4

 24 Fig. 2.4 Dimensionless growth rate σ (positive branches only) and the dimensionless angular frequency ω as a function as the dimensionless wavenumber, with the bending rigidity D = 10 -3 . Black and red curves indicate the behaviors of σ and ω for k < kc and k > kc , respectively. Circles show the positions of cutoff wave number kc , maximum wave number km and maximum growth rate σm .

Fig. 2 .

 2 Fig. 2.5 (a) Dimensionless growth rate σ (positive branches only) as a function as the dimensionless wavenumber, for k < kc , Eq. 2.26. The different lines correspond to different values of the bending rigidity D. Dotted-dashed line: D = 10 3 , kc ≃ 3.1 10 -2 . Dashed line: D = 10 0 , kc ≃ 8.4 10 -1 . Dotted line: D = 10 -3 , kc ≃ 1.2 10 1 . Solid line: D = 10 -6 , kc ≃ 1.3 10 2 . (b) Dimensionless angular frequency ω as a function as of k, Eq. 2.25. Note:

Fig. 2 . 6

 26 Fig. 2.6 Cut-off wavenumber kc (a), maximum wavenumber km (b), maximum angular frequency ωm (c) and maximum temporal growth rate σm (d) as functions of the bending rigidity D. Circles: data obtained numerically from the equations. Solid lines: asymptotic expressions (2.28-2.31).

Fig. 2 .

 2 Fig. 2.7 (a) Experimental values of the angular frequency as a function of the rescaled wind velocity. The best linear fit for the paper sheet (triangles) is shown in solid line (slope 1.9). The dotted line is the best linear fit (slope 2.6) for the plastic sheets (circles for the thick sheet (58 µm) and squares for the thin one (30 µm)). (b) Experimental values of the wavenumber as a function of the rescaled wind velocity. The best fit by a 2/3-power law gives a multiplicative factor 0.78. All these data are for ∆L = 4 cm.

Fig. 2 . 8

 28 Fig.2.8 Same as Fig.2.7, for the paper material only but with different sheet lengths: ∆L = 4 cm (circles), ∆L = 8 cm (squares) and ∆L = 12 cm (triangles). The solid lines are the same as in Fig.2.7 and correspond to the limit of small ∆L (linear regime).

Fig. 2 .

 2 Fig. 2.10 Rescaled longitudinal profiles of the wave amplitude. Symbols: red circles (∆L = 4 cm) green circles (∆L = 8 cm) and blue circles (∆L = 12 cm). All these data correspond to the paper sheet. Data for different wind velocities have been gathered.

Fig. 2 .

 2 Fig. 2.11 A piezoelectric flag in a uniform axial flow and the equivalent circuit of a piezoelectric pair connected with a parallel RL circuit [99].

  -Gerasimenko (abbreviated as 67P) is a Jupiter-family comet, originally from the Kuiper belt. It was first observed on photographic plates in 1969 by Soviet astronomers Klim Ivanovych Churyumov and Svetlana Ivanovna Gerasimenko, after whom it is named. As in Fig.3.1a, 67P has two lobes, a small one and a large one, with a thick neck connecting in between. The current orbit is shown in Fig.3.1b, and the latest fly-by over the perihelion was on 13 August 2015. There are 19 distinct regions on 67P, with each named after an Egyptian deity[START_REF] Thomas | The morphological diversity of comet 67p/churyumov-gerasimenko[END_REF]. Parameters of 67P and its current orbit are listed in Table3.1.67P was the destination of the European Space Agency (ESA) Rosetta mission (Fig.3.1b). Rosetta is a space probe launched on 2 March 2004 from the Guiana Space Centre in French Guiana. On 6 August 2014, it reached 67P, performed a series of manoeuvres and entered orbit on 10 September 2014. Rosetta's lander, Philae (Fig.3.1c), touched down on its surface on 12 November 2014, becoming the first spacecraft to land on a comet nucleus.Along with Philae, Rosetta is performing detailed investigations of 67P, and the mission continues to return data from the spacecraft in orbit and from the lander in the comet's surface as of 2015. On 30 September 2016, the Rosetta spacecraft ended its mission by landing on the comet in its Ma'at region.

Fig. 3 .

 3 Fig. 3.1 (a) Comet 67P by Rosetta's OSIRIS narrow-angle camera on 3 August 2014 from a distance of 285 km, with he image resolution of 5.3 meters/pixel. (b) Orbit of Comet 67P and the relative paths of the Rosetta probe and it. (c) Philae touched down on its surface. Credit: ESA/Rosetta/MPS.

3 . 3 )Fig. 3 .

 333 Fig. 3.2 (a) View of the comet's neck (Hapi) region by OSIRIS narrow-angle camera dated 18 September 2014, i.e. long before perihelion, showing evidence of local gas-driven transport producing dune-like ripples, with a wavelength λ ≃ 5 m. Superimposed yellow marks: position of the ripples from a photo dated 17 January 2016, i.e. long after perihelion (Fig. S1) providing evidence for their propagation. (b) Photograph of ripplemarks in 'Maftet' region. Image taken on 05 March 2016, when Rosetta was 20.3 km from 67P, with a resolution of 0.36 m/pixel. (c) Photograph of ripples at 'Hatmehit' region. Image taken on 13 April 2016, when Rosetta was 109.2 km from 67P, with a resolution of 1.98 m/pixel. Credit: ESA/Rosetta/MPS.

Fig. 3 . 3

 33 Fig. 3.3 Picture of the comet and its close coma. Red line shows the contour of the comet. Green line shows the contour of the vapor halo at the resolution of the instrument. Image taken on 18 February 2016, when Rosetta was 35.6 km from the comet, with a resolution of 3.5 m/pixel. Photo credit: ESA/Rosetta/MPS.

  is the solar radiation flux received by the comet at latitude ϕ, where ψ E ≃ 1360 W/m 2 is the radiation flux received from the sun at η E = 1 astronomical unit (au). η is the heliocentric distance of the comet, which is a known function of time along the comet's orbit. φ encodes the day-night alternation following φ(t) = max[cos(2πt/Γ d ), 0]. The heat flux, computed at the comet's surface by

  .18) k B = 1.3810 -23 J/K is the Boltzmann constant, d w ≃ 0.34 nm is water molecule size and m ≃ 3 10 -26 kg is the mass of a water molecule. Similarly, the averaged heat flux writes:

Fig. 4 . 2

 42 Fig.4.2 Time evolution of the vapour temperature (left axis) and corresponding thermal velocity (right axis) just above the comet's surface, calculated along the comet's orbit around the sun (Fig.3.1b). Time is counted with respect to the zenith, at perihelion. Bold orange lines: envelopes of the daily variations (inset) at perihelion, emphasising the maximum and minimum values. Inset: Zoom on the time evolution during one comet rotation at perihelion. The day/night alternation is suggested by the background grey scale.

Fig. 4 . 3

 43 Fig. 4.3 Time evolution of the vapour density (left axis) and corresponding mean free path (right axis) just above the comet's surface. Curve conventions are the same as in Fig. 4.2.

Fig. 4 .

 4 Fig. 4.4 (a) Time evolution of the vapour pressure just above the comet's surface, calculated along the comet's orbit around the sun. Curve conventions are the same as in Fig. 4.2. (b) Schematic of the outgassing process (blue) and the resulting winds (red arrows) driven by strong pressure gradients from illuminated to shadow areas.

Fig. 4 .

 4 Fig. 4.5 (a) Outgassing flux qm as a function of the comet's heliocentric distance η. Solid line: prediction of the model. Symbols: data from the literature, as specified in legend. Three symbols are used for data from Bertaux et al. (2014) corresponding to three different years: star (2009), inverted triangle (2002) and circle (1996). (b) Vapor flux at one single comet day (perihelion), where θ = 0 corresponds to the noon. Therefore it is night for 67P for θ = π/2 ∼ 3π/2, when q m is negligible.

Fig. 5 . 1

 51 Fig. 5.1 Estimation of the grain size on the surface of 67P. A Auto-correlation function C(δ) (red circles) com-puted from the photograph of the comet's granular bed, taken by Philae just be-fore its touch down in the Ma'at region (B), where large boulders and rocks have been excluded. The resolution of the picture is 9.5 mm/pixel. Photo credit: ESA/Rosetta/Philae/ROLIS/DLR. The correlation is compared to that computed with pictures of calibrated aeolian sand from the Atlantic Sahara (green square, lower axis, is expressed in units of the grain diameter) taken in the laboratory (Methods). The best collapse of the correlation functions is obtained for a mean grain diameter d ≃ 9.7 mm on the comet. C Histogram of grain size d computed from the photo-graph of the comet's granular bed shown in panel D taken by Rosetta just before its touch down in the Ma'at region. The best fit by a log-normal distribution, shown in red, gives a mean grain diameter d ≃ 38 mm.

Fig. 5 . 3

 53 Fig. 5.3 Dependence of the threshold shear velocity u t on the grain diameter d. The best fit of experimental measurements (symbols) by theoretical predictions gives a cohesive diameter d m ≃ 10 microns. Data from Yalin & Karahan, Hydraul. Div., Am. Soc. Civil Eng 105, 1433 (1979) [147].

Fig. 5 .

 5 Fig.5.4 (a) Time evolution of the velocity ratio u * /u t , calculated along the comet's orbit around the sun. Time is counted with respect to the zenith, at perihelion. Bold orange lines: envelopes of the daily variations (inset), emphasizing the maximum and minimum values. Inset: Zoom of the evolution of u * /u t during one comet day, at perihelion. The day/night alternation is suggested by the background grey scale. Wind is above the transport threshold in the afternoon (counted positive, schematized here and in next figures by a yellow circle with an arrow pointing downwards) and in the morning (counted negative, yellow circle with an arrow pointing upwards). (b)|u * | (black curve) and u t (blue curve) at perihelion.

Fig. 5 . 5

 55 Fig.5.5 Dependence of the threshold shear velocity u t with the grain diameter d at perihelion, for afternoon conditions. The minimal velocity above which sediment transport takes place is computed from the force balance on a grain between hydrodynamic drag, bed friction and Van der Waals cohesive forces. The threshold increases above d = 1 m due to gravity and below d = 1 mm due to cohesion. In between, u t is almost constant and on the order of 50 m/s due to the large mean free path of the vapor ℓ ≃ 3 cm (Fig.4.3). Yellow mark: range of observed grain sizes (Fig.5.1).

[ 153 - 156 ,Fig. 5 . 7

 15315657 Fig. 5.7 Saturation length L sat in units of d as a function of the flow velocity at a grain size above the surface u(d) rescaled by the grain settling velocity V fall . Data are obtained for various experimental conditions: grains in oil (circles), in water (squares), and in water-glycerin solution (triangles); color codes for the grain size from 100 microns (red) to 830 microns (violet). Black solid line: L sat /d ≃ 24.

5 10 4

 4 s. This time must be compared to the total time during which sediment transport takes place during a revolution around the sun, which is around 10 6 s (0.7% of the revolution period), i.e. 20 times larger. The ripples therefore have enough time to emerge and mature during one comet revolution.At the neck (Hapi) region, the observed emergent ripple wavelength λ is around 7 m. This corresponds to a grain size d = 4 mm. The sediment flux estimated from (5.20) is on the order of 4 10 -6 m 2 /s. The dispersion relation (Fig.6.4) indicates that the most unstable

Fig. 6 . 2

 62 Fig.6.2 Relaiton betwwen λ/L sat and the rescaled saturated length L sat u * /ν predicted from the most unstable mode (red curve), and the approximated scaling λ/L sat ≈ (L sat u * /ν) -2/5 ( black curve).

Fig. 6 . 4

 64 Fig. 6.4 Dispersion relation: dimensionless growth rate (a) and propagation speed (b) as functions of the rescaled wavenumber kν/u * , computed at perihelion for d = 4 mm, with a saturation length L sat /d = 24. This corresponds to the neck (Hapi) region, where the observed emergent ripple wavelength λ is around 7 m. The corresponding most unstable mode (red dot) is at kν/u * ≃ 0.06. Vapor viscosity and shear velocity are respectively ν ≃ 5 m 2 /s and u * ≃ 70 m/s, respectively. With a reference sediment flux ≃ 4 10 -6 m 2 /s, the growth rate of this mode is σ m ≃ 5.2 10 -3 (u * /ν) 2 ≃ 5 10 -6 s -1 . Mature ripples at a wavelength of 18 m (kν/u * ≃ 0.024) propagate at a velocity c ≃ 0.18 u * /ν ≃ 10 -5 m/s, i.e. over ≃ 10 m for the time during sediment transport occurs ≃ 10 6 s.

Fig. 7 .

 7 Fig. 7.1 (a) A global view photo, created by four images taken 450,000 km away from the Pluto by the New Horizons spacecraft on July 14 2015, showing features as small as 2.2 km. Credit: NASA/JHUAPL/SwRI. (b) New Horizons spacecraft from NASA's New Frontiers program. (c) Pluto's orbit (Yellow curve) and New Horizons spacecraft's trajectory (red curve). (b) and (c) are from the Web.

Fig. 7 . 2

 72 Fig. 7.2 Phase diagrams of four geologically relevant substances: nitrogen, ethane, carbon dioxide and water. Shaded regions indicate where the liquid phase is stable. Also shown are the surface pressure and temperature conditions for Pluto, Titan, Mars and the Earth. It is clear that liquids are expected at the surfaces of Titan and Earth, but not Mars or Pluto. Conversely, on Mars and Pluto, direct conversion from solid to gas (sublimation) or vice versa is likely. Figure from Schenk and Nimmo (2016) [164].

Fig. 7 .

 7 Fig. 7.3 (a) A highest-resolution image showing that the mountains end abruptly at the shoreline of the Sputnik Planum, where the soft, nitrogen-rich ices of the plain form a nearly level surface, broken only by the fine trace work of striking, cellular boundaries and the dune-like patterns on the surface of the plain's ices. This view is about 50 miles wide, and the top of the image is to Pluto's northwest. Credit: NASA/JHUAPL/SwRI. (b) An extended colour image shows the rounded and bizarrely textured mountains, where some crescentic dunes can be observed (framed by a red rectangle). Image was taken by the Ralph/Multispectral Visual Imaging Camera (MVIC) on July 14, 2015, with a resolution around 510 m/pixel.

Fig. 8 . 1

 81 Fig. 8.1 Departure shot of Pluto by New Horizons, showing Pluto's atmosphere backlit by the Sun. The blue color is close to what a human eye would have seen, and is caused by layers of haze in the atmosphere.

Fig. 8 . 2

 82 Fig. 8.2 Pressure and temperature in Pluto's lower atmosphere. (a) Pressure. (b) Temperature. These profiles were retrieved from radio occultation data recorded by the REX instrument onboard New Horizons. Each graph shows results at both entry (red line with circles) and exit (blue line with triangles), situated on opposite sides of Pluto. The profiles are most accurate at the surface, where the uncertainties in pressure and temperature are 1 mbar and 3 K, respectively. Temperature fluctuations at altitudes of >20 km are caused by noise; no gravity waves were detected at the sensitivity of these measurements. The dashed line indicates the saturation temperature of N 2 . Figure from [166].

1 ) 8 . 2 Thermo-hydrodynamics of Pluto's atmosphere 8 . 2 . 1

 182821 with A 0 = 1.24 10 1 , A 1 = -8.074 10 2 , A 2 = -3.926 10 3 , A 3 = 6.297 10 4 , A 4 = -4.633 10 5 , A 5 = 1.325 10 6 (Table5and Figure21in[START_REF] Fray | Sublimation of ices of astrophysical interest: A bibliographic review[END_REF]). One therefore obtains p sat = 8 Pa for T 0 = 40 K, and the corresponding saturated mass density ρ sat = mP sat /(k B T ) = 6.8 10 -4 kg/m 3 . Thermal processes Following the procedure in Section 4.2, we compute the thermal diffusion in the nucleus of Pluto, and find a diurnal thermal penetrating length δ d ≃ 0.2 m, which means that a few tens of cm below the surface, the day-night alternation has no influence on the temperature field, and a seasonal penetrating length δ y ≃ 22 m. Energy balance at the Pluto surface is governed by Eq. 4.8. Neglecting the sublimation due to the almost saturated state, it simplifies into:(1 -Ω N )ψ = ΣεT 4 0 + J 0 . (8.2) Together with Eqs. 4.10 and 4.11, Eq. 8.2 is solved at the revolution scale of Pluto, including the day-night alternace. As you could see in Fig. 8.3, the computed temperature is around the measured 40 K, with very small day-night variation.

Fig. 8 . 3

 83 Fig. 8.3 Time evolution of the surface temperature of Pluto, calculated along the Pluto's orbit around the sun (Fig. 7.1c). Time is counted with respect to the zenith, at perihelion. (A) Envelopes of the daily variations at the revolution scale, emphasising the maximum and minimum values. (B) Zoom on the time evolution during one Pluto rotation at perihelion.

Fig. 8 . 4

 84 Fig. 8.4 Time evolution of the atmosphere pressure near Pluto's surface. We take p 0 = p sat (T 0 ). Curve conventions are the same as in Fig. 8.3.

Fig. 8 . 5

 85 Fig. 8.5 Time evolution of the atmosphere density near Pluto's surface. We take ρ 0 = ρ sat (T 0 ). Curve conventions are the same as in Fig. 8.3.

Fig. 8 . 6

 86 Fig. 8.6 Time evolution of the the friction velocity u * at Pluto's surface. Curve conventions are the same as in Fig. 8.3.

Fig. 9 . 2 A

 92 Fig. 9.2 A given point x receives light from a portion of the surface. The rays determining the limits of this portion are either tangential to the surface at the point considered (a) or at the point of emission (b). These conditions determine x a and x b by Eqs. 9.21-9.23. Figures from [67].

Fig. 9 .

 9 Fig. 9.3 (a) Bounds of the integral (9.20) given by Eqs. 9.21-9.23. (b) Integral I(η) (solid line) giving the illumination profile for a sinusoidal surface, compared to the function 1 -cos(k x) (dotted line). Figures from [67].

  58 and 9.25, it gives

( a ) 3 .

 a3 The particluar solution is computed, (b) The general solution is computed associated to the missing variables at the bottom: S t (0), S n (0) and Y (0); Similarly to Step 2, We cross downward from the top to the center of the domain to find another solution. Firstly, the particular solution is computed, and then the general solution is computed by adjusting the missing variables at the top: U(η D ), S n (η D ) and Θ(η D ); 4. Then we match the two solutions at the center of the domain. The matched solution satisfies the boundary conditions at both sides of the domain, from which we are able to compute the dispersion relation and other secondary parameters interested.

= 37

 37

3 10 - 6 9 . 4 .

 694 and send other effects to 0: = 0, -= 0 and Σ = 0. The dispersion relation and instability diagram are computed by varying + , as shown in Fig. The results indicate that heat flux has a destabilising effect and induces instabilities in the range where ∈ (10 -4 , 10 -1

Fig. 9 .

 9 Fig. 9.4 (a) Dispersion relation computed at different + . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. = 3.3 10 -6 , = 0, Σ = 0 and = 0 for all the curves; and + = 5.5 10 -14 for black ones, + = 3.0 10 -12 for blue ones, + = 1.6 10 -10 for red ones. Dotted line shows the position of the maximum growth rate. (b) Stability diagram in + --1 plane. Dashed line is the marginal stability curve, and solid curve indicates the position of the most unstable mode.

3 10 - 6 ,

 6 + = 3.0 10 -12 , -= 0 and Σ = 0. The dispersion relation is then computed varying , as shown in Fig. 9.7. The results indicate that the non-zero creates a new range of unstable wave numbers at large wavelength where the growth rate is constant. The cut-off wavenumber increases linearly with the blowing velocity . We can see that the dispersion relation keeps the same shape at large wave numbers, at least when is small (Fig. 9.7a).

Fig. 9 . 5

 95 Fig. 9.5 Perturbed fields showing the instability due to heat diffusion and convection. Diffusion coefficient field (a), velocity field (b) and temperature field (c) are plotted both in iso-coutours and in colours. The values are decreasing from red to violet for both colours and iso-coutours.

Fig. 9 .

 9 Fig. 9.6 (a) Dispersion relation computed at different . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. + = 3.0 10 -12 , -= 0, Σ = 0 and = 0 for all the curves; and = 3.3 10 -7 for black ones, = 3.3 10 -6 for blue ones, = 3.3 10 -5 for red ones. (b) Stability diagram in --1 plane. Dashed line is the marginal stability curve, and solid curve indicates the position of the most unstable mode.

Fig. 9 . 7

 97 Fig. 9.7 Dispersion relations computed at different . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. + = 3.0 10 -12 , -= 0, Σ = 0 and = 3.3 10 -6 for all the curves. In panel (a) = 0 for black ones, = 10 -6 for red ones, = 10 -5 for blue ones, = 10 -4 for green ones; in panel (b) = 10 -3 for red ones, = 10 -2 for blue ones, = 10 -1 for green ones.

Fig. 9 . 8 Fig. 9 . 9

 9899 Fig. 9.8 Dispersion relation computed at different Σ . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. + = 3.0 10 -12 , = 3.3 10 -6, -= 0 and = 0 for all the curves; and Σ = 0 for grey ones, Σ = 10 -8 for red ones, Σ = 10 -7 for black ones, Σ = 10 -6 for violet ones, Σ = 10 -5 for blue ones, Σ = 10 -3 for green ones.

Fig. 9 .

 9 Fig. 9.10 (a) Dispersion relation computed at different . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. Σ = 0, ψ = 10 -10 , -= 0 and = 0 for all the curves; and = -5.5 10 -5 for orange ones, = 6.9 10 -5 for green ones, = 7.6 10 -5 for blue ones, = 8.3 10 -5 for red ones. (b) Stability diagram in --1 plane. Blow the dotted line, it is scale-free unstable zone. Dashed line is the marginal stability curve, and solid curve indicates the position of the most unstable mode.

Fig. 9 .

 9 Fig. 9.11 Dispersion relations computed at different . Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. ψ = 10 -10 , -= 0, Σ = 0 and = 7.6 10 -5 for all the curves; = 0 for red ones, = 10 -4 for black ones, = 10 -3 for blue ones, = 10 -2 for green ones.

Fig. 9 .

 9 Fig. 9.13 Dispersion relation computed at different -. Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite valueσ. ψ = 10 -10 , = 7.6 10 -5Σ = 0 and = 0 for all the curves. In panel (a), -= 0 for blue curves, -= 10 -12 for red ones, -= 10 -11 for green ones, -= 10 -10 for orange ones. In panel (b), -= 0 for blue curves, -= -4 10 -12 for violet ones, -= -6 10 -12 for grey ones, -= -8 10 -12 for green ones, -= -10 -11 for red ones and -= -10 -10 for orange ones.

Fig. 9 .

 9 [START_REF] Lighthill | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF], dispersion relations for different flux distributions are computed with = 3.3 10 -6 , corresponding a wind velocity u * = 0.1 m/s . The value ε = 0 is the pure solar radiation situation. One obtains instability over all the wave number without selection (black curve in Fig.9.15). For this kind of scale-free dispersion relation, we expect a nonlinear pattern coarsening, which will lead to larger and larger wavelengths[START_REF] Politi | When does coarsening occur in the dynamics of onedimensional fronts?[END_REF]. This is probably not the case on Pluto. Moreover, it can not explain the direction of the patterns changing from place to place (Fig.9.14). If a small fraction of the heat flux comes from the atmosphere ( ε = 0.25), one can see the instability being selective in wave number,

Fig. 9 .

 9 Fig. 9.15 Dispersion relation computed at different ε for Pluto. Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite value.= 3.3 10 -6 for all the curves, ε = 0 for black one, ε = 0.25 for orange ones, ε = 0.5 for blue ones, ε = 0.75 for green ones, and ε = 1 for red ones. The dotted line indicates the locating most unstable mode.

Fig. 9 .

 9 Fig. 9.16 Dispersion relation computed at different for Pluto. Solid curves and dashed ones are respectively the rescaled growth rate σ and its opposite value. ε = 0.5 for all curves, = 3.3 10 -6 for blue ones, = 3.3 10 -5 for red ones, = 3.3 10 -4 for green ones, and = 3.3 10 -3 for black ones.

Fig. 10 . 1

 101 Fig. 10.1 Selection of the ripple wavelength and propagation speed. (a) Ripple wavelength λ as a function of the wind shear velocity u * . (b) Ripples propagation speed c as a function of u * . Experimental data are from[START_REF] Andreotti | Aeolian sand ripples: experimental evidence of fully developed states[END_REF], and the simulation ones are from[START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF].

Fig. 10 . 2 A

 102 Fig. 10.2 A particular grain trajectory over a flat (a) and a rippled (b and c) bed. Red curves are the trajectories and the black solid curves are the profile of the bed surfaces. The reference flat bed surface are noted as dashed black lines in (b) and (c).

(Fig. 10 . 3

 103 Fig. 10.3 Modulus (a) and phase (b) of the factor involved in Eq. 10.2 as functions of ℓ/λ. Solid line: geometrical contribution (Eq. 10.4), taking into account the small variations of α ↓ . Dotted line: contribution of the modulated wind. Dashed line: sum of these two contributions. Data from the numerical simulation by Orencio Durán with the model similar to [131].

Fig. 10

 10 Fig. 10.4 Modulus (a) and phase (b) of flux modulation rate . Dots are computed from the full simulations[START_REF] Durán | Direct numerical simulations of aeolian sand ripples[END_REF], with u * /u t = 4, and with increasing wave numbers k labelled from red to violet. Analytical prediction (Eq. 10.18) corresponding to the geometrical contribution, with an adjustable overall prefactor on the order of 2.

  condition P(ℓ)dℓ = 1, this leads to the following constrain on : (s) ds = 0. (10.23)

  ) (kℓ) e -ikℓs (s) dsdℓ.(10.24) where we have changed variable and taken s = (xx ↓ )/ℓ. With the expression of P(ℓ) (Eq. 10.1) and changing variable with r = kℓ, we can write Eq. 10.24 as (s) e -irs dsd r.(10.25) Considering the definition of Fourier transform, one further obtains,

  e -irs dsd r = -ik ln (ℓ s /ℓ d )

Fig. 10 . 5

 105 Fig. 10.5 Complex growth rate Ω d as a function of rescaled wave number k ℓ s . Computed from Eq. 10.26 withgiven by (10.18) and given by (10.29), and with a = 1, s = 0 and ℓ d = 0.

2 Fig. 10 . 6

 2106 Fig. 10.6 Growth rate Ω s as a function of rescaled wave number k ℓ s /b. Computed from Eq. 10.31 with s = 0 and ℓ d = 0.

Fig. 10 . 7

 107 Fig. 10.7 Dispersion relation measured from numerical simulations. Ripple growth rate as a function of the wavenumber k for different wind speeds (red: u * /u th = 3, green: u * /u th = 4 , blue: u * /u th = 5 ). The solid line is the best fit by Eq. 10.34. Figure from [131].

Fig. 10 . 8

 108 Fig. 10.8 Experimental evidence of a modulated saltation flux. Images obtained after analysis of a fast movie recorded in a wind tunnel. The ripple is vertically illuminated with a laser sheet. Images are corrected from the heterogeneity of the light intensity. Figure from [131].

Fig. 10 . 9

 109 Fig. 10.9 Sketch of the experimental setup for aeolian sand ripples. The sand tank and the wind tunnel are connected to the downstream air flow, induced by an axial fan through a honeycomb and a contraction. Sand bed is labelled in dirt-colored. Dotted dashed line shows the central line of the tunnel. x, y and z are denoted as the longitudinal, transverse vertical and horizontal directions, respectively.

Fig. 10 .

 10 Fig. 10.10 Aeolian sand ripples emerging from a flat bed in the wind tunnel, with the arrow showing the wind direction. Wind velocity at the central of the tunnel is V ≈ 4.0 m/s

2 .

 2 The development of ripples on 67P can be ascribed to the outgassing across a porous surface granular layer, as well as the surface vapor flows induced by the vapor pressure gradient due to the day-night alternace. Most vapour is emitted from sub-surface ice and flows out at a much smaller velocity, resulting in the presence of a thicker atmosphere on the comet surface. The surface flow estimated by the turbulent boundary layer theory is responsible for the sediment transport for the grains at the centimeter scale. Expected temporal and spatial scales are obtained by making the cometary case an analogue to the terrestrial subaqueous bedload. More generally, bedforms reflect the characteristics of the bed and the flow they originate from, making them a good indicator of the physical mechanisms at work. Comets thus provide an opportunity to better understand erosion and accretion processes on planetesimals, with implications for the open question of how these bodies can grow from the meter to the kilometer scale[START_REF] Johansen | Rapid planetesimal formation in turbulent circumstellar disks[END_REF][START_REF] Jansson | Formation of pebble-pile planetesimals[END_REF].

Cd

  Drag coefficient. C s Heat capacity of the solid ice. C ν Heat capacity of the vapour. C Bulk heat capacity of 67P's nucleus. D Bending rigidity of flexible sheets. E Grain Young modulus.F drag Drag force exerted on a grain.

F↑FL

  Dimensionless function related to the first order correction of phi. Ice surface fraction. I Effective thermal inertia of 67P's nucleus. L 0 Distance between two masts. L sat Saturated length. Mixing length. M P Mass of Pluto. M c Mass of comet 67P. N The number of transported grains per unit area. P(ℓ) Distribution of the hop length. R P Radius of Pluto. R c Effective spherical raduis of comet 67P. R l Effective spherical raduis of large lobe comet 67P. R n Raduis of the neck of comet 67P. R s Effective spherical raduis of small lobe comet 67P. R 0 t Transitional Reynolds number for the homogeneous case. R t Transitional Reynolds number.S c Surface area of comet 67P.

Sn

  Dimensionless function related to the first order correction of τ zz -p. S t Dimensionless function related to the first order correction of τ xz . S x x Dimensionless function related to the first order correction of τ x x .

WYℓηγκλ

  Dimensionless function related to the first order correction of z velocity. Dimensionless function related to the first order correction of heat flux in z direction. ∆ Crater function. Γ d Rotation period of comet 67P, and Pluto. Γ y Revolution period of comet 67P, and Pluto. Ω Complex growth rate. Θ Dimensionless function related to the first order correction of temperature field. Υ 0 Ratio of u 0 and V 0 th . α ↓ Arrival angle of the grain on the bed. D Dimensionless bending rigidity. ω Dimensionless angular frequency. σ Dimensionless growth rate. k Dimensionless wave number. β t Turbulent Prandtl number. β ν Viscous Prandtl number. δ i Thickness of the inner turbulent layer. Hop length. Heliocentric distance of 67P. Adiabatic expansion coefficient.γ Surface tension of the grain material.κ c Thermal diffusivity of 67P's nucleus. Karman constant. Wave length.

ν t Turbulent viscousity.νωψστcdgk

  Molecule viscousity. Angular frequency. φ b Bed volume fraction. Distribution of grains volume arriving per unit time in a interval with a certain hop length. ρ i Density at the sublimation interface. ρ 0 Base state density. ρ P Bulk density of Pluto. ρ c Bulk density of comet 67P. ρ sat Saturated density. Growth rate. Basal shear stress. ϵ Specific energy. ϕ Erosion rate. ξ Solid ice surface profile. ζ Two-dimensional shape of the sheet. Propagation speed. d m Cohesive size. Grain size. Gravity acceleration. k B Boltzmann constant. k s Thermal conductivity of the solid ice. k ν Thermal conductivity of the vapour. Wave number. p i Pressure at the sublimation interface. p 0 Base state pressure. p sat Saturated vapour pressure. q(x, t) Particle flux. q sat Saturated particle flux. r a Aphelion distance of comet 67P, and Pluto. r p Perihelion distance of comet 67P, and Pluto. u * shear velocity. u r Radial velocity components in inner layer. u t Threshold velocity. u x Vapour velocity in x direction.u z Vapour velocity in z direction. u θ Azimuthal velocity components in inner layer.

Table 2 .

 2 1 Values of the mass per unit surface m and the bending rigidity D of the three different sheets used in the experiments.

		Material	m	D	h
				kg/m 2	Nm	µm
		Paper	2.5 × 10 -1 7.6 × 10 -3 90
		Plastic	8.6 × 10 -2 1.5 × 10 -4 58
		Plastic	4.9 × 10 -2 4.2 × 10 -5 30
	are given in Tab. 2.1, namely the mass per unit surface m and the bending rigidity D. This
	last quantity was measured by determining the deflection under gravity of horizontally
	clamped strips of various lengths. Interestingly, D changes by two orders of magnitude
	from the thin plastic sheet to the paper, which allows us to investigate the instability over a
	large range of parameters. A fast camera (Phantom Miro M340, Appendix A) is employed
	to record the motion of the sheet through the transparent ceiling of the tunnel. The camera
	is operated with a spatial resolution of 2560 pixel × 320 pixel, at 1.5 pixel/mm, and an
		Y
		O	X	Fast Camera
		Z			Control Laptop
	Air Flow
	V m/s	Wind Tunnel
		Mast			Flexible Sheet	Mast

exposure time of 400 µs. The sample rate is typically between 1000 Hz and 2500 Hz during the measurements.

Table 3 .

 3 1 Parameters of 67P and its current orbit.

	Item	67P
	Mass M c	1.0 × 10 13 kg
	Bulk density ρ c	470 kg/m 3
	Rotation period Γ d	12.4 hours
	Revolution period Γ y	6.44 years
	Perihelion distance r p	1.24 AU
	Aphelion distance r a	5.68 AU

  .[START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF] We finally solve (4.44) and (4.45) to obtain u * as well as δ i . Note that the above equations are only valid if the thickness of the turbulent boundary layer is larger than that of the viscous sub-layer, i.e. when δ i ≳ 10ν/u * .

  56. Its value is remarkably insensitive to p c , and is always around Υ 0 = 0.11, corresponding to a Mach number ≃ 0.15. So far, we have presented the full description of the coma's thermo-hydrodynamics, and the problem is then implemented in an ideal spherical geometry, considering solar radiation modulated by both the seasonal scale and the day-night alternation. Eqs. 4.5, 4.6, 4.8, 4.10, 4.11, 4.17, 4.51, 4.54 and 4.55 are solved iteratively for thermodynamics. It should be noted that the saturated vapour density ρ sat is only the function of the temperature and

Table 6 .

 6 1 Ripple crest-to-crest distance measured on pictures of different regions of 67P. N + 1 is the number of successive ripple crests identified on the picture. λ is the average value of their distance. t is the time to perihelion(13 Aug. 2015), counted positive (negative) after (before) it.

	Photo name

Table 7 .

 7 1 Parameters of Pluto and its orbit.

	Item	Values
	Mass M P	1.3 × 10 22 kg
	Bulk density ρ P	1885 kg/m 3
	Radius R P	1185 km
	Surface gravity g	0.625 m/s
	Escape velocity v e	900 m/s 2
	Rotation period Γ d	153.282 hours
	Revolution period Γ y	247.74 years
	Perihelion distance r p	29.6 AU
	Aphelion distance r a	49.3 AU
	Perihelion solar radiation ψ p	1.6 W/m 2
	Aphelion solar radiation ψ a	0.6 W/m 2

Table 8 .

 8 1 Parameters of nitrogen (N 2 ) and nitrogen ice (N s 2 ) .

	Item

  2 ) i 1 + e -ir + 2a 3 a π 1 -e -a 2 [1 + 2a 2 ] (1 -erf(ia)) -i4a 2 (10.30)

				e -ir -1 r	d r
	=	π 2a π 2	erf(a) + ia 2	2 F 2 1, 1; 3/2, 2; -a 2
		+ s	π

3

a π 1 + e -a 2 (1 -erf(ia)) -

Table A .

 A 2 Technical data of fast camera Phantom Miro M340

	Item	Value
	Resolution	2560 × 1600
	Sensor	25.6 mm × 16.0 mm
	Frame rate	800 fps
	Operation temperature	10 ∼ -40 • C
	Size	190 mm × 90 mm × 100 mm
	Power supply	100 -240 VAC power supply, Rechargeable battery
	Weight	1.4 kg (without CineFlash, lens or battery)

appliqué à la récupération d'énergie.La deuxième partie porte sur une analyse théorique de la formation de rides géantes sur la comète 67P, récemment observées par la sonde Rosetta. Nous montrons comment le dégazage de vapeur se produit au travers d'une couche poreuse granulaire superficielle et comment l'alternance jour/nuit conduit à des gradients de pression gigantesques qui engendrent des vents thermiques de surface. Ces motifs apparaissent comme étant les analogues de rides qui se forment à la surface de lit sableux dans un écoulement visqueux. L'analyse de stabilité linéaire du problème permet de prédire quantitativement l'émergence de ces rides à la longueur d'onde et à la vitesse de propagation observées. Cette description fournit un outil robuste et fiable pour décrire les processus d'érosion et d'accrétion dans l'évolution des petits corps. Dans la troisième partie, nous proposons un modèle pour l'apparition de motifs de sublimation sur Pluton, tels que ceux observés sur Sputnik Planum. La formation et l'évolution de ces motifs proviennent de la sublimation/condensation différentielle de la glace d'azote. Nous montrons que l'atmosphère de Pluton possède des propriétés (température et pression) peu variables en espace et en temps. Nous analysons les différents mécanismes d'instabilité en compétition et concluons à un mécanisme original, basé sur le mélange et le transport de chaleur dans l'atmosphère, plutôt qu'au mécanisme des pénitents, basé sur l'auto-éclairement de la surface de glace. Enfin, nous avons étudié théoriquement l'instabilité de formation des rides éoliennes en considérant les trajectoires des grains résonantes avec le relief. Cette modélisation prend en compte de manière simple et effective les effets collectifs du transport de sédiments. Le modèle est validé à partir de simulations numériques existantes, elles mêmes calées sur des expériences contrôlées.

Acknowledgements

I would like firstly to thank my supervisors, Profs. Bruno Andreotti and Philippe Claudin, for their unwavering support throughout my PhD process in PMMH lab at ESPCI. They

Pluto and the New Horizons mission

Pluto is a dwarf planet in the Kuiper belt, a ring of bodies beyond Neptune (Fig. 7.1a), discovered by the American astronomer Clyde Tombaugh in 1930. It was originally considered the ninth planet from the Sun, yet reclassified as a member of the new "dwarf planet" category in 2005. Pluto has a moderately inclined relative to the ecliptic (over 17°) and moderately eccentric (elliptical) orbit (Fig. 7.1c), which leaves a small region of Pluto's orbit lies nearer the Sun than Neptune's. Pluto's last perihelion transit was on May 8, 1989. Parameters of Pluto and its orbit are listed in Table 7.1.

Pluto was the primary mission to perform a flyby of NASA's New Horizons spacecraft, which is a part of NASA's New Frontiers programme. New Horizons was engineered by Johns Hopkins University and Southwest Research Institute in US, and launched on January 19 2006 (Fig. 7.1b,c). After nearly ten years, three-billion-mile journey, the spacecraft began its approach phase to Pluto on January 15, 2015, and it flew 12,500 km above the surface of Pluto on July 14 2015, making it the first spacecraft to explore the dwarf planet.

Having completed its flyby of Pluto, it is now on the way to a secondary mission to fly by and study one or more other Kuiper belt objects, expected to take place in January, 2019, when it is 43.4 AU from the Sun.

Rythmic patterns on the surface of Pluto

During the flyby, New Horizons made detailed measurements and observations of Pluto and its moons, and lots of investigations have been stimulated on atmosphere [START_REF] Singer | On the provenance of pluto's nitrogen (n 2 )[END_REF][START_REF] Gladstone | The atmosphere of pluto as observed by new horizons[END_REF],

on geology [START_REF] Moore | The geology of pluto and charon through the eyes of new horizons[END_REF], on space enviroment [START_REF]Pluto's interaction with its space environment: Solar wind, energetic particles, and dust[END_REF][START_REF] Mccomas | Pluto's interaction with the solar wind[END_REF], on surface compositions [START_REF] Dombard | Pluto's polygons explained[END_REF][START_REF] Grundy | Surface compositions across pluto and charon[END_REF],

Chapter 8

Pluto's atmosphere

In this chapter, the dynamics of Pluto's atmosphere is modelled in a simple way, to serve the basics/reference for the subsequent discussions on the physical modelling and on the estimation of parameters. More recent and detailed discussions can be found in [START_REF] Gladstone | The atmosphere of pluto as observed by new horizons[END_REF] and [START_REF] Bertrand | Observed glacier and volatile distribution on pluto from atmosphere-topography processes[END_REF]. In the following, we first give a general description, and then discuss the thermal and hydrodynamical properties of Pluto's atmosphere .

A general description

The atmosphere of Pluto is the tenuous layer of gases surrounding Pluto. Its existence has been studied since 1980s by way of earth-based observation of occultations of stars by Pluto [START_REF] Elliot | The recent expansion of pluto's atmosphere[END_REF][START_REF] Elliot | Changes in pluto's atmosphere: 1988-2006[END_REF][START_REF] Olkin | Evidence that pluto's atmosphere does not collapse from occultations including the 2013 may 04 event[END_REF] and spectroscopy [START_REF] Yelle | Pluto and Charon[END_REF], and was convincingly confirmed by observations in 1988 [START_REF] Stern | Pluto' in Encyclopedia of the Solar System[END_REF]. In 2015, it was revisited from a close distance by the spacecraft New Horizons [START_REF] Gladstone | The atmosphere of pluto as observed by new horizons[END_REF][START_REF] Stern | The pluto system: Initial results from its exploration by new horizons[END_REF]. From a recent observation by New Horizons (Fig. 8.1), we see the atmosphere is present all over the surface. This is because of the larger gravity, which is only 10 times smaller that that on the earth. More over the escape velocity is larger than the thermal velocity, which is on the order of 200 m/s, so there is almost no leak there, and atmosphere is close to the saturated state (Table 7.1). The main component of the atmosphere of Pluto is nitrogen (N 2 ), with minor amounts of methane (CH 4 ) and carbon monoxide (CO).

From the latest measurements by New Horizons, the content of methane is 0.25% [START_REF] Stern | The pluto system: Initial results from its exploration by new horizons[END_REF]; regarding the carbon monoxide, there are Earth-based estimates 0.05-0.075% in 2015 [START_REF] Gurwell | Detection of atmospheric co on pluto with alma[END_REF]. All kinds of gas are vaporized from their ices on Pluto's surface [START_REF] Stern | The pluto system: Initial results from its exploration by new horizons[END_REF][START_REF] Stern | Pluto' in Encyclopedia of the Solar System[END_REF], and we consider only nitrogen and nitrogen ice (termed as N s

2 ) in the following analysis due to the dominated proportion, whose parameters are listed in Table 8.1.

Chapter 9

Physical model of sublimation dunes

Most natural patterns form by linear instability, in which disturbances of different wavelengths are selectively amplified, such as aeolian dunes and subaqueous ripples on the Earth and other planets like Mars [START_REF] Charru | Sand ripples and dunes[END_REF][START_REF] Bourke | Extraterrestrial dunes: An introduction to the special issue on planetary dune systems[END_REF][START_REF] Claudin | Physical processes causing the formation of penitentes[END_REF]. Moreover, from the modeling of Pluto's atmosphere in the previous chapter, we see rather gentle variations in the parameters at both rotation and revolution time scales. Taking this two factors, we propose here a model for sublimation dunes based on the linear instability analysis by considering a flat solid surface which is slightly out of equilibrium with its vapour, no other phase in the atmosphere, submitted to a homogeneous and time independent wind. The analysis is two dimensional, with x the horizontal wind direction and the z-axis vertical, oriented upwards (Fig. 9.1).