Antoine Renaud

Louis-Philippe Nadeau

Antoine Venaille
email: antoine.venaille@ens-lyon.fr

Periodicity disruption of a model quasi-biennial oscillation

The quasi-biennial oscillation (QBO) of equatorial winds on Earth is the clearest example of the spontaneous emergence of a periodic phenomenon in geophysical fluids.In recent years, observations have revealed intriguing disruptions of this regular behaviour, and different QBO-like regimes have been reported in a variety of systems. Here we show that part of the variability in mean flow reversals can be attributed to the intrinsic dynamics of wave-mean flow interactions in stratified fluids. Using a constant-in-time monochromatic wave forcing, bifurcation diagrams are mapped for a hierarchy of simplified models of the QBO, ranging from a quasilinear model to fully nonlinear simulations. The existence of new bifurcations associated with faster and shallower flow reversals, as well as a quasiperiodic route to chaos are reported in these models. The possibility for periodicity disruptions is investigated by probing the resilience of regular wind reversals to external perturbations.

Earth's equatorial stratospheric winds oscillate between westerly and easterly mean flow every 28 months. These low-frequency reversals known as quasi-biennial oscillations are driven by high-frequency waves emitted in the lower part of the atmosphere, and supported by the presence of stable density stratification [START_REF] Baldwin | The quasi-biennial oscillation, Reviews of Geophysics[END_REF]. It is an iconic example of the spontaneous emergence of a periodic phenomenon in a turbulent geophysical flow [START_REF] Vallis | Atmospheric and oceanic fluid dynamics second edition[END_REF], with analogues in other planetary stratospheres [START_REF] Dowling | Planetary science: Music of the stratospheres[END_REF], in laboratory experiments [START_REF] Plumb | The instability of a forced standing wave in a viscous stratified fluid: A laboratory FIG. 3. Critical slowing down. a. Hovmöller diagrams of the mean-flow[END_REF][START_REF] Sémin | Nonlinear saturation of the large scale flow in a laboratory model of the quasibiennial oscillation[END_REF], as well as in idealized numerical simulations [START_REF] Wedi | Direct numerical simulation of the Plumb-McEwan laboratory analog of the QBO[END_REF][START_REF] Couston | Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves[END_REF]. In recent years, increasing attention has been given to the robustness of these regular reversals to external wave forcing and perturbations. Disruptions of this type of oscillations have been observed both in the Earth's atmosphere [START_REF] Osprey | An unexpected disruption of the atmospheric quasi-biennial oscillation[END_REF][START_REF] Newman | The anomalous change in the QBO in 2015-2016[END_REF] and in Saturn's atmosphere [START_REF] Fletcher | Disruption of Saturn's quasiperiodic equatorial oscillation by the great northern storm[END_REF]. In addition, a variety of oscillatory regimes, including non-periodic ones, have been reported in direct numerical simulations of a stratified fluid forced by an oscillating boundary [START_REF] Wedi | Direct numerical simulation of the Plumb-McEwan laboratory analog of the QBO[END_REF] or driven by an explicitly resolved turbulent convective layer [START_REF] Couston | Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves[END_REF]. Non-periodic oscillations have also been reported in global circulation model simulations of the solar interior and Giant planets [START_REF] Rogers | Angular momentum transport by gravity waves in the solar interior[END_REF][START_REF] Showman | Atmospheric Circulation of Brown Dwarfs and Jupiter and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations[END_REF].

Until now, the non-periodic nature of the reversals were interpreted as the system's response to transient external variations.

For example, the non-periodic disruption of the Earth's QBO and Saturn's QBO-like oscillation have been attributed to the response of equatorial stratospheric dynamics to extratropical perturbations [START_REF] Osprey | An unexpected disruption of the atmospheric quasi-biennial oscillation[END_REF][START_REF] Newman | The anomalous change in the QBO in 2015-2016[END_REF][START_REF] Fletcher | Disruption of Saturn's quasiperiodic equatorial oscillation by the great northern storm[END_REF]. Also, the existence of nonperiodic regimes in direct numerical simulations of stratified flows has been related to the time variability of the underlying turbulent convective layer [START_REF] Couston | Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves[END_REF][START_REF] Showman | Atmospheric Circulation of Brown Dwarfs and Jupiter and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations[END_REF]. Here, we show that the non-periodic nature of the reversals is a fundamental characteristic of stratified fluids by revealing the existence of a vast diversity of oscillatory regimes obtained using a simple steady monochromatic forcing. We further demonstrate that this rich intrinsic variability effectively controls part of the system's response to a transient external variation. Periodicity disruptions are more easily triggered and are increasingly lengthened when the system approaches a bifurcation point.

Model. The simplest configuration capturing the dynamics of the quasi-biennial oscillation (QBO, Fig. 1a) is given by a vertical 2D section of a stably stratified Boussinesq fluid, periodic in the zonal (longitudinal) direction, and forced by upward propagating internal gravity waves. This wave forcing is typically generated by an oscillating bottom boundary meant to represent the effect of tropopause height variations on the stratosphere. The evolution of the horizontally averaged zonal velocity, u, is governed by a simplified version of the momentum equation

∂ t u -ν∂ zz u = -∂ z u w , ( 1
)
where ν is the kinematic viscosity, z is the upward direction, and u w is the Reynolds stress due to velocity fluctuations around the zonal average. In weakly nonlinear regimes, this stress is carried by internal gravity waves, and any process damping the wave amplitude leads to a transfer of momentum from the waves to the mean-flow through the Reynolds stress divergence. Wave properties are also affected by the mean-flow and this interplay results in a complex coupled system.

To close the dynamical system, one needs to compute the Reynolds stress in [START_REF] Baldwin | The quasi-biennial oscillation, Reviews of Geophysics[END_REF]. In this study, the wave field is simulated either by taking into account all nonlinear interactions between waves and mean-flow (hereafter "nonlinear 2D model") or by considering a simplified closure that neglects wave-wave interactions, together with a WKB approach [START_REF] Plumb | The interaction of two internal waves with the mean flow: Implications for the theory of the quasibiennial oscillation[END_REF]. The latter approach (hereafter "quasilinear 1D model") has proven to be successful in explaining the spontaneous emergence of low-frequency periodic flow reversals [START_REF] Vallis | Atmospheric and oceanic fluid dynamics second edition[END_REF][START_REF] Plumb | The interaction of two internal waves with the mean flow: Implications for the theory of the quasibiennial oscillation[END_REF] and synchronisation with an external forcing [START_REF] Rajendran | Synchronisation of the equatorial QBO by the annual cycle in tropical upwelling in a warming climate[END_REF]. By assuming horizontally averaged dynamics, this quasilinear model is much simpler than the original flow equations but has nevertheless a large number of degrees of freedom since an infinite number of vertical oscillatory modes are possible. The 1D model is thus a natural starting point to investigate how periodic reversals are destabilized when the forcing strength is increased.

We consider a standing wave pattern with wavenumber k and frequency ω, forcing a stratified fluid with buoyancy frequency N , for which the background stratification is maintained by Newtonian cooling with damping rate γ. Together, the Newtonian cooling γ, and the viscosity ν, damp the wave amplitude over a characteristic e-folding length Λ = αkc 4 /(νN 3 ), where c = ω/k is the zonal phase speed, and where α = νN 2 /(νN 2 + γc 2 ) is the ratio of viscosity to Newtonian cooling in wave damping. Another essential parameter of the problem is the effective Reynolds number Re = F 0 Λ/ (cν), where F 0 = (u 0 w 0 ) r.m.s. is the wave forcing strength at the bottom boundary. This wave forcing strength further sets a characteristic time scale of low-frequency flow reversals T = cΛ/F 0 [START_REF] Vallis | Atmospheric and oceanic fluid dynamics second edition[END_REF]. Supplemental Material provides details on the simulations as well as estimates of the key parameters of the the quasilinear and nonlinear models, as well as for the Earth's stratosphere.

The parameter range used in the simulations is close to that used in the pioneering work on the subject [START_REF] Plumb | The interaction of two internal waves with the mean flow: Implications for the theory of the quasibiennial oscillation[END_REF][START_REF] Holton | An updated theory for the quasi-biennial cycle of the tropical stratosphere[END_REF], and presented in standard textbooks on geophysical fluid dynamics [START_REF] Vallis | Atmospheric and oceanic fluid dynamics second edition[END_REF]. Notice that the effective Reynolds number is based on an eddy viscosity, meant to represent the turbulent eddy motion at scales smaller than the internal gravity waves and used as a subgridscale parameterization for turbulence in coarse-grained climate models. The actual Reynolds number of the atmosphere based on the kinematic viscosity of the air is higher by many orders of magnitude than the effective Reynolds number. A self-consistent theory for the QBO would require to infer the eddy viscosity from the knowledge of the actual Reynolds number and other problem parameters, but this conundrum has up to now be out of reach. Here we follow a common practice in geophysical fluid dynamics that amounts to: (i) use an eddy viscosity to describe bifurcations occurring under an increase in forcing amplitude, and (ii) test the robustness of these bifurcations in more complex members of the hierarchy of geophysical flow models [START_REF] Dijkstra | Nonlinear climate dynamics[END_REF].

Bifurcation diagrams. To map the bifurcation diagram of the quasilinear 1D model, we performed a large number of simulations spanning effective Reynolds numbers between Re = 2 and 330, covering roughly the relevant range for the Earth's stratosphere (Table 1). For sufficiently low values of Re, the system has only one attractor: a stable point at u = 0. A first bifurcation occurs above the critical value Re c1 ≈ 4.25/(1+α) [START_REF] Yoden | A new look at equatorial quasibiennial oscillation models[END_REF], for which the zonally averaged velocities are attracted towards a limit cycle [START_REF] Plumb | The interaction of two internal waves with the mean flow: Implications for the theory of the quasibiennial oscillation[END_REF][START_REF] Yoden | A new look at equatorial quasibiennial oscillation models[END_REF] corresponding to horizontal mean-flow reversals and downward phase propagation (Fig. 1b). This period-1 cycle arguably reproduces the salient features of the observed QBO before the disruption event of 2016 (Fig. 1a). Figure 1d shows a bifurcation diagram plotted for increasing Reynolds numbers. A second bifurcation from periodic to quasi-periodic regimes occurs above the critical value Re c2 . Additional bifurcations occur at higher Reynolds numbers, with transitions to frequency-locked regimes, and chaotic regimes. The term 'frequency-locking' is often used where a nonlinear oscillator forced at some frequency exhibits, as a dominant response, an oscillation at the forcing frequency. By extension, we use this term here to describe synchronisation between oscillating modes of the dynamical system. As Re increases, new oscillating modes appear in the vertical structure of the mean flow. For example, a unique frequency is observed at all heights for the period-1 limit cycle shown in Fig. 1b, while faster reversals are observed in the lower levels for the frequency-locked regime shown in Fig. 1c. Ultimately, in chaotic regimes, the superposition of these modes yields a fractal-like structure of nested flow reversals (Fig. S1 in Supplemental Material). Such regimes with faster reversals in the lower layers have also been reported in direct numerical simulations driven by a convective boundary layer [START_REF] Couston | Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves[END_REF][START_REF] Showman | Atmospheric Circulation of Brown Dwarfs and Jupiter and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations[END_REF].

The quasiperiodic regime occurring at Re > Re c2 is embedded with a complicated set of frequency-locked regimes (Fig. 1d). The global structure of the bifurcation diagrams is better appreciated by considering, in Fig. 1e, the two-dimensional parameter space spanned by the Reynolds number Re and the parameter α. This figure shows a range of parameters where frequency locked regions are organized into a sequence of staircases, qualitatively similar to Arnold's tongues [START_REF] Arnold | Small denominators. I. Mapping the circle onto itself[END_REF].

Transition to chaos in a similar 1D quasilinear model was reported in Ref. [START_REF] Kim | Gravity wave-driven flows in the solar tachocline[END_REF], which focused only on the purely viscous case, α = 1, with other boundary conditions relevant for the solar tachocline. In fact, this behavior occurs generically in nonlinear systems, with numerous examples in hydrodynamics [START_REF] Swinney | Observations of order and chaos in nonlinear systems[END_REF].

In the case of internal gravity wave streaming, frequency locked states organized into Arnold's tongues were found when the 1D quasilinear model is coupled to an external low frequency forcing mimicking seasonal forcing [START_REF] Rajendran | Synchronisation of the equatorial QBO by the annual cycle in tropical upwelling in a warming climate[END_REF], which is reminiscent of synchronisations phenomena in models of El Nino Southern Oscillations [START_REF] Tziperman | El Niño chaos: Overlapping of resonances between the seasonal cycle and the Pacific oceanatmosphere oscillator[END_REF][START_REF] Jin | El Niño on the devil's staircase: annual subharmonic steps to chaos[END_REF]. Here, by considering a simple monochromatic forcing, and by covering the full parameter space Re -α, we bring to light an unforeseen intrinsic dynamical structure of the underlying quasilinear model.

The 1D quasilinear model is a highly truncated version of the original flow equations. It is thus crucial to see whether the aforementioned bifurcations occur in Navier-Stokes simulations of the fully nonlinear dynamics, including both wave-mean and wave-wave interactions. In Fig. 2a, we performed more than 200 two-dimensional numerical simulations to build a diagram similar to the one obtained with the 1D quasilinear model. These simulations show that the route to chaos is robust to the presence of nonlinear interactions between waves and mean-flow, with transitions from periodic solutions (Fig. 2b) to quasiperiodicity (Fig. 2c), to frequency locking (Fig. 2d) and eventually to chaos. However, significant differences from the quasilinear case are observed in the nonlinear simulations, where bifurcations occur at different effective Reynolds numbers, and where new dynamical regimes emerge. For instance, the large region of period-3 frequency locking obtained in the 1D quasilinear model (Fig. 1c) is replaced by a thin region of period-2 frequency locking for which the symmetry U → -U is broken (Figs. 2a and2d).

Response to external perturbations. By considering a fixed monochromatic wave forcing, we show above that quasi-periodicity arises naturally at steady state in the stratified fluid. This fixed forcing contrasts however with the actual QBO signal, which is driven by timevarying wave forcing and extra-tropical perturbations. In the following, we investigate how the presence of a bifurcation point influences the resilience of a given period-1 Phase space trajectory for Re -1 = 0.009, projected on a 3D space defined by velocities at three different heights: (u1, u2, u3) = (u(z = 0.5Λ, t), u(z = 1.5Λ, t), (z = 3Λ, t)). c. Same as b, but using Re -1 = 0.0068. d. Same as b, but using Re -1 = 0.0058.

QBO-like oscillation to external variability by considering the effect of a time-dependent perturbation superimposed on its reference monochromatic wave forcing.

We first consider the effect of a time-dependent pulse in wave forcing strength, F 0 , mimicking the reported sudden increase in wave activity at the equator in the winter preceding the observed periodicity disruption of 2016 (see Supplemental Material for details on the perturbation). From a dynamical point of view, this perturbation suddenly drives the system out of its limit cycle, until it eventually relaxes back to its original period-1 oscillation over a characteristic time τ . Figures 3a and3b show examples of transient recovery periods for two values of the effective Reynolds number using the nonlinear model. In each case, the time evolution of the mean-flow displays short eastward-flow structures sandwiched between broader westward wind patterns. These higher vertical modes of oscillations, frequently excited in transient disrupted regimes, share qualitative similarities with the periodicity disruption observed in 2016.

Figure 3d shows that the characteristic timescale for recovery diverges as the system approaches the bifurcation point Re c2 . We found similar responses to a pulse in zonal mean momentum, and for the spin-up of the system from a state of rest. The recovery time's divergence is observed both in the quasilinear model and the nonlinear model. The swift increase in recovery timescale observed as the system approaches a bifurcation point is a generic feature of dynamical systems often referred to as "critical slowing down" [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF][START_REF] Kuehn | A mathematical framework for critical transitions: Bifurcations, fast-slow systems and stochastic dynamics[END_REF]. In the climate system context, critical slowing down has proven to be useful in detecting early warnings of a bifurcation point [START_REF] Lenton | Early warning of climate tipping points[END_REF].

Conclusions and perspectives. Our study demonstrates that erratic mean flow reversals are recovered with a simple monochromatic wave-forcing, provided that the forcing strength is sufficiently large. This suggests that similar states previously observed with more complex forcing [START_REF] Couston | Order Out of Chaos: Slowly Reversing Mean Flows Emerge from Turbulently Generated Internal Waves[END_REF][START_REF] Showman | Atmospheric Circulation of Brown Dwarfs and Jupiter and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations[END_REF] can partly be attributed to the intrinsic dynamics of stably stratified fluids, rather than to fluctuations of the forcing itself. The quasiperiodic route to chaos found both in our quasilinear and fully nonlinear simulations reveal that increasing the forcing strength leads to the excitation of fast and shallow bottom-trapped modes nested in deeper and slower vertical modes. These fast bottom-trapped reversals are also excited during the transient response to an external perturbation. Most important, our results have crucial implications for the interpretation of the variability of a QBO-like oscillation: (i) the existence of a second bifurcation is robust in a hierarchy of models (suggesting that it may exist for the atmosphere), (ii) the proximity to this second bifurcation has a strong effect on the response of the oscillation to external perturbations, and consequently (iii) the intrinsic variability of a given oscillation is key to interpret its response to external perturbations.

Several aspects of actual planetary flow such as seasonal forcing, rotation, meridional circulation and twoway coupling between stratospheric and tropospheric dynamics are omitted in the simplified flow models considered in this letter. The interplay between intrinsic modes of variability and these additional features will need to be addressed in future work, but we expect that the the existence of a second bifurcation as well as the critical slowing down approaching this bifurcation will be robust through the whole hierarchy of geophysical flow models with stable stratification. Exploratory 3D simulations with rotation -presented in supplementary materialscomfort our 2D results.

1D quasilinear simulations 2D Nonlinear simulations Stratosphere Re = F0Λ/ (cν)

15

-350 50 -600 2 -400 α = νN 2 / νN 2 + γc 2 0 -1 0.25 0 -0.3 P r = ν/κ ∞ 740 N.A. Fr = c/ (ΛN ) Fr → 0 0.06 0.1 =ω/N ω/N → 0 0.1 10 -4 -10 -3 =ωτu 0 0.1 N.A. =H/Λ 3.5 4.1 1.5 TABLE I. Dimensionless parameters.
For the stratospheric values, we considered F0 = 3 -10 × 10 -3 m 2 s -2 , c = 25 ms -1 , γ = 0.5 -1.5 × 10 -6 s -1 and N = 2.2 × 10 -2 s -1 based on [S3]. We considered Λ ∼ 10 km based on the QBO observations (see fig. 1a in the letter). We consider ν = 0.01 -0.3 m 2 s -1 corresponding to turbulent vertical eddy-diffusivity measured in the lower stratosphere (see e.g. [S2, S5]). The closure used in the 1D quasilinear model reduces the number of dimensionless parameters down to three by assuming ω/N → 0 (hydrostatic approximation) and a low Froude number limit F r → 0 (WKB approximation).

Newtonian cooling are both increased in the 20 upper grid layers.

Poincaré sections. For each combination of parameters (Re, α), experiments are first spun-up over a time t e = 1500T where T = cΛ/F 0 . This time is sufficient for the system to reach its attractor. To combine the information of more than 10 6 simulations into a single bifurcation diagram, we first select two vertical levels: z 1 near the surface and z 2 aloft. Resuming the simulation at statistical equilibrium (t > t e ), we store the values of u(z 2 ) that intersects u(z 1 ) = 0 in the set

O Re,α = { u (z 2 , t) | u (z 1 , t) = 0 } . ( S4
)
The simulations are stopped once 200 values are stored (i.e. after 200 reversals of the lower-level mean-flow u (z 1 )).

For each simulation associated with couples of parameters (Re, α), we build an histogram of the values stored in (S4), using 1000 bins in the range [-c, c]. Histograms corresponding to all values of Re for a fixed α = 0.6 are drawn horizontally in figure 1c using a binary colour-map.

To collapse the information of the Poincaré sections into a 2D bifurcation diagram (α, Re -1 ), we compute the ratio of populated bins to the total number of reversals for each histogram in the set (S4). This ratio with values in [0, 1] provides an empirical estimate of each histogram's distribution and allows for an extensive classification of the different dynamical regimes (see Fig. 1d in the letter)

Recovery from a perturbation. We consider perturbations to a given period-1 QBO-like oscillation. Three types of external perturbations are considered at t = t p : (i) a pulse in wave amplitude (representing a sudden increase of the underlying tropical convection) (ii) a body force acting directly on the mean flow (representing a reorganization of the mean flow due to extratropical perturbations) (iii) a reboot of the oscillations from a state of rest (the recovery time is then equivalent to the spin-up time).

Perturbation (i) is modeled using a time-dependent momentum flux amplitude F 0 in Eq. ( S1) for the quasilinear model and in Eq. ( S3) for the nonlinear simulations:

F 0 (t) = F 0,p 1 + 9e -1 2 10 t-tp T qbo 10 , (S5)
where F 0,p is a constant forcing amplitude corresponding to a periodic regime with period T qbo . Perturbation (ii) is represented in the quasilinear model by an additional body forcing term F bulk in the r.h.s. of the mean flow equation (1):

F bulk (z, t) = 10 T qbo e -1 2 10 t-tp T qbo 10 e -800 z-zp zmax 2 , (S6)
where z p = 0.2z max sets the height. Results are insensitive to the specific choices of z p . For all three types of external perturbations, the system is driven away from its steady state period-1 limit cycle and then freely recovers back to the cycle. To estimate the recovery timescale, we first introduce the running mean-square transition from periodic to quasiperiodic oscillations. Results of Fig. S4 show that the value of Re c2 is converging for a resolution δ z = 3.5/60, corresponding to the reference resolution used in this work.

u 2 (t) = 1 T qbo H t+T qbo /2 t-T qbo /2 H 0 u 2 (z, t ) dzdt , (S7)
As far as critical slowing down is concerned, the response of the system approaching the bifurcation from periodic to quasi-periodic state will be robust to higher resolutions, as the threshold Re c2 and the nature of the bifurcation remains the same. However, the details of the response, including the vertical structure of transient oscillations and the prefactor of the recovery time power law may be affected by a change in resolution.

3D NONLINEAR SIMULATIONS WITH ROTATION

In this section, we solve the 3D Navier-Stokes equations with rotation approximated by an equatorial beta-plane. The horizontal momentum equation in (S2) now writes

∂ t u + u.∇u -βyv = -∂ x φ + νδ 2 u -γ u (u -u 0 ) ∂ t v + u.∇v + βyu = -∂ y φ + νδ 2 v , ( S9
)
where the velocity vector field is now 3D with u = u î + v ĵ + w k. β is the Rossby parameter. Let us denote L y the length of the added meridional dimension and δ y the associated resolution. We consider an horizontal aspect ratio L y /L x = 1 and resolution ratio δ y /δ x = 1, with free-slip lateral boundary condition at y = ±L y /2. We explore a weak rotation case, for which the equatorial Radius of deformation L d = N Λ/β is much larger than the meridional extension of domain: L d /L y = 96. All other parameters are identical to the 2D nonlinear simulations, including the forcing, constant along the y direction. The chosen initial condition breaks the meridional invariance. Fig. S5 shows that bifurcation from a periodic regime (panel a) to a quasiperiodic regime (panel b) occurs when the Reynolds number is increased from Re = 125 to Re = 250. This demonstrates that the intrinsic variability observed in the 1D quasilinear model is robust to the presence of 3D wave-wave interactions.

In the Earth's stratosphere, the parameter L d /L y is smaller than one, and equatorial waves are trapped along the equator over a typical scale of the order of the deformation radius. Future work is needed to explore the robustness of the bifurcation point Re c2 in the presence of strong rotation.

FIG. 1 .

 1 FIG. 1. Bifurcations in the 1D quasilinear model. a. Observations of the Earth's atmospheric quasi-biennial oscillations. Hovmöller diagram of monthly averaged zonal winds measured by radiosonde in the lower stratosphere above Singapore (1.4 • N) [27]. b. Hovmöller diagram of the mean-flow u(z, t) in stationary regime for α = 0.6 and Re -1 = 0.06, using the 1D quasilinear model. c. Same as b, but using Re -1 = 0.025. See Supplemental Material for similar time-height plots in other regimes. d. Bifurcation diagram for α = 0.6, obtained for each value of Re -1 by considering the value of u at two different heights z1 and z2, and then by plotting u(z2) at times when u ( z1). The dashed red line corresponds to the first bifurcation (from rest to period-1) occurring at Rec1. The dashed blue line corresponds to the second bifurcation (from period-1 to quasiperiodic) occurring at Rec2. e. Bifurcation diagram in parameter space α, Re -1 . The coloured field is an empirical estimate of the area of the attractor projected in panel d. Low values (in white) correspond to QBO-like regions (period-1) and frequency locked regions (rational numbers). The vertical black line at α = 0.6 corresponds to the bifurcation diagram plotted in panel c.

FIG. 2 .

 2 FIG. 2. Bifurcations in Navier-Stokes simulations. a. As in Fig 1d, but showing the bifurcation diagram obtained with the nonlinear simulations, using α = 0.6. Selected values of Re (marked in orange, red and purple) correspond to panels b, c, and d. The dashed blue line corresponds to the second bifurcation (from period-1 to quasiperiodic) occuring at Rec2. b.Phase space trajectory for Re -1 = 0.009, projected on a 3D space defined by velocities at three different heights: (u1, u2, u3) = (u(z = 0.5Λ, t), u(z = 1.5Λ, t), (z = 3Λ, t)). c. Same as b, but using Re -1 = 0.0068. d. Same as b, but using Re -1 = 0.0058.

FIG. S2 .

 S2 FIG. S2. Bifurcations in the 1D quasilinear model. a. A Poincaré section is shown for varying values of Re -1 and α = 0.6 (see Methods). b. Projection of phase-space trajectory in a 3D space (u1, u2, u3) = (u(z = 0.1Λ), u(z = 1.5Λ, t), u(z = 3Λ, t)) for Re -1 = 0.059. c. Same for Re -1 = 0.045. d. Same for Re -1 = 0.025. e. Same for Re -1 = 250. f. Hovmöller diagram of the mean-flow u (z, t) for Re -1 = 0.059. Time is rescaled by T = cΛ/F0. The velocity u ranges from to -c (blue) to +c (red). The horizontal dotted lines highlight the height z = 0.1Λ, z = 1.5Λ and z = 3Λ, associated with the 3D projections plotted in panels b to e. g. Same for Re -1 = 0.045. h. Same for Re -1 = 0.025. i. Same for Re -1 = 0.0004.

FIG. S3 .

 S3 FIG. S3. Additional bifurcation diagrams of the 1D quasilinear model. a. Poincaré sections for each value of Re -1 using α = 0 (see Methods). b. Same for α = 1/3. c. Same for α = 2/3. d. Same for α = 1.

FIG. S4 .

 S4 FIG. S4. Resolution dependence in the 1D quasilinear model a Poincaré sections for each value of Re -1 with α = 0.6. The spatial resolution used is δz = 3.5/15. b. Same for δz = 3.5/30. c. Same for δz = 3.5/60. d. Same for δz = 3.5/120. e. Same for δz = 3.5/240.

* antoine.venaille@ens-lyon.fr

Supplemental Material for 'Periodicity disruption of a model quasi-biennial oscillation' Antoine Renaud 1 , Louis-Philippe Nadeau 2 , Antoine Venaille 1 1 Univ Lyon, Ens de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France 2 Institut des Sciences de la Mer de Rimouski, Université du Québec à Rimouski. Rimouski, Québec, Canada (Dated: January 15, 2019)

METHODS AND DETAILS ON NUMERICAL SIMULATIONS 1D quasilinear simulations. Using a static Wentzel-Kramers-Brillouin (WKB) approximation to compute the wave field for a given mean-flow [S4], the wave-induced Reynolds stress u w is parametrised by the formula given in Eq. (S1).

This formula is derived under the hydrostatic balance assumption (valid in the limit k|c ± u|/N → 0) and the weak damping assumption (valid in the limit γ/(k|c ± u|)

1 and νN 2 /(k|c ± u| 3 ) 1). Assuming that the characteristic vertical length for u is Λ, then the small parameter needed in the WKB approach is the Froude number F r = c/(ΛN ) → 0. In practice, the different assumptions are most certainly violated. However, this set of equations has long been recognized as a useful model to probe the salient features of QBO reversals.

We solve numerically using a centered second-order finite difference method with grid size δ z = H/60, and a second-order Adams-Bashforth scheme with time-step δ t = 0.005T ; T = cΛ/F 0 . A no-slip condition is used at the bottom boundary, z = 0, and a free-slip condition is used at the upper boundary, z = H. Singularities in Eq. (S1) appear when u = c (critical layers). These singularities are treated as follows: at a given height z = z c , if the absolute value of u reaches locally a value higher than c, then the corresponding exponential is set to zero for all z ≥ z c . The definition and value of each of the model's dimensionless numbers are given in table I.

2D nonlinear simulations. The fully nonlinear simulations are conducted using the MIT general circulation model [S1] solving the 2D Navier-Stokes equations under the Boussinesq and hydrostatic approximations

where u = uî + w k is the velocity field, u h is its projection on the horizontal plane (x, y); b = g (ρ 0 -ρ) /ρ 0 is the buoyancy; ρ is the density and ρ 0 is a reference density; g is the gravitational acceleration; φ = P/ρ 0 + gz; P is the pressure; ν is the viscosity coefficient; κ is the buoyancy diffusion coefficient; γ u and γ are the rates at which the momentum and buoyancy are linearly restored to the reference profiles u 0 and b 0 , respectively. The domain is a Cartesian grid, periodic in the zonal direction, with zonal length L x = 2π/k and height H . The horizontal and vertical resolutions are respectively δ x = L/26 and δ z = H/200. A free-slip condition is used at the bottom boundary, while a free-surface condition is used at the top.

The zonal momentum equation is forced at the bottom boundary using a linear velocity relaxation γ u = δ b /τ u , where τ u is a relaxation timescale, and δ b is a delta function equal to 1 for the bottom grid-point and 0 for all other vertical levels. In this last grid-point, velocity is relaxed to a zonally periodic standing wave pattern

where F 0 controls the wave momentum flux amplitude at the bottom. This forcing is thought to generate a standing internal gravity wave field while enforcing an effective no-slip condition for the mean flow u. Buoyancy is relaxed to the linear profile b 0 = N 2 z. To avoid any wave reflection at the upper free-surface, the vertical grid spacing and the where T qbo is the period of the limit cycle. At steady state equilibrium, this running mean-square has a constant value u 2 ∞ . Assuming a pulse shorter than the period of the limit cycle (see Fig. 3c in the letter), occurring at time t p , the recovery timescale is then defined by

We reproduced figure 3d of the letter in logarithmic scale in order to exhibit the power-law like scaling of the recovery timescale as the system approaches Re c2 . It proved very difficult to deduce a precise value for the critical exponent as the uncertainty on the value of Re c2 echoes on it. However, the critical exponent remains close to -1.

VERTICAL FLOW STRUCTURE IN DIFFERENT REGIMES AND EFFECT OF RESOLUTION

Bifurcations in the quasilinear model. In order to develop intuition for the underlying dynamics of the bifurcation diagrams of Fig. 1d, we show in Fig. S2 phase space trajectories (panels b-e) and hovmöller diagrams of the mean-flow (panels f-i) for four selected values of the Reynolds number. Shown are examples for a period-1 limit cycle (panels b and f), a quasiperiodic oscillation (panels c and g), a frequency locked oscillation with frequency ratio 1/3 (panels d and h), and a chaotic oscillation (panels e and i).

Additional bifurcation diagrams obtained for different values of α are shown in Fig. S3. Although sharing a common qualitative structure, each bifurcation diagrams show distinct interesting features. For example, Fig. S3b shows that the upper quasiperiodic region vanishes almost entirely when α approaches 1/3. Fig. S3d (α = 1) shows a spontaneous breaking of the symmetry U ↔ -U occurring in one of the frequency locked states (Re -1 ∼ 0.045), while all the frequency-locked regimes preserve this symmetry at lower values of α in panels a, b and c.

Effect of the resolution in the quasilinear model. In order to test the robustness of the quasilinear model results to resolution, we show in Fig. S4 five bifurcation diagrams for which the vertical resolution has been successively doubled from δ z = 3.5/15 to 3.5/240. Panel c corresponds to the reference resolution used in Fig. 1d. Results show a strong dependence on vertical resolution, in particular for the structure of the embedded frequency locked regimes. However, the essential feature relevant to the periodicity disruption is the second bifurcation point Re c2 , marking the


O Bokhove
email: o.bokhove@leeds.ac.uk

M Kelmanson

T Kent
email: t.kent@leeds.ac.uk

G Piton

J.-M Tacnet

Using 'flood-excess volume' to assess and communicate flood-mitigation schemes

Flood-excess volume: definition and use

Flood-Excess Volume (FEV): the volume V e of water in a flood-peak at a river gauge location due to river levels exceeding a relevant threshold h T , chosen such that, for any river level h above h T , flooding occurs. Best approximation:

V e ≈ N m k=1 Q( hk ) -Q T ) ∆t. • h: water level [m] • Q: discharge [m 3 /s] • Q T = Q(h T ): threshold discharge • T f = N m ∆t: flood duration
• GOAL: to quantify and communicate the efficacy of various floodmitigation measures in a straightforward and concise manner. • IDEA: calculate the FEV for a flood event of interest and express as the capacity of a 2m-deep square 'flood-excess lake' with side-lengths O(1km). • OUTCOME: a graphical tool that (i) contextualises the magnitude of the flood relative to the river valley and (ii) facilitates quick and direct assessment of the contribution and value of various mitigation measures. • express the FEV as a 2m-deep 'flood-excess lake' of sidelength 2.15km: FEV≈ 9.34Mm 3 ≈ (2150 2 × 2)m 3 .

• given a calculation (or estimate) of potential flood storage volume and associated cost of each mitigation measure, the 'flood-excess lake' can be partitioned accordingly and overlaid with a cost per 1% of FEV mitigated Hypothetical scenario for Leeds' flood alleviation scheme II (FASII+), comprising flood walls and flood-storage sites:

2b. Case study II: River Calder, Boxing Day 2015

Exploratory flood-alleviation scheme comprising (i) temporary storage in reservoirs, (ii) upscaled 'leaky' debris dams, and (iii) tree planting:

• takes into account uncertainty in storage capacity;

• draw-down and control of reservoirs has great potential;

• major upscaling of leaky dams can have a significant and cost-effective impact; • mean FEV mitigated is 50%: more measures (e.g., flood walls) required to offer full protection.

Right: woody-debris 'leaky' dams, an example of Natural Flood Management (NFM), are being installed in tributaries of the River Calder. Photo courtesy Robin Gray (Pennine Prospects).

2b. Case study III: River Brague, October 2015

An alternative to raising flood walls: giving room to the river (GRR) increases the river width in order to increase the discharge capacity for a similar water depth. • New 'GRR' rating curve reduces FEV for same h T (left)

• GRR is assessed favourably alongside flood walls and retention measures in an exploratory flood-mitigation analysis for the River Brague (right)

Highlights

• a complementary way of classifying flood events, to be used either prior to or in tandem with more detailed hydrodynamic numerical modelling • a new protocol to optimise assessment of mitigation scenarios, including cost-effectiveness analyses, for the benefit of policy makers • encourages evidence-based decision-making for flood-mitigation schemes

2a.

  Case study I: River Aire, Boxing Day 2015 Q: what fraction of the FEV is reduced, and at what cost, by a suite of mitigation measures?

  -room-to-the-river transverse profile b 1

We thank the Environment Agency for providing data and reports for the UK cases, and SIAQUEBA and DDTM06 for the Brague case. Work funded by EPSRC [grant nos. EP/M008525/1 and EP/P002331/1] and the EU H2020 project NAIAD [grant no. 730497].


N Jamet

R Prud'homme

L Gottesdiener

MODELING OF DISCONTINUITIES IN TWO-PHASE COMPRESSIBLE FLOWS 1

Keywords: multiphase flow, interface, suspension, shock wave, asymptotic method

The weak shock wave in a dilute suspension and the downstream relaxation zone are considered as examples of interfaces in multiphase flows. Thermodynamics of irreversible processes is used to analyze momentum and energy transfer between phases. The singularity near the sonic regime is studied by means of asymptotic expansions. Interfacial constitutive relations are then written for interfaces which obey to linear irreversible thermodynamics.

INTRODUCTION

Transport phenomena in two-phase flows have been studied since many years [START_REF] Fortier | Mécanique des suspensions[END_REF][START_REF] Marble | Dynamics of dusty gases[END_REF], Nigmatulin 1989[START_REF] Soo | Multiphase fluid dynamics[END_REF]) especially in gas solid suspensions, but also in liquid vapor non reacting and reacting flows [START_REF] Kuentzmann | Aérothermochimie des suspensions Mémoire des Sciences Physiques[END_REF][START_REF] Barrère | Equations fondamentales de l'aérothermochimie[END_REF][START_REF] Delhaye | Jump conditions and entropy sources in two-phase systems[END_REF][START_REF] Williams | Combustion Theory[END_REF]. Interfacial transfers through discontinuities of these systems are not so well known and most of the time empirical expressions are not available due to small particle size and numerical calculations are difficult to perform because of the various (hydrodynamic, molecular and particle) scales of the problem. Such discontinuities appear in emulsions and foams as well as in sprays. Generally, phase transitions between pure liquid or liquid suspensions are considered as interfaces. But shock waves in pure gases or in sprays and dusty gases, occurring for instance in motors of supersonic aircrafts or in rocket nozzles, can also be modeled as interfacial media. Such generalized interfaces have their own constitutive relations. We will illustrate this with the example of shock waves and relaxation zones in dusty gases of which basic equations and spatial structure are well known especially in stationary case. To do this we will use the same way as for premixed flames [START_REF] Prud'homme | Modélisation interfaciale des flammes minces[END_REF], for some boundary layers [START_REF] Prud | Fluides hétérogènes et réactifs: écoulements et transferts[END_REF] as well as for other generalized interfaces [START_REF] Gatignol | Mechanical and thermodynamic modelling of fluid interfaces[END_REF]. The interfacial balance equation for a given quantity of which the value per unit mass is  in the bulk and S  at the interface can be written in the form [START_REF] Gatignol | Mechanical and thermodynamic modelling of fluid interfaces[END_REF])

S a // a S S S S S a S a W J ] )] W V ( J [ [ W dt ) ( d S                              (1)
1 Unpublished paper, which has been written in 1999 at Laboratoire de Modélisation en Mécanique. 2 PhD at the University Pierre et Marie Curie, defended in 1998. 3 Presently Emeritus Research Director -E-mail : roger.prud_homme@upmc.fr 4 Assistant Professor presently retired.

2

where  and a  are respectively the bulk and interfacial densities, V  is the material velocity in the bulk, dt d S W 

is the material derivative associated to an interfacial velocity S W  of which the tangential component to the interface is equal to the material tangential interfacial velocity

// S V 
and the component normal to the interface is the normal displacement velocity of the interface   w ,   being the unit normal to the interface orientated from side I to side II (Figure 1), J  and a J  are respectively the bulk and interfacial fluxes and S  is the interfacial production per unit mass,

I II f f ] ] f [ [  
is the jump of f through the interface and

                           1 S is the tangential divergence operator.   S W    w S V  // S V  S I II Figure 1 Interfacial velocities
Problems of wave propagation in two-phase media have been revisited since a few years. Some authors have studied gas/solid suspensions [START_REF] Sharma | Weak discontinuities in a non equilibrium flow of a dusty gas[END_REF], Margulies & Schwartz 1994) whereas others have focused on liquid/vapors systems that can be found in many industrial applications [START_REF] Jackson | An equation set for non-equilibrium two-phase flow and an analysis of some aspects of shocking, acoustic propagation, and loss in low pressure wet steam[END_REF], Young & Guha 1991).

We consider here the case of the propagation of a normal stationary shock wave in a two-phase medium (gas/solid suspension). The study will be made in the frame linked to the wave [START_REF] Jamet | Modélisation des discontinuités dans les écoulements diphasiques compressibles[END_REF]. In this case, equation (1) reduces to

S a S ] )] W V ( J [ [             
(2)

As volumetric fraction of particles is negligible compared to fluid volumetric fraction and size of particles is much larger than width of gaseous shock wave, it can be deduced that the internal structure of the wave is not altered by the particle cloud.

The two-phase shock wave is split in a gaseous shock wave and a relaxation zone between phases [START_REF] Rudinger | Some properties of shock relaxation in gas flows carrying small particles[END_REF].

The proposed method refers to Thermodynamics of Irreversible Processes (TIP) to express when possible, momentum and energy transfers by mean of phenomenological coefficients [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]. When dissipation is due to both viscous effects and heat conduction, the generalized forces are A v and A T and the corresponding generalized fluxes are J v and J T , then the entropy production is of the form

T T J A J A    v v
(3) and when linearized irreversible thermodynamics is valid, we will obtain in each case relations of type

J v = L v A v , J T = L T A T (4)
In all cases we need to obtain the expression of the interfacial entropy production  a .

Two ways are possible: by writing directly the entropy balance law or by integrating the bulk entropy production  through the interfacial layer. For the last one we must suppose that the interfacial layer can be considered, at small scale, as a continuous medium in the three dimensions of space. This is not always possible for shock waves.

For sake of simplicity, we consider a two-phase medium with rigid and spherical solid particles of same size of about one micron with no interaction between them and for the fluid phase a perfect compressible gas like air. The two phases are considered as a continuum (size of particles much larger than mean free path of molecules in gas). Applying the classical method used for interfacial zones [START_REF] Ishii | Thermo-fluid dynamic theory of two phase flow[END_REF][START_REF] Barrère | Equations fondamentales de l'aérothermochimie[END_REF][START_REF] Gatignol | Mechanical and thermodynamic modelling of fluid interfaces[END_REF] to this suspension, one gives interfacial equations for each phase and for the whole medium. In addition to interfacial equations of mass, momentum and energy, one writes entropy production rate at this interface.

Because the first part of the two-phase shock wave is a gaseous shock, we will study first of all the case of a stationary gaseous shock wave in a one-dimensional ideal gas flow.

Next, the two-phase shock wave in a dilute suspension will be considered. The jump relations between supersonic and subsonic flows are easily obtained. For very weak shock waves, the mach number M 0 is unity. We are near equilibrium, nevertheless the theory of linear non equilibrium thermodynamics can not be applied. On the contrary, it will be shown that the two-phase relaxation interface downstream the shock agrees with linear TIP, at the condition of deducing accurately the generalized forces and fluxes from the entropy production. This transonic situation presents a singularity when M 0 2 -1 is small just as of the mass ratio K of dispersed and fluid phases. An asymptotic method is applied to analyze this zone where there is some difficulties to obtain separately the thermal and viscous interfacial dissipations.

For relatively weak shock waves with higher Mach number, but not too far from equilibrium, the problem becomes non-singular and it is then possible to apply directly classical forms of interfacial equations and to use the linearized theory of TIP to close the system and express transfer phenomena with phenomenological coefficients that rely thermodynamic fluxes and forces linearly. Here the thermal and viscous effects appear separately.

It can be remarked that several scales of work must be considered to analyze interfaces in multiphase flows (Jamet, Prud'homme & Gottesdiener 1997): First, the molecular scale is used to determine constitutive relations for dispersed phase and carrying fluid; then equation of state and constitutive relations are established. These laws come from an integration at a larger scale. Finally it gets coefficients of viscosity, diffusion and conduction. The second scale or mesoscopic scale deals with the particle and its fluid environment. This former is then a macroscopic object flowing in a fluid and macroscopic Navier-Stokes equations are used to describe its evolution. With simplifying assumptions like Stokes one, drag law and heat transfer equation for the particle can be deduced. The third scale, called macroscopic 1, is the scale where two-phase medium is a quasi-continuum. When writing balance laws and especially entropy law, phenomenological relations are established and coefficients can be deduced of results obtained at mesoscopic scale. Finally the fourth scale, or macroscopic 2, where a thin slice of two-phase medium is assimilated to a surface. At this scale we also write balance laws that lead to phenomenological equations with coefficients that are similar to exchange coefficients. The results obtained at the third scale lead to the determination of these coefficients. In this paper, we will focus on the fourth scale and then we must first work at the third scale to finally express interface transfer coefficients.

THE GASEOUS SHOCK WAVE AS AN INTERFACE

Navier-Stokes equations applied to a stationary one-dimensional flow

Let us consider the one-dimensional flow of a simple gas. The balance equations of this flow are particular cases of the general balance equations, and the entropy production is deduced. The two generalized forces and the two conjugate generalized fluxes, which appear in the resulting entropy production formula, correspond to heat conduction and viscous dissipative phenomena. Closure relations can be derived from linearized TIP (de Groot 1963), and the following system is obtained

                          T T v v xx xx 2 2 xx 2 x T 1 q x u T T ) entropy ( 0 x T q x s u t s ) energy ( 0 x T u x u p x q x 2 u e u t 2 u e ) momentum ( 0 x T x p x u t u ) mass ( 0 x u t J A J A                                                                (5)
where u is the velocity, p the pressure, T the temperature, e and s the internal energy and entropy per unit mass, T xx the viscous stress, q the heat flux, and  the entropy production. One has, for the generalized forces

  x u , x T 1 v T       A A (6)
and, for the generalized fluxes

T , q xx T v T   J J (7)
The linearized relations of TIP (Thermodynamics of Irreversible Processes) are then

v v T 2 T T 3 4 , T k A J A J    (8)
where k and  are the heat conduction and viscosity coefficients. Introducing these constitutive relations in the previous system (5), one obtains the one-dimensional Navier-Stokes equations.

Applied to the stationary case, these equations give

                0 x d dx du u d 3 4 x d u p d x d x d T d k d x d 2 u e u d 0 x d x d u d d 3 4 x d p d x d u d , 0 x d u d 2 2               (9)
The entropy equation becomes

  2 2 2 x d u d T 3 4 x d T d T k x d T d T k x d d x d s u d                                 (10)
A first integration of (9) leads to

E x d u d u 3 4 u p x d T d k 2 u e m I x d u d 3 4 p u m , m u 2                      (11)
where m  , I and E are constants. For given values of these constants, we obtain generally two non-dissipative solutions (du/dx=dT/dx=0), that we call (0) and (I). One of the solutions is subsonic and the other supersonic. Both solutions become identical when the constants are chosen in such a manner that the resulting flow is sonic.

Discontinuity in a stationary one-dimensional flow

In the stationary case, when the upstream flow is supersonic, the associated downstream flow can be subsonic and the two flows are separated by a discontinuity, which is a stationary shock wave (Figure 2). The unit mass flow rate m  , the momentum I and the total energy E keep their values across the discontinuity. The balance equations of the shock are obtained (balance equations are also obtained for discontinuities of shock type by several authors [START_REF] Germain | Mécanique des milieux continus[END_REF][START_REF] Germain | Mécanique des milieux continus[END_REF][START_REF] Jaumotte | Chocs et ondes de choc[END_REF][START_REF] Barrère | Equations fondamentales de l'aérothermochimie[END_REF])

) energy total ( 0 ] ] 2 u h [ [ m ), momentum ( 0 ] ] p [ [ ] ] u [ [ m ), mass ( 0 ] ] m [ [ ] ] u [ [ 2           (12)
These equations have the usual form of interfacial equations (2). From these conservation equations, one deduce the entropy balance equation and the Clausius-

Duhem inequality 0 ] ] s [ [ m a     (13)
For an ideal gas, the entropy production becomes

     v p ln c m v a  (14)
where v ~= 0 / I , p =p I /p 0 , =c p /c v =Cte. The result (14) has been deduced from (13) using the expression of the entropy for an ideal gas s= c v Log p +c p Log v +Cte. The solution of the stationary shock wave equations is given, for the ideal gas, by the Hugoniot theory [START_REF] Jaumotte | Chocs et ondes de choc[END_REF]. It is easy to find

Interfacial

2 2 0 2 0 2 0 2 0 2 0 ) 1 ( ) 1 M )( 1 M 1 )( 1 ( 2 1 v p T ~1 ) 1 M 1 ( M 2 1 p ~) 1 M 1 ( 1 2 1 u ṽ ~                      (15) with M 0 =u 0 /a 0 , a 0 =(rT 0 ), r=R/M, 0 I 0 I T T T , u u u ~ 
, R the universal gas constant and M the molar mass of the gas. For weak shock waves (near-equilibrium shocks) one put

M 0 2 -1=<<1 (16)
Then, the entropy production becomes in first approximation

    3 2 v a 1 3 1 2 c m          (17)
This last result is well known: the entropy production through a weak shock wave is proportional to the cube of (M 0 2 -1).

Remarks

-The foregoing equations are interfacial equations and have the form already given in Introduction. Here the material interface is at rest (the chosen frame is attached to the shock and both upstream and downstream flows are uniform).

-None assumption was made on the validity of the N.S. equations inside the shock (the interfacial layer).

-Thermal conduction and viscous effects do not appear separately in the entropy source formulas ( 14) or ( 17). If it was the case, this would lead to a separation between the thermal force and the momentum flux terms

] ] [ [ and ] ] [ [ u T
. Indeed, we obtain a global expression of  a which is, near equilibrium, in

 3 =(M 0 2 -1) 3 , when ] ] u [ [ and ] ] T [ [ are in . -If we take the generalized flux A S = ] ] u [ [ , we have A S = ] ] u [ [ =-((+1)/2)  (18)
The conjugate generalized force is then

        2 S 5 v 2 3 v S a S 1 3 1 16 c m 1 1 8 c m A A J                    ( 19
)
A vS corresponds to a stress T xxS . One sees that the generalized force is proportional to the square of the generalized flux. This is a non-linear irreversible phenomena and, consequently, the obtained interfacial constitutive relation is out of the field of linearized TIP.

Structure of a shock wave

The internal structure of the shock can be studied assuming that the Navier-Stokes equations are valid inside the wave. In the frame of the kinetic theory of gases [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF][START_REF] Chapman | The mathematical theory of non uniform gases[END_REF][START_REF] Smolderen | Structure des chocs et théorie cinétique des gaz[END_REF], it is shown that the Navier-Stokes equations are obtained and are valid for the small Knudsen numbers (the Knudsen number is the ratio of the free mean path by the macroscopic length). It is also established that the thickness of a shock wave is of order of magnitude of the mean free path. But the Navier-Stokes equations are valid to describe the shock structure only of weak shock waves. This last case has been studied by [START_REF] Smolderen | Structure des chocs et théorie cinétique des gaz[END_REF] who has made approximations and obtained analytical formulas. The system ( 11) is solved and this gives an infinity of solutions. For a given Mach number M 0, only one of these solutions is compatible with the downstream solution (15). The found solutions are put in the plane (u,T) in Figure 3. By integration across the interfacial layer 0I, the entropy production ( 9) gives again the value of  a found in ( 13) and ( 14). If we need to obtain separately the thermoconductive and viscous effects we have to calculate, using ( 6)

            I 0 xx av I 0 aT du T T and T 1 d q (20)
This calculation must give again

 a = aT + av (21)
In the case of linearized TIP, we deduce from ( 10)

       I 0 av I 0 2 aT du x d u d T 3 4 , dT dx dT T k ( 22
)
We remark that, even though the value of  a do not depend of the internal processes inside the shock wave, each part  aT and  av of  a will depend of the bulk constitutive relations in the interfacial layer.

With the approximation of [START_REF] Smolderen | Structure des chocs et théorie cinétique des gaz[END_REF], we have: h+u 2 /2E=Cte. Then dT/duu/c p and

2 2 p 2 v T M Pr' 1 c T u 4 k 3       
, where Pr'=4  c p /3 k.

In the case of weak shock waves, the Mach number is near unity and it can be written Pr'

1 v T av aT         (23)
Notice that the expressions of  av and  aT can be directly obtained from the following equation that [START_REF] Smolderen | Structure des chocs et théorie cinétique des gaz[END_REF] deduced from the above approximation

    I 0 u u u u 1 Pr' 1 1 2 1 x d u d m u 3 4           (24)
Indeed with this approximation

                 I 0 I 0 u u I 0 u u av u d T u u u u u 1 Pr' 1 1 2 m u d x d u d T 3 4 
and for uTu 0 T 0 , in the transonic case

      3 2 v av 3 I 0 av 1 3 1 2 c m 1 Pr' 1 1 u u 1 Pr' 1 1 12 m                      (25)
For heat transport we have from ( 22

                   (26)
The founded values of  av and  aT satisfy (21) if  a is given by ( 17).

SHOCK WAVE IN A DUSTY GAS

Structure of a two phase shock wave

Let us consider now a weak stationary shock wave in a two-phase medium. Upstream (0) and downstream (II) of this shock wave, which is assumed to be planar and normal to the flow, the velocity is uniform, and equilibrium conditions are achieved between the gas and condensed particles (which are assumed to be solid small spheres). A hypothesis is often used to simplify the situation: one assumes that the real two-phase shock is a gaseous shock followed by a gas/particle relaxation zone (Figure 4). This assumption has been well verified for dusty gases by [START_REF] Marble | Dynamics of dusty gases[END_REF]. Just after the gaseous shock wave 0I, the velocity and temperature of the dispersed phase remains unchanged. At the end of the relaxation zone, one finds again the equilibrium downstream values of the two-phase shock wave 0II.

A two-phase medium at equilibrium is analogous to a simple gas, with adequate calorimetric coefficients. Then, it is easy to determine the downstream values of the two-phase shock wave. These results will be summarised in the following.

The entropy production a  through the two-phase shock wave 0II, is the sum of the entropy production  a of the gaseous shock wave 0I and the entropy production  a of the relaxation layer III. Then, the entropy production  a can be deduced by simple difference between the entropy productions of the two shocks. The sound velocities are distinct for the gas, which is assumed to be ideal, and for the suspension [START_REF] Kuentzmann | Aérothermochimie des suspensions Mémoire des Sciences Physiques[END_REF], of which the dispersed phase is assumed to be incompressible (it is a solid phase), with a specific heat equal to the constant pressure specific heat of the gas phase. However, for dilute suspensions, these two sound velocities are close to one another. The same situation occurs for chemically relaxing flows [START_REF] Prud'homme | Sur les écoulements transsoniques avec relaxation chimique ONERA[END_REF][START_REF] Prud'homme | Study of transsonic flow with chemical reactions by the small perturbation method. The one and two-dimensional problem[END_REF].

The corresponding Mach numbers are different but not by much. It is well-known that this problem is singular in the transonic zone, where the two Mach numbers are near unity, although it can be solved with the help of the asymptotic expansion theory.

The calculations are summarised in the following.

Simplified bulk balance equations for a dilute suspension

In the present problem, we assume that heat conduction and momentum transfer are limited to exchanges between gas and particles and we consider a simplified model of dilute suspension [START_REF] Marble | Dynamics of dusty gases[END_REF][START_REF] Kuentzmann | Aérothermochimie des suspensions Mémoire des Sciences Physiques[END_REF][START_REF] Prud | Fluides hétérogènes et réactifs: écoulements et transferts[END_REF]. The volume occupied by the dispersed phase can be neglected and, for a onedimensional continuous flow, one has (see section 4 for detailed considerations)

                                                                              ) energy total ( 0 x u p x 2 u e u 2 u
                               ) phase dispersed
                                                    T F , Q , u u , T 1 T 1 u u T F T 1 T 1 Q 0 x s u s u t s s v T p v p T v v T T p p p p p p p J J A A J A J A (29)
In these equations, ( p )designates the particles of the dispersed phase,  p and  are the mass of particles per unit volume of the suspension and the gas density, u p and u are the velocities of the two phases. The only pressure p is that of the gas, e p and e are the specific energy per unit mass of each phase, F is the force exerted by the gas on the dispersed phase per unit volume and Q is the heat transferred from the gas to the particles per unit volume of the suspension. The linearized relations of Irreversible Thermodynamics may be written

             T p p p v p p p l T T u Q l u u u F ( 30
)
where l v and l T are reference lengths for momentum and heat exchanges. These reference lengths are connected to the phenomenological coefficients of formulas (4) by

v p p v T p p p T L T / u l , L / T T u l     .
Entropy production in (29) has the form (3) and here A T and A v result of an integration of three-dimensional quantities through thermal and viscous boundary layers surrounding each particle. The reference lengths can be determined by a study at mesoscopic scale. In Stokes regime one has [START_REF] Kuentzmann | Aérothermochimie des suspensions Mémoire des Sciences Physiques[END_REF])

k 3 a c l , 9 a 2 l s ps 2 T ps 2 v       
, where  is the particle radius (the particles are assumed spherical),  ps the specific density of the particles, s c their mean specific heat,  , k the mean heat conductivity and viscosity of the carrier gas and a the sound velocity of the gas.

For a stationary flow, mass conservation gives

Cte m u , Cte m u p p p        
, and we put

Cte K m m u u p p p        (31)
The other equations ( 27) -(30) become, in the frame of linearized TIP As it is written above, we will assume that the gas phase is an ideal gas and that the dispersed phase is made of solid non deformable particles. The energy per unit mass of particles is: e p =c s T p , and specific heat c s is equal to the constant pressure specific heat c p of the gas.

                              

Dusty gas at equilibrium

At equilibrium, we have T p =T and u p =u . As it can be seen in ( 32), it is possible to define an effective ideal gas for the suspension writing

p p p p m m m , T r p p , e e e , u u , T T , , p p                      We deduce       2 2 2 2 2 2 v v p p M K 1 a / u M , K 1 a a K 1 r r , K 1 K 1 , K 1 K 1 c c c c , K 1 m m , K 1                             (33)

The two-phase stationary shock wave as an interface

One can now write the equations of the shock wave 0II considered as a whole. Using the equivalence between the suspension at equilibrium (states 0 et II) and an ideal gas, the equations for the shock 0II become:

          0 0 0 0 0 2 2 0 0 0 2 2 1 2 1 2 u u m with u u u p u p u p u p u II II p g II II II II II II              , : (34)
The following usual relations may be written:

    2 0 2 0 2 0 2 0 0 II 2 0 2 0 M K 1 c / u M , , p / p p with 1 M 1 2 1 p ~1 M 1 2 1 1 u ~                          0 II v / v v v (35)
The entropy production rate is here

                       II 0 0 II v a p p ln c m  (36)
For weak shock waves the small parameter corresponding to ( 16) becomes

1 1 M 2 0     (37)
and gives then for the entropy source

    3 2 v a 1 3 1 2 c m          (38)
Equation ( 38) is similar to equation ( 17). The foregoing remarks stay valid. In particular, in (38), there is no separation between thermal and viscous dissipations. For small mass fraction K of particles, the Mach number 0 M is near the Mach number M 0 . Then, the shock waves 0I and 0II are weak at the same time. We have

  K K 1       
It appears then a singularity for the magnitude orders of  and K. At a given small value of K, the behavior of the flow will depend of the relative orders of magnitude of  and K. The asymptotic analysis of the transonic regime can be performed with the small parameter . From (33) we deduce

   3 2 a a K K 1 1 K 2 1 K 1 K 1                              (39) Putting p 1 K    
, with p>0, one has

p 1 K K    
, and there are three possibilities :

1) 1>p>0 we can take for example p=1/2 2) p=1 3) p>1

we can take p=2

The following approximate formulas are deduced in each of these cases

                  3 3 2 v a a 3 2 2 2 2 v a a 2 2 2 v a a K 1 3 1 2 c m ) 3 K O 3 3 K K 1 3 1 2 c m ) 2 K O K 3 1 3 1 2 c m ) 1                                           (40)
It can be noticed that the cases 1) and 3) are particular cases of 2). Indeed, if K0 or K<<<<1, 2) gives 1). If <<K<<1, 2) gives again 3). Equation of case 2) is then the most significant approximate equation of the problem.

In first approximation, one has then

      3 2 v a K 1 3 1 2 c m            (41)

The interfacial behavior of the relaxation zone

Let us now consider the analysis of the relaxation zone III (Prud'homme, Jamet & Gottesdiener 1997[START_REF] Jamet | Modélisation des discontinuités dans les écoulements diphasiques compressibles[END_REF]. It is assumed that the shock 0II results from a gaseous shock 0I (through which the particles remain unchanged) followed by a relaxation zone III. The entropy production of the relaxation zone is

0 a a a       (42)
Two curves which give  a from (42) as function of M 0 , are put in Figure 5. One is obtained with the exact formulas ( 14) and ( 36) and the other with the approximate ones, deduced from ( 17) and ( 41).

In this model of shock wave in a suspension, we have assumed that a gaseous shock wave was already present. Consequently, the only values of  to consider are positive because the gaseous Mach number must be larger than unity. In addition, the Clausius-Duhem inequality must be applied to the shock waves 0I and 0II, and to the relaxation interface III. Consequently all the corresponding entropy productions must be non-negative. Let us consider now the generalized forces and fluxes. Equations ( 15) and ( 35) allow the connection between the parameter  and the various velocity jumps. One has

  K 1 u 2 u u , K 1 u 2 u u , 1 u 2 u u 0 II I 0 II 0 0 I 0                 
and the entropy production of the relaxation interface can be rewritten

1 II 0 I 0 2 II 0 2 I 0 S        A (44)
The conjugate generalized force is then

  S 2 0 2 v S K u 1 c m 2 A J      (45)
Here, we can see that the linearized TIP remains valid for the relaxation interface if we choose for A S the definition (44): the flux-force relation ( 45) is well linear.

For the limiting case 1), one observes that K<< and u II u I A S (u 0 -u I )/2 and  a 3K 2 << a .

In the case 3), K>> and A S (u 0 -u II )/2 and  a  3 K 3 .

Figure 5 Reduced entropy production v a c m /   as function of distance to sonic case M 0 =1 [START_REF] Jamet | Modélisation des discontinuités dans les écoulements diphasiques compressibles[END_REF]

Internal structure of the relaxation zone

Evolutions of the flow variables inside the relaxation layer can be determined numerically [START_REF] Marble | Dynamics of dusty gases[END_REF], and Figure 6 shows the obtained results for M 0 =1.6, far from the transonic regime. Entropy production due to the momentum and heat exchanges between gas and particles is given by ( 29) and one has

                           II I p p II I a dx u u T F T 1 T 1 Q dx (46)
Using ( 27) -( 29) one obtain well

  ] ] s [ [ m dx u u T F T Q ], ] s [ [ m dx T Q II I II I p p p p                ( 47
)
but this is not a separation between heat conduction and viscous effects. This separation is given by

                  II I II I p av p aT dx u u T F , dx T 1 T 1 Q (48)
and the ratio  aT / av is depending of the chosen constitutive relations in the bulk of the suspension. With the linearized relations (30) and for a stationary flow, one obtains

            II I II I 2 p v av 2 p p T p aT dx u u T l m K , dx T T T T l c m K   (49)
or, utilizing (32)

                      II I II I 2 p v av 2 p p T p aT dx dx u d T l m K , dx dx T d T T l c m K   (50)
From ( 32) we deduce the relation

I , 0 , E 2 u T c m 2 p                ( 51
)
We recall that T pI =T p0 =T 0 , u pI =u pO =u 0 and that T pII =T II , u pII =u II . For weak shock waves, and as in Section 2.3, if one admits that, in first approximation, relation ( 51) remains valid inside the whole relaxation zone, one has

p 2 p p p c u du dT , Cte m E 2 u T c , dx u d dx u d , dx T d dx dT         (52)
We obtain for the ratio of the entropies of the bulk

  2 0 v T p 2 v T p 2 v 2 p T p v T M 1 l l Tp T T c u l l T T dx u d T l m K dx dT T T l c m K                      (53) Assuming 1 M T T 2 0 p
 and with a ratio v T l l equal to Pr 2 3 for spherical mono-dispersed particles [START_REF] Kuentzmann | Aérothermochimie des suspensions Mémoire des Sciences Physiques[END_REF], where Pr is the mean Prandtl number of the gas, one obtains an evaluation of the relative dissipation effects inside relaxation zone of a weak wave in a dusty gas

  Cte Pr 2 1 3 v T av aT          (54)
From ( 43)-( 45) and ( 54) it is then possible to calculate the phenomenological coefficients L vS and L TS of the relaxation interface. For such an interface one has, in the same manner as for the volume in ( 29)

vS vS TS TS a J A J A   
. We take for A vS the value of A v given by ( 44) and J vS =L vS A vS . One find easily for Pr=2/3

  2 0 v p vS u 1 c m 2 L      (55)

LINEARIZED TWO-PHASE INTERFACES

First we set conditions of the problem and write interface equations for the suspension. The thermodynamic theory is then recalled. Interfacial transfer coefficients of momentum and energy are finally expressed.

In the second part we apply this method to the case of stationary shock wave propagation in a suspension and we treat relaxation zone behind shock wave as an interface. The comparison of the resulting data with those obtained previously in the transonic zone by asymptotic analysis permits in particular to deduce appropriate values for the Interfacial transfer coefficients.

General assumptions and balance equations

General assumptions

More details than in section 3 are given here on the considered simplified suspension.

Fluid phase is considered as a perfect gas except in a viscous boundary layer around each particle. In the same way thermal conduction is neglected except in a thermal boundary layer around particles. The various properties of fluid like dynamic viscosity, thermal conduction coefficient and specific heats are taken constant. Inertia, gravity, electrostatic and van der Waals forces due to Brownian motion of the particles are neglected. Thermal exchanges except conduction are not taken into account.

Particles are mono-dispersed, spherical and rigid and there are neither collisions nor hydrodynamic interactions between them. At last, mass transfer between gas and particles is not considered. For numerical examples, one takes micrometric carbon particles. As mean free path of oxygen or nitrogen molecules is about 5 10 -8 m (at normal pressure and temperature conditions) and size of particles is more than 10 times larger, each phase may be considered as a continuum and we make assumption that mean values for each variable can be defined everywhere in the flow. Every variable relative to solid phase will be subscript with a ( p ) and no subscript is affected to fluid variable. Specific density of particle and fluid are written  p s and  s and solid and fluid densities in the suspension are written  p and . One has =(1-c) s with c the particle concentration in suspension (we shall take = s as c is very small) and  p =c p s =m n p where m=(4/3)  3  p s mass of a solid particle and n p number of particles per unit volume of suspension.

The perfect gas verifies the equation of state: p= r T where p is pressure, T fluid temperature and r=R/M where R is the molar universal gas constant and M is the gas molecular weight; e=c v T is the fluid specific internal energy with c v specific heat at constant volume assumed as constant and s=(e+p/-g)/T is the specific entropy with g the free enthalpy per unit mass of the fluid. In the same way: e p =c s T p and s ps =(e ps -g ps ) for dispersed phase.

Moreover, Reynolds number based upon relative velocity between gas and particles is small compared to unity i.e.

1 d u u Re p       
, so we shall take stokes law for the drag. For the bulk, with the above assumptions, the equations ( 23)-( 25) remain valid in the case of a one dimensional flow.

Interfacial balance equations

For an interface one need to have balance equations of the type (1) in the more general case, or (2) in the present stationary case. The production term will vanish for conservative quantities. Then, applying (2) to total mass, momentum and energy balance of dilute simplified suspension lead to the following relations for a planar interface between two sides of a one dimensional flow (Figure 7)

                                  0 ] ] u p 2 u e m 2 u e m [ [ 0 ] ] u m u m [ [ 0 ] ] m m [ [
] ] s m s m [ [ a p p       (57)
Interfacial equations can be written for each phase separately. In the bulk it is generally more easy to write at first mass, momentum and internal energy balance laws for dispersed phase because they can be directly deduced from the basic equation obtained for one particle by multiplication by the number of particles [START_REF] Prud'homme | Notions de base sur les équations des brouillards[END_REF]. We will use the same method for interface.

It gets for dispersed phase

          a p p a p p p Q ] ] e [ [ m F ] ] u [ [ m 0 ] ] m [ [    (58)
where characteristic exchange terms F a and Q a of the interface have been introduced in (58) with the assumption that these terms keep the same form as in the bulk. The kinetic energy and total energy of the dispersed phase can be deduced in the usual form

pS a a 2 p p p pS a 2 p p u F Q ] ] 2 u e [ [ m u F ] ] 2 u [ [ m       (59)
if, and only if, we take for the interfacial velocity u pS the mean value

2 u u u pII pI pS   (60)
The balance equations for gaseous phase are deduced from ( 56) and ( 58) by subtraction

                  pS a a 2 a u F Q ] ] u p [ [ ] ] 2 u e [ [ m F ] ] p [ [ ] ] u [ [ m 0 ] ] m [ [    (61)
The interfacial fluid velocity is not determined a priori and can be chosen. By analogy with (60) we set for u S the mean value of u

2 u u u II I S   (62)
Let us finally give the expression of the interfacial entropy production rate. Assumption has been made that disequilibrium between phases was small, so equations can be linearized. Expression of entropy jump is then (using Gibbs relation)

S p p S S S T ] ] e [ [ ] ] s [ [ ] ] [ [ T p T ] ] e [ [ ] ] s [ [    v ( 63
)
where v is gas volume per unit mass   1

v

. Interfacial parameters T S , T pS and p S are mean values of temperatures and pressure that will be set below. Then it gets for the entropy production

a pS S 2 2 p p a S S p p S S S a Q T 1 ] ] [ [ p m ] ] u p [ [ ] ] 2 u [ [ m ] ] 2 u [ [ m Q T 1 T ] ] e [ [ m ] ] [ [ T p T ] ] e [ [ m                            v v      ( 64
)
Some terms of this relation can be transformed as follow with the help of ( 59) and ( 61). This lead to the intermediate result

    ] ] u [ [ p ] ] pu [ [ u ] ] p [ [ u u F T 1 T 1 T 1 Q S S pS S a S S pS a a                 ( 65
)
where pressure terms can be rearranged

  0 2 p p p u u ] ] u [ [ p ] ] pu [ [ u ] ] p [ [ II I S I II S S              (66)
by the above definition of the interfacial pressure taken equal to the mean pressure.

Finally the interfacial entropy source becomes

  pS S a S S pS a a u u F T 1 T 1 T 1 Q              ( 67
)

Thermodynamical theory

We recall now some results of the theory of linearized Thermodynamics of Irreversible Processes (TIP). The interfacial rate of entropy production  a is a sum of terms. TIP tells us that every term can be considered as the product of a generalized force by a generalized flux where a force is the cause of the phenomena and a flux is the result [START_REF] De Groot | Thermodynamics of irreversible processes[END_REF]. In other definition [START_REF] Gatignol | Mechanical and thermodynamic modelling of fluid interfaces[END_REF] a force is unchanged if the time is changed in its opposite and a flux is changed in its opposite in the same operation. We will work here with the first definition and consider that heat transfer is due to temperature difference between particle and gas and that viscous force is caused by their velocity difference. Heat flux Q a and drag F a appear then as generalized fluxes. Besides, when disequilibrium is small, fluxes and forces can be connected linearly as follow

                                   pS S S
u u L T 1 T 1 L T F u u L T 1 T 1 L Q (68)
where the L i jS are interfacial phenomenological coefficients. Onsager showed that the matrix of coefficients L i j is symmetric when there is no magnetic field [START_REF] De Groot | Thermodynamics of irreversible processes[END_REF]. Finally the Curie principle tells us that coupling between different tensorial orders phenomena is not possible. Here the two studied phenomena are respectively scalar (heat transfer towards the particles) and vectorial (drag) so that the relations between generalized fluxes and forces are simplified

  pS S vS S a S pS TS a u u L T F , T 1 T 1 L Q             (69)
The above coefficients L TS =L 11S and L vS = L 22S traduce exchanges in interface and will be calculated for the particular case of "relaxation interface" inside two-phase shock wave.

      .
Finally thermodynamic fluxes are given as functions of thermodynamic forces as in section 4.1, and then Q a and F a can be written as follow Finally, the system of mass, momentum and energy jump equations can be rewritten introducing the interface phenomenological coefficients which appears in the right hand side of (74e) has been neglected in this linearized theory. Equations ( 74d) and (74f) can be written in a form analogous to the one of system (32)

    pS S vS a pS S TS a u u F , T T Q     L L (
                                     ) f ( T T ] ] e [ [ m ) e ( ] ] u [ [ p T T ] ] e [ [ m ) d ( u u ] ] u [ [ m ) c ( u u ] ] p [ [ ] ] u [ [ m ) b ( 0 ] ] m [ [ ] ] u [ [ ) a ( 0 ] ] m [ [ ] ] u [ [ pS S
TS pS S p vS pS S p l T T ] ] T [ [ , l u u ] ] u [ [     (75)
with the non-dimensional coefficients

TS p s TS vS p vS m c l , m l L L     .
Coefficients L vS and L TS (or l vS and l TS ) can be determined using the shock relations and (74d) and (74f). It gets for L vS pI I

I p II p p vS u u u u m 2     L ( 76
)
where u pI =u p0 =u 0 and u pII =u II . This leads to

      2 0 2 0 p 2 0 2 0 p vS M 1 M 1 K 2 1 1 2 m M 1 K 1 M 1 K 2 1 1 2 m                       L (77) When M 0   it gets    K 2 1 K 1 1 2 m l p vS 1 vS               L (78)
Values of this reduced limit coefficient are given in Coefficient variation with K is rather small but one can say that, as K represents fraction of particles in the flow, when K goes up, momentum transfer is higher so coefficient p vS m  L evolves in the same way. For M I higher than 0.8, a quick divergence of coefficient is observed, divergence that is quicker when K is higher (Inversely l vS tends to zero). This divergence prevents us from characterizing momentum exchange in the whole zone. This problem of transonic divergence has been treated above. It will appear also with the interfacial heat transfer coefficient as we will see bellow and we need to match the linear and non linear expansions.

M M 2 1 1 M 2 I 2 I 2 0         (
Let us now determine coefficient l TS . One has from (74f)

I p I pI II p s p TS T T T T c m 2     L ( 80
)
where e p has been replaced by c s T p with c s =c p = r /(-1), T S -T pS by (T I -T pI )/2 and where T pI =T p0 =T 0 and T pII =T II .This leads to

        2 0 2 0 2 0 2 0 2 2 s p TS M 1 1 M K 1 1 M 1 M K 1 K 2 1 1 c m 2                  L (81)
This coefficient L TS can be written in a non-dimensional form by dividing by the unit mass flow rate of the particles multiplied by their specific heat and we obtain then l TS

- 1 . When M 0   one gets        2 2 s p TS 1 TS K 2 1 K 1 K 1 1 2 c m l                L ( 82
)
Values of this reduced limit coefficient are given in table 2 for =1.4.

In the general case, if we consider state I as the reference state for the relaxation interface I  II, we must express L TS as a function of M I Conclusions for this coefficient are almost the same that for coefficient l vS -1 that is to say that l TS -1 is nearly constant (2) until M I reaches a value of about 0.75, beyond this value the coefficient grows quickly and finally diverges when M I tends to unity. The fact that l TS -1 grows up with K seems logical too because when K gets higher, fraction of particles is more important and so energy transfer is increased. Table 2 Limiting values of the reduced interface transport coefficient for internal energy

Reduced transport coefficient

Matching with transonic results

Let us compare the values of l vS obtained by the present linear analysis and those deduced from the coefficient L vS in the transonic region.

We have, in the transonic region and with the linearized theory, for

  1 M K 2 0  <<1 and K<<1 2 1 1 M K M 2 1 K 1 2 1 2 1 l 1 _ 2 0 2 0 vS                        ( 84
)
The obtained values are the same.

It must be remarked that (44) gives the good value of the interfacial generalized force for momentum A vS inside and outside the transonic region. Generalized force for temperature A TS can be written in the same manner

                                    A ( 85
)
The main difference between linearized and non-linearized theory is that, in linearized theory applied to interface of relaxation, the generalized forces of (69)

                I 0 S pS TS 0 I pS S vS T 1 T 1 2 1 T 1 T 1 , 2 u u u u A A
are only depending on initial differences between particles and gas, while in transonic analysis generalized forces (44) and ( 85) are also functions of the final values (at state II). But this difference between these analysis disappears since   1 M 2 0  becomes sufficiently large (but not too much to remain in near equilibrium situation) compared to the particle concentration K. Then, we can take for phenomenological coefficients the values obtained by linearized theory for large Mach numbers that come from (78) and ( 82 It can be noticed also that the entropy production obtained in the linearized theory is not in accordance with those obtained by the asymptotic method in the transonic region for p=1. But if we make a truncation of (40-2) and write

      2 2 v a 3 3 K 1 3 1 2 c m              ( 87
)
the obtained entropy production agrees with those of the linearized theory [START_REF] Jamet | Modélisation des discontinuités dans les écoulements diphasiques compressibles[END_REF]) as shown in Figure 10.

CONCLUSION

Weak shock waves in a pure gas and in a dusty gas have been considered as interfaces and studied with the help of Thermodynamics of Irreversible Processes.

Only the relaxation zone in dusty gas is a matter for linear irreversible thermodynamics. For Mach numbers very near unity, there is a singularity due to the assumption of very small concentrations of particles and we have performed asymptotic expansion to solve this problem . For larger Mach numbers a linearized theory of interfaces is valid.

Generalized forces and fluxes have been deduced from the expression of entropy production and phenomenological coefficients L vS and L TS have been determined for interfacial heat and momentum transfer between particles and carrier gas. These interfacial transfer coefficients are both proportional to unit mass flow rate of the dispersed phase with proportionality coefficients depending on interfacial temperatures, isentropic coefficient of gas, specific heats and concentration of particles. The values of non-dimensional coefficients l vS and l TS are both near 1/2. These results are available for dilute suspensions and weak stationary shock waves.

The theory has to be extended in the other cases and the effects of wave front curvature and non-stationary flow must be analysed. The influence of larger particle concentrations is certainly a problem to be solved. Shock phenomena are only examples of generalized interfaces in multiphase media and other situations should be studied.

Figure 2

 2 Figure 2 Interface and interfacial zone a)Interfacial zone (at small scale) b) Interface of discontinuity (the same as a) but at larger scale)

Figure 3

 3 Figure3Structure of the shock wave[START_REF] Smolderen | Structure des chocs et théorie cinétique des gaz[END_REF]

Figure 4

 4 Figure 4 Structure of the shock in a two-phase Flow velocity

  the entropy balance equation and the entropy production

Figure 6

 6 Figure 6 Spatial structure of shock relaxation zone M 0 =1.6, =1.4, K=0.25, l T =l v[START_REF] Marble | Dynamics of dusty gases[END_REF]

Figure 7

 7 Figure 7 Planar interface in a one-dimensional flow

  79) Reduced transport coefficient (per unit mass flow rate of the dispersed phase) p vS m  L appears to be only function of the Mach number M I and the mass flow rates ratio K, for a given value of . It is put in Figure 8, in function of M I for a set of values of K. This coefficient remains more or less constant (2) for a Mach number M I less than 0.8 and for all values of parameter K.

Figure 8

 8 Figure 8 Reduced interface transport coefficient for momentum as a function of Mach number upstream the relaxation interface

  only function of the Mach number M I and the mass flow rates ratio K, for a given value of . It is put in Figure 9, in function of M I for a set of values of K.

Figure 9

 9 Figure 9 Reduced interface transport coefficient for internal energy as a function of Mach number upstream the relaxation interface

Figure 10

 10 Figure 10 Interfacial entropy production: linearized and formula (85).

table 1

 1

	for =1.4

This last expression does not vanish, neither for u 0 =u I (which gives  a =0) nor for u 0 =u II (which gives 0 a  

). Then the generalized flux to consider can be neither u 0 -u I nor u 0 -u II . At equilibrium, both generalized forces and fluxes must vanish together with the entropy production [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]. In fact,  a never vanishes for any value of the upstream Mach number M 0 (Figure 5). That means that the interfacial relaxation is always out (even if it is near) equilibrium. The generalized flux is certainly a combination of the various velocity jumps which appear in (43). A convenient generalized force may be

Application to the relaxation zone inside a two-phase shock

Let us consider the stationary shock in a one-dimensional two-phase flow defined in section 3. Focus will be made on the relaxation zone between states I and II, zone of few millimeters width for micrometric particles and considered here as an interface. This system sustains two disequilibria, the first, internal to the gaseous shock between 0 and I, and the second downstream this shock between I and II. In order to apply results from section 4.1 to this problem, one has to determine values of variables at states I and II. This has been made in sections 2.2 (system 15) for the conditions in state I and in section 3.4 for the conditions in state II. One recalls that for the jump 0I one gets from ( 15)

with M 0 =u 0 /a 0 , a 0 =(rT 0 ), r=R/M,

, R the universal gas constant and M the molar mass of the gas, and that for the jump 0II the same relations are valid with the values of the equivalent ideal gas described in section 3.3

As relaxation zone is very thin, it is considered as an interface and jump relations associated to balance equations can be written as above to express interfacial entropy production rate, if the previous assumptions are valid. We then have