Colin Leclercq

Benoît Pier

Julian F Scott

Absolute instabilities in eccentric Taylor-Couette-Poiseuille flow

The effect of eccentricity on absolute instabilities (AI) in the Taylor-Couette system with pressure-driven axial flow and fixed outer cylinder is investigated. Five modes of instability are considered, characterised by a pseudo-angular order m, with here |m| 2. These modes correspond to toroidal (m = 0) and helical structures (m = 0) deformed by the eccentricity. Throughout the parameter range, the mode with the largest absolute growth rate is always the Taylor-like vortex flow corresponding to m = 0. Axial advection, characterised by a Reynolds number Re z , carries perturbations downstream, and has a strong stabilising effect on AI. On the other hand, the effect of the eccentricity e is complex: increasing e generally delays AI, except for a range of moderate eccentricites 0.3 e 0.6, where it favours AI for large enough Re z . This striking behaviour is in contrast with temporal instability, always inhibited by eccentricity, and where left-handed helical modes of increasing |m| dominate for larger Re z . The instability mechanism of AI is clearly centrifugal, even for the larger values of Re z considered, as indicated by an energy analysis. For large enough Re z , critical modes localise in the wide gap for low e, but their energy distribution is shifted towards the diverging section of the annulus for moderate e. For highly eccentric geometries, AI are controlled by the minimal annular clearance, and the critical modes are confined in the vicinity of the inner cylinder. Untangling the AI properties of each m requires consideration of multiple pinch points.

Introduction

The flow between rotating cylinders has attracted attention since the end of the 19 th century, starting with the experiments of Couette (1888b,a) and [START_REF] Mallock | Determination of the viscosity of water[END_REF], and the landmark work by [START_REF] Taylor | Stability of a viscous liquid contained between two rotating cylinders[END_REF], who first predicted theoretically the threshold for centrifugal instability. Taylor characterised centrifugal effects using a non-dimensional number appropriate in the limit of small clearance, d = (b -a) ≪ a, with a and b the inner and outer cylinder radii (see figure 1). In this paper, a wide gap geometry with radii ratio η = a/b = 0.5 will be considered, and centrifugal effects will be conveniently measured by an azimuthal Reynolds number Re Ω = aΩd/ν, with Ω the inner cylinder rotation rate and ν the kinematic viscosity.

Adding axial flow to this system, one obtains a simple prototype for the study of pattern formation in real open flows. The effect of axial advection can also be measured by a Reynolds number Re z = W d/ν, based on the mean axial velocity W . First theoretical predictions of the absolute instability (AI) threshold were obtained by Tsameret & Steinberg (1991a) with a criterion based on a one-dimensional Ginzburg-Landau equation (with coefficients determined by two-dimensional numerical simulations), and then by [START_REF] Babcock | Noise-sustained structure in Taylor-Couette flow with through flow[END_REF][START_REF] Babcock | Stability and noise in Taylor-Couette flow with through-flow[END_REF], using the full set of hydrodynamic equations and a saddlepoint criterion [START_REF] Briggs | Electron-stream Interactions with Plasmas[END_REF][START_REF] Bers | Space-time evolution of plasma instabilities-absolute and convective[END_REF]) that will be discussed in §2.3. They showed that upon crossing the AI threshold, periodic self-sustained vortices appear, in contrast with the irregular patterns emerging from noise amplification in the convectively unstable régime. In these papers and subsequent work (Tsameret & Steinberg 1991b;[START_REF] Babcock | Stability and noise in Taylor-Couette flow with through-flow[END_REF][START_REF] Lücke | Amplification of molecular fluctuations into macroscopic vortices by convective instabilities[END_REF][START_REF] Swift | Effects of thermal noise in Taylor-Couette flow with corotation and axial through-flow[END_REF][START_REF] Babcock | Noise amplification in open Taylor-Couette flow[END_REF][START_REF] Tsameret | Absolute and convective instabilities and noise-sustained structures in the Couette-Taylor system with an axial flow[END_REF], effort was dedicated to identifying the noise sources (inlet noise versus thermal noise) which sustain the convective instability (CI). These studies were restricted to small axial Reynolds numbers Re z , typically below 4, and it was found that the most unstable (fastest growing) perturbations were in the form of propagating Taylor vortices. For higher values of axial advection, [START_REF] Takeuchi | Numerical and experimental investigation of the stability of spiral Poiseuille flow[END_REF] and [START_REF] Ng | On the linear stability of spiral flow between rotating cylinders[END_REF] had previously shown numerically (and also experimentally for the former reference) that critical modes consisted in propagating helical vortices, with helicity opposite that of the basic flow, and with azimuthal order m increasing with Re z . However, the concept of AI was not widespread in fluid mechanics back then, and these studies were restricted to CI. Theoretical prediction of AI of helical modes was investigated only recently (Hoffmann 2004[START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF]. In these papers, the authors studied the effect of axial through-flow on the spatio-temporal properties of toroidal and helical wavepackets with angular orders |m| 2, for Re z 20. For |m| 1, it was shown that critical azimuthal Reynolds numbers Re Ω are higher for AI than for CI, and that the difference between the two thresholds increases with Re z . For |m| = 2 and a stationary outer cylinder, AI was found to occur in a closed region of the Re z -Re Ω plane, considering only the saddle point originating at the critical conditions for CI with Re z = 0 (detailed in §2.4). However, the authors mentioned other saddle points expected to destabilise these modes in other regions of parameter space and which will be taken into account in the present article. More recent work on AI in the Taylor-Couette-Poiseuille flow concerned the effect of radial flow at the inner cylinder, representative of filtration devices [START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF]. In this analysis, it was shown that axisymmetric modes become absolutely unstable for inward radial flow, while helical modes with helicity identical to that of the basic flow dominate at high enough Re z , for outward radial flow.

When the two cylinder axes do not coincide, axisymmetry is broken and the stability properties of the flow are modified. Eccentricity is generally measured by the nondimensional distance between the two cylinders e = c/d (see figure 1). Adding eccentricity to the Taylor-Couette flow with axial advection, one obtains a basic model for annular mud flows in oil-well drilling, or lubrication flows present in high-speed journal bearings. In the first case, mud is injected in a rotating drillstring, and flows back to the surface through the annular domain between the drillstring and the rock face, with several engineering functions: carry the rock cuttings out, lubricate, prevent inflow of formation gases and wellbore collapse, etc. [START_REF] Escudier | Fully developed laminar flow of purely viscous non-newtonian liquids through annuli, including the effects of eccentricity and inner-cylinder rotation[END_REF][START_REF] Guo | Applied Drilling Circulation Systems -Hydraulics, calculations, models[END_REF]. For deep wells, the drillstring inevitably bends along its axis, on a typical length scale much larger than the well diameter. As a result, a parallel-flow assumption is reasonable, and the flow can be locally described as a Taylor-Couette-Poiseuille flow between eccentric cylinders. In high-speed turbomachinery, a similar configuration is found: oil is contained in eccentric journal bearings for lubrication purposes, and a pressure gradient is imposed along the shaft to evacuate damaging impurities [START_REF] Sep | Journal bearing with an intensive axial oil flow-experimental investigation[END_REF]. Aside from its fundamental interest, these industrial applications motivate the present analysis. In both applications, transition to complex hydrodynamic régimes would result in increased frictional losses, detrimental to the system efficiency. If the basic flow advection is weak compared to the rotation rate, hydrodynamic resonance may occur and the entire flow would bifurcate to an undesired self-sustained oscillatory state. This specific behaviour, called absolute instability, is particularly 'dangerous', because it does not require a permanent forcing: once the instability is triggered, it will propagate in both the downstream and upstream directions, and amplify using energy from the basic flow. On the other hand, convective instabilities correspond to wavepackets propagating only in the downstream direction: in the absence of forcing, the system eventually relaxes to its initial state at any fixed location, after perturbations have been 'blown away' from the source. The most temporally amplified perturbations are given by a classical temporal stability analysis, and such a study was recently carried out for this flow [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]. It was shown that the physics is essentially similar to the axisymmetric case [START_REF] Takeuchi | Numerical and experimental investigation of the stability of spiral Poiseuille flow[END_REF][START_REF] Ng | On the linear stability of spiral flow between rotating cylinders[END_REF], with propagating toroidal vortices replaced by helical structures of increasing azimuthal complexity as Re z is increased. Eccentricity deforms the critical modes, but does not introduce new instabilities to the problem. The effect of eccentricity is stabilising for all values of Re z , and this result is interpreted as a consequence of the reduction of centrifugal effects in the basic flow. Indeed, as eccentricity increases, the azimuthal flow rate decreases for a fixed inner cylinder rotation rate, resulting in weaker driving of the instability. This is a consequence of the appearance of a recirculation region in the wide gap which does not contribute to the net azimuthal flow rate. Outside this zone, in the vicinity of the inner cylinder, the flow resembles an axisymmetric Taylor-Couette flow with clearance scaling with the inner gap d(1 -e). The reduction of the 'effective' clearance ratio δ = d/a with eccentricity, or increase in 'effective' radius ratio η, also explains stabilisation (see [START_REF] Diprima | The stability of a viscous fluid between rotating cylinders with an axial flow[END_REF] for the effect of η on the Taylor-Couette flow). To date, it is, to the authors' knowledge, the only available theoretical study of eccentric Taylor-Couette-Poiseuille flow. The only known series of experiments were performed by Coney & Mobbs (1969-70); [START_REF] Coney | Taylor vortex flow with special reference to rotary heat exchangers[END_REF]; [START_REF] Younes | The hydrodynamic stability of spiral flow between eccentric rotating cylinders[END_REF]; [START_REF] Younes | Hydrodynamic stability of the flow between eccentric rotating cylinders with axial flow: torque measurements[END_REF]; [START_REF] Mobbs | The Taylor vortex regime in the flow between eccentric rotating cylinders[END_REF]; [START_REF] Coney | The effect of Taylor vortex flow on the radial forces in an annulus having variable eccentricity and axial flow[END_REF] and show good agreement with our a posteriori predictions, despite small discrepancies attributed to finite-length effects. For a brief review of other theoretical and experimental results on eccentric Taylor-Couette flow on the one hand, and axisymmetric Taylor-Couette-Poiseuille flow on the other hand, we refer to [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF].

The present paper extends this previous linear stability analysis by considering the case of AI. In §2, the linear stability framework is presented. The governing equations and numerical methods are briefly described, and the main properties of the basic flow and normal modes are recalled. In §2.3, the methods used to investigate AI, based on the [START_REF] Briggs | Electron-stream Interactions with Plasmas[END_REF][START_REF] Bers | Space-time evolution of plasma instabilities-absolute and convective[END_REF] pinching criterion, are described. In §3, results are presented for the 5 modes of instability with angular orders |m| 2, which include the fastest growing temporally unstable modes for Re z 50.

Linear stability framework

In the following, the geometry will be described using the ratio 0 < η = a/b < 1 between the inner and outer cylinder radii a and b (see figure 1), and the eccentricity 0 e = c/(b -a) < 1, based on the distance c between centers, divided by the clearance d = b -a. The gap varies azimuthally between d(1 -e) and d(1 + e). The radii ratio will be fixed at the value η = 0.5 throughout this paper. Rotation and axial advection will be quantified using the two Reynolds numbers given in introduction: Re Ω = aΩd/ν and Re z = W d/ν, with Ω the inner cylinder rotation rate, W the basic flow mean axial velocity and ν the kinematic viscosity.

The velocity u will be made non-dimensional with the rotation speed V ≡ aΩ. The clearance d will be taken as the reference length scale L. Finally, the pressure p will be in units of P ≡ ρV 2 , with ρ the density of the fluid. All equations and physical quantities will be written in non-dimensional form, using V , L and P .

Basic flow

The velocity u can be decomposed into a component w parallel to the axis e z , and a component u ⊥ = u -we z in a plane perpendicular to the axis. The axial flow is driven by a pressure gradient G in the z-direction. Denoting the in-plane pressure gradient as ∇ ⊥ p = ∇p -Ge z , the incompressible Navier-Stokes equations read:

(∂ t + u • ∇) u ⊥ w = - ∇ ⊥ p G + Re -1 Ω ∇ 2 u ⊥ w , ∇ • u = 0,      (2.1)
with impermeability and no-slip boundary conditions on the fixed outer cylinder and on the inner cylinder, whose rotational velocity is 1. In-plane and axial derivatives can be separated, using convenient notations:

u ⊥ • ∇ ⊥ ≡ u • ∇ -w∂ z , ∇ 2 ⊥ ≡ ∇ 2 -∂ 2 zz and ∇ ⊥ •u ⊥ ≡ ∇•u-∂ z w.
Assuming an axially invariant flow, the problem is two-dimensional and u ⊥ becomes independent of w:

∂ t u ⊥ + u ⊥ • ∇ ⊥ u ⊥ = -∇ ⊥ p + Re -1 Ω ∇ 2 ⊥ u ⊥ , ∇ ⊥ • u ⊥ = 0.
(2.2)

Basic flows Q ≡ (U, P ), denoted with capital letters, are defined as axially invariant, steady solutions of (2.1). Such solutions are found by integrating forwards in time equation (2.2) until convergence of U ⊥ is attained, and then solving for the corresponding axial velocity W , given by:

U ⊥ • ∇ ⊥ W = -G + Re -1 Ω ∇ 2 ⊥ W. (2.3)
Equations are expressed using locally orthogonal, body-fitted coordinates (ξ, φ), with -1 ξ 1 and 0 φ < 2π the pseudo-radial and pseudo-azimuthal coordinates respectively (see figure 2). In this modified bipolar coordinate system, a Fourier-Chebyshev pseudospectral projection method is implemented, with N φ = 2K φ + 1 Fourier modes, and N ξ Gauss-Lobatto collocation points. For more details on the numerical procedure, the reader is referred to [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF].

For an axisymmetric flow, the basic in-plane motion results from diffusion of axial vorticity from the rotating inner cylinder to the fixed outer cylinder. In cylindrical coordinates (r, θ), it takes the well-known form U ⊥ = (0, Ar + B/r), with A and B two constants depending on the geometry. For low eccentricities, the result is quite similar, as can be seen in figure 3(1a). However, for higher eccentricities, a low-speed recirculation region forms in the wide gap (figure 3(2a)). For the relatively high value of Re Ω = 500 presented here, small recirculation can already be seen for e = 0.2, whereas in [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF], figure 3, it was not present for Re Ω = 100 and appeared around e ≈ 0.3 for that lower value of Re Ω .

In the axisymmetric case, W is independent of U ⊥ , and the axial flow is very similar to a parabolic Poiseuille flow, with small corrections due the annular geometry. As eccentricity is increased, W decreases in the small gap, because of viscous effects, and most of the volume flux passes through the wide gap (see figure 3(2b)). Distortion also occurs, due to coupling with U ⊥ , and the peak velocity is no longer in the symmetry plane. For high rotation rates, the nonlinear interaction term U ⊥ • ∇ ⊥ W can locally dominate the viscous term Re -1 Ω ∇ 2 ⊥ W , and there is significant transport of W by in-plane components.

Normal modes

Let q ′ ≡ q-Q be three-dimensional perturbations of small amplitude superimposed onto the two-dimensional basic flow, and satisfying the linearized Navier-Stokes equations with no-slip boundary conditions. Because of temporal and axial invariance of the basic flow, perturbations are sought in the form of normal modes

(1a) (1b) (2a) (2b)
q ′ = q(ξ, φ) exp i(kz -ωt) + c.c., (2.4)
where c.c. denotes the complex conjugate. In a general framework, k is the complex axial wavenumber and ω is the complex frequency. As usual, ω r ≡ Re(ω) is the temporal frequency and ω i ≡ Im(ω) is the temporal growth rate. Equivalently, k r ≡ Re(k) is the wavenumber of the mode, and -k i ≡ -Im(k) is the spatial growth rate. Injecting the modal form (2.4) into the linearized Navier-Stokes equations with boundary conditions, one obtains a problem of the form (A -iωB)q = 0, where (A, B) are two linear operators, with A depending on k. Expressions for A and B are given in [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF], in the modified bipolar coordinate system. The wavenumbers k and frequencies ω satisfying this problem for non-zero q define the dispersion relation D(k, ω) = 0. Using the same Fourier-Chebyshev decomposition as for the basic flow, the linear problem is converted into a generalized eigenvalue problem for ω and q that can be solved numerically, using LAPACK (www.netlib.org/lapack) or ARPACK++ [START_REF] Lehoucq | ARPACK users' guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF] routines. For more information on the numerical procedure, the reader is referred to [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF].

RH2 RH1 TV LH1 LH2 m = -2 m = -1 m = 0 m = 1 m = 2
In this previous study, a temporal stability analysis was carried out to predict the fastest growing perturbations with k real and ω complex. It was found that among the large set of temporal modes ω(k), the most unstable ones were in the form of deformed toroidal vortices for low Re z , or complex helical structure for higher Re z . The modes were labelled according to a pseudo-azimuthal-integer-wavenumber m, or angular order. The labelling was made in accordance to the axisymmetric case, where normal modes can be written as q ′ = q(r) exp i(kz + mθ -ωt) in polar coordinates (r, θ) (see figure 2). Restricting attention to positive k, because of symmetry arguments to be discussed in the next paragraph, positive values of m (resp. negative) correspond to helical structures winding clockwise (resp. counter-clockwise) around the inner cylinder, and were called left-helical (resp. right-helical), or LH|m| (resp. RH|m|) modes. The case m = 0 corresponds to the classical toroidal Taylor vortex flow (TV). By following these modes as eccentricity is continuously varied, one obtains the corresponding pseudo-angular-order m for e = 0. Figure 4 shows the structure of modes m = -2, . . . , 2 for the classical Taylor-Couette flow.

Note that the symmetry Π 0 ≡ (m, ω, w) → (-m, -ω ⋆ , -w) (with ⋆ denoting the complex conjugate) between RH and LH in figure 4 is broken when axial flow is added, or when k is complex. Indeed, by taking the complex conjugate of the axisymmetric modal form, the general symmetry Π 1 ≡ (k, m, ω) → (-k ⋆ , -m, -ω ⋆ ) appears, also valid for e = 0. By considering the mirror image of the system (z → -z), one obtains another symmetry: Π 2 ≡ (Re z , k, w) → (-Re z , -k, -w). Combining Π 1 and Π 2 , one gets

Π 3 ≡ (Re z , k, m, ω, w) → (-Re z , k ⋆ , -m, -ω ⋆ , -w).
(2.5)

Setting Re z to zero and k real in Π 3 , one recovers Π 0 . In the general case, because of Π 3 , one can choose to study only m 0, or only Re z 0, without loss of generality. Π 1 also indicates that it is possible to restrict computations to k r 0.

Absolute instability threshold

Absolute instability occurs when the impulse response wavepacket is temporally growing at any fixed axial position z. The asymptotic dynamics of the wavepacket is dominated by the normal mode which satisfies the pinching criterion of [START_REF] Briggs | Electron-stream Interactions with Plasmas[END_REF][START_REF] Bers | Space-time evolution of plasma instabilities-absolute and convective[END_REF] and has the largest temporal growth rate. This mode has zero group velocity ∂ω/∂k(k 0 ) = 0 for the complex absolute wavenumber k 0 . This condition indicates the presence of a saddle point of ω(k) at k 0 . This saddle point is associated to a branch-point singularity at the complex absolute frequency ω 0 = ω(k 0 ). Additionally, the saddle point must comply with causality and result from the collision between downstream-and upstream-propagating spatial branches, respectively denoted k + (ω) and k -(ω). The flow is absolutely unstable if and only if the most unstable pinch point has positive absolute growth rate ω 0,i . Otherwise, the flow is either stable or convectively unstable. In the latter case, the wavepacket grows while being advected away from the impulse location so the system eventually relaxes to its initial state at any axial position. For a comprehensive review of these concepts, the reader is referred to [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF]; [START_REF] Huerre | Open shear flow instabilities[END_REF]; [START_REF] Chomaz | Global instabilities in spatially developing flows: non-normality and nonlinearity[END_REF].

The border of the absolutely unstable domain is determined by following all the neutrally stable saddle points in parameter space. This is done by performing Newton-Raphson iterations at each point, varying simultaneously k r , k i and Re Ω until |ω i |, |Re(∂ω/∂k)| and |Im(∂ω/∂k)| are all below 10 -6 . Estimated values for the independent variables are obtained by linear extrapolation with respect to the parameter being varied, e.g. e, Re Ω or Re z . For |m| = 2, critical curves display folds, and it is necessary to implement a continuation scheme based on an arclength variable [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF].

However, not all saddle points are valid and only the ones satisfying the pinching criterion are relevant. In order to discard invalid saddle points, extensive tests are carried out, where the two spatial branches k(ω) coalescing at k 0 are tracked numerically as ω i is increased from ω 0 . The saddle point is a genuine pinch only when the spatial branches separate into the upper and lower half k-planes for large enough ω i . Indeed, causality demands that this be true for ω i > ω i,max , where ω i,max is the maximum temporal growth rate for real k. Spatial branches are obtained by numerically inverting the relation ω(k) with a Newton-Raphson iteration.

Validation

To validate the numerical procedure, critical curves in the axisymmetric case were computed and compared with literature results. For m = 0 and 1, [START_REF] Pinter | Spiral and Taylor vortex fronts and pulses in axial through flow[END_REF] provide the coefficients of fourth-order polynomials f m fitting their data in the range -20

Re z 20, with step δRe z = 1. The same procedure was applied here, and our calculated values gm at the same points were fitted by polynomials g m . To compare our results, the residual 20 -20 |f m (Re z ) -g m (Re z )| 2 between the two fits was divided by the residual 20 -20 |g m (Re z ) -gm (Re z )| 2 between our fit and our calculated values. For m = 0 and 1, this ratio is respectively 1.16 × 10 -2 and 1.10 × 10 -2 , showing agreement with the authors' calculations. For m = 2, only graphical data was available, and figure 5 shows excellent agreement with the 'island' of instability found by [START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF] in Re Ω -Re z plane. In their analysis, the authors considered only the saddle points originating at the real critical wavenumber k c of temporal instability with Re z = 0, as will be explained now. For Re Ω above the temporal instability threshold Re Ω,c , the medium is unstable, and the growing part of the wavepacket is bounded by two spatio-temporal rays referred to as leading and trailing fronts, respectively z/t = V + and z/t = V -with V -< V + . The fronts are defined by the conditions ∂ω/∂k(k and ω i (k ± * ) -V ± k ± * ,i = 0 (see [START_REF] Huerre | Open shear flow instabilities[END_REF] for more details). At critical conditions for temporal instability, Re Ω = Re Ω,c and k = k c , the constraint ∂ω i /∂k = 0 for k real defines two degenerate fronts propagating at the group velocity V max of the most rapidly amplified temporal mode: V ± = V max = ∂ω/∂k(k c ). For Re Ω just above Re Ω,c , V + = V - so the fronts are properly defined and the now complex wavenumbers k ± * are close to k c so V ± ≈ V max . In general, V max > 0, so the wavepacket is advected downstream and the flow is only convectively unstable. But when Re Ω is further increased, one front may eventually change propagation direction, which translates into the saddle-point condition ∂ω/∂k = 0 defining the AI threshold. In [START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF], only the two fronts bounding the convectively unstable wavepacket at Re Ω slightly above Re Ω,c and Re z = 0 were considered. These specific fronts were followed as Re Ω and Re z were varied, and the AI boundary in figure 5 corresponds to points where one of the fronts was stationary. Saddle points corresponding to other stationary fronts were ignored in figure 5, even though the authors mentioned the existence of more. Note that validation for m 0 is sufficient because of the Π 3 -symmetry (2.5).

In this study, only values of |m| 2 will be considered and a small number of Fourier modes K φ = 8 is deemed satisfactory, as can be seen in table 1. A higher number of collocation points N ξ = 16 is however required for accuracy at large Re Ω .

Results

A parametric study has been performed within the ranges e 0.7, 0 Re z 60 and 0 Re Ω 500. Within these bounds, modes m = 0, 1, 2 are always the most temporally unstable, except for a small range e 0.3, 50 Re z 60, where m = 3 has the largest temporal growth rate [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]). It will be assumed that the absolute growth rate of these modes will be higher than that of |m| 3. However, right helical modes RH1 and RH2 will be retained in the analysis, as RH are known to be more absolutely unstable than TV and LH in some cases (e.g. high Re z , outward radial flow, e = 0, cf. [START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF]). Because we are considering both positive and negative m, it is unnecessary to consider negative Re z , because of the Π 3 -symmetry (2.5).

Reference saddle point

The bifurcation to (pseudo-)toroidal vortices without axial flow is a steady one, ω r (m = 0) = 0, so CI and AI thresholds, respectively denoted here Re Ω,c and Re Ω,c-a , coincide in this case. For modes m = 0, CI occurs through Hopf bifurcations at Re z = 0, so AI only occurs above a higher threshold: Re Ω,c-a > Re Ω,c . For these modes, Re Ω,c-a (Re z = 0) is found by locating the saddle point with k 0 closest to k c , the real critical wavenumber of CI. This neutral saddle point corresponds to a stationary front of the impulse response wavepacket for Re Ω just above Re Ω,c , as explained in detail in section §2.4. For m = 1, 2, the stationary front is the trailing one: V -= 0, k 0,i < 0. For negative m, k 0,i > 0 because of the Π 3 -symmetry (2.5), and the stationary front is the leading one: V + = 0. The present paragraph defined the reference saddle point for each m, obtained for Re z = 0 and e = 0. These saddle points are systematically followed in parameter space to define critical curves of AI. However, as will be seen in the next paragraph, other saddle points are also relevant to the spatio-temporal dynamics and must be considered.

Multiplicity of saddle points

Pinch points corresponding to other stationary fronts can be identified using the geometric method described in §3.2 of [START_REF] Juniper | The effect of confinement on the stability of two-dimensional shear flows[END_REF]. The same approach is used here to study the dispersion relation of mode m = 0 in the presence of strong axial advection, Re z = 60, for three eccentricities around e = 0.3. Figure 6, similar to figure 2 in Juniper ( 2006), shows isocontours of ω i (k) for complex values of k. The thick black line indicates a contour in the complex k-plane including the steepest descent paths of all genuine pinch points, here s 1 and s 2 (s 3 and s 4 are 'spurious' k -/k -saddle points). The impulse response can be obtained at any time using a classical inverse Laplace transform formula (see [START_REF] Huerre | Open shear flow instabilities[END_REF] for instance), which is easily evaluated at large time using this integration contour. The asymptotic response is indeed dominated by the pinch point of largest temporal growth rate, which here is s 2 . With this geometrical approach, genuine pinch points and invalid saddle points are easily identified.

If the eccentricity is varied by a small amount, the nature of the saddle points may change: valid saddle points may become invalid and vice versa. This is illustrated in figure 7: in case (b), s 1 and s 2 are the two pinch points, but in case (a), only s 1 is a pinch, and in case (c), s 1 , s 2 and s 3 are all three valid. Since the growth rate of each saddle point also varies with the control parameters, the dominant pinch may either be s 1 , s 2 or s 3 , and a careful analysis is required when parameters are varied.

Maps similar to figure 6 are drawn for each value of m, and a large number of saddle points are identified each time. Saddle points with growth rate close to the referencesaddle-point 's are systematically followed in parameter space. Over the whole range of parameters, three different 'pinch points' are found to be relevant for TV, two for LH1, LH2 and RH2, but surprisingly, just one for RH1, despite numerous candidate saddle points.

Critical azimuthal Reynolds number

In figure 8, the critical azimuthal Reynolds number is represented as a function of eccentricity, for Re z = 0, 10, . . . , 60. We start by describing the curves in terms of m. Solid lines correspond to m 0, and dotted lines are for m < 0. Changes of critical saddle point are indicated with a filled (resp. open) circle for m 0 (resp. m < 0). The most important result is that for any value of e and Re z , m = 0 always has the lowest critical Re Ω , followed by |m| = 1, and finally |m| = 2. This observation gives credence to the assumption that modes with |m| 3 can be ignored in the analysis. Curves 8(b) and 8(c) also prove the importance of considering both positive and negative m, as LH are not always more absolutely unstable than RH. Indeed, for low eccentricities, RH1 is slightly more unstable than LH1, but the converse is true for high eccentricities. The dynamics is even more subtle for |m| = 2, as LH2 are generally more unstable than RH2, except for a small range of eccentricities that varies with Re z .

The effect of axial advection is to stabilize all the modes. For high enough Re z , the critical Re Ω seems to increase almost linearly with Re z . The rate of increase is much stronger for low eccentricities than for high eccentricities, regardless of the value of m, as already mentioned. Critical Re Ω are typically one order of magnitude higher than Re z , which means that the inner cylinder must be rotated much faster than the mean axial velocity to have self-sustained oscillations. Indeed, axial flow prevents AI by carrying perturbations downstream while rotation amplifies them.

The effect of eccentricity is more complex than that of Re z . For low eccentricities, Re Ω increases slowly for m = 0, 1, but decreases for all other modes. For high enough e, all the curves have the same shape: Re Ω decreases before reaching a minimum and then increases again beyond this minimum. TV and LH1 display another similarity: they switch critical saddle point between 0.2 e 0.4, for high enough Re z . This change of saddle point coincides with the change in sign of the slope: critical Re Ω increases with e for the first saddle point, but decreases for the second one. As was already mentioned in paragraph §3.2, TV even changes pinch point twice for Re z = 60. On the other hand, the critical curves for RH1 are smooth, because they are obtained by continuously following a single saddle point. Finally, two saddle points define the critical curves of LH2 and RH2. For low values of Re z and e, the AI threshold is quite complex, and the curves display folds. Folds in the critical curves mean that there are finite ranges of AI in parameter space, surrounded by CI. This unusual behaviour has already been pointed out by [START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF] in the concentric case, but it seems important here to underline the fact that for high enough Re Ω , all the modes eventually become absolutely unstable, regardless of any 'island' of AI occurring at lower Re Ω (cf. figure 5).

Critical absolute wavenumber k 0,r

Figure 9 shows the evolution of the absolute wavenumber k 0,r , spatial growth rate -k 0,i and frequency ω 0,r associated with the dominant mode m = 0, at critical conditions. As before, curves are plotted as functions of e, for Re z = 0, 10, . . . , 60. The absolute wavenumber k 0,r (figure 9(a)) evolves in different ways below and above e ≈ 0.3. Below e ≈ 0.3, critical modes have longer wavelengths as Re z increases, spanning up to 6 times the clearance for Re z = 60 and e ≈ 0.3. When e is high enough, the trend is the opposite, and critical modes have shorter wavelengths as Re z increases. Below e ≈ 0.3, k 0,r is almost constant, or slightly decreasing with e, whereas above e ≈ 0.3, it is clearly increasing with e. For large enough e, the critical wavelength seems to be controlled by the smaller clearance d(1 -e). Small discontinuities in k 0,r around e ≈ 0.3 indicate a change of critical saddle point.

Comments on the curves obtained for other m, although not displayed here for clarity, can be made. First, curves of k 0,r for LH1 are very similar to those of TV, except that LH1 only changes saddle point once for Re z = 60 instead of twice for TV. Ignoring the change of saddle point, they are also qualitatively similar to those of RH1. For all m, k 0,r always increases with e for high enough eccentricities, showing that all modes scale with the small gap at critical conditions. For |m| = 2, large values of k 0,r up to 8 and more are obtained for low eccentricites as well, provided Re z is high enough. This behaviour is not found for other modes, where k 0,r is always between 1 and 3.5 when e 0.3. For LH2 and RH2, a large discontinuity of axial wavenumber is observed upon switching saddle point. For example, k 0,r of LH2 varies from about 7.5 to about 3.5 for Re z = 60 and e ≈ 0.6.

Critical absolute spatial growth rate -k 0,i

The absolute spatial growth rate -k 0,i (figure 9(b)) measures the 'steepness' of the stationary front of the impulse response wavepacket. For e 0.3, it is slightly increasing with e. It varies quickly for Re z 20, and then slowly varies in the range 2.5 -k 0,i 4 for higher Re z . For e 0.3, -k 0,i increases with e, reaches a maximum value, and then decreases again. For high enough e, -k 0,i increases steadily with Re z . Similar trends are noticed for LH1 and RH1, with comparable ranges of values. For LH2 and RH2 however, the curves are quite different. For low eccentricities, -k 0,i increases significantly with Re z , whereas it is almost constant for other m. Therefore, extreme front steepness occurs for LH2 at low e, with -k 0,i > 12 for Re z = 60.

For RH1 and RH2 at low Re z , -k 0,i can be negative over the whole range of eccentricities. Physically, this means that the stationary front is the leading one in this case V + = 0 (cf. section §2.4), and that the most temporally unstable RH wave has a negative group velocity.

Finally, a common feature of all m is that -k 0,i has a maximum as a function of e.

3.6. Critical absolute frequency ω 0,r

The absolute frequency ω 0,r (figure 9(c)) is given here for reference. After nonlinear saturation of the instability, self-sustained oscillations with frequency close to ω 0,r are expected for a supercritical transition. Therefore, the values of ω 0,r can be used as a good estimate of the hydrodynamic resonance frequency of the flow, and may be of interest for engineering applications. The trends of the curves are very similar to those obtained for the spatial growth rate. The frequency range is shifted towards higher values as m increases. For RH1 and RH2, ω 0,r can be negative (always the case for RH2), indicating that the absolute phase speed c 0 ≡ ω 0,r /k 0,r of the mode is negative in this case. Finally, discontinuities in the absolute frequency occur for all modes but RH1, because of changes in critical saddle point. Discontinuities are not clearly visible for m = 0 and the change of saddle point is 'smooth', as will be discussed later in §4.3.

3.7.

Absolute temporal growth rate ω 0,i maps Critical curves in §3.3 indicate the AI domain for each m. However, they do not indicate which mode will be the most absolutely unstable if Re Ω is above two or more thresholds. Indeed, the mode which bifurcates first as Re Ω is increased does not necessarily have the highest absolute growth rate ω 0,i for larger driving. Figure 10 represents isocontours of ω 0,i in Re Ω -Re z plane for m = 0, 1, 2, and e = 0, 0.2, 0.7. Curves for m = -1, -2 can be recovered upon applying the Π 3 -symmetry (2.5).

The first conclusion is that m = 0 remains the most unstable mode over the whole range of parameters. However, for e = 0.7, isocontours of m = 1 are very close to those for m = 0, and one may expect LH1 to become more absolutely unstable than TV for larger eccentricities and large Re z .

Isocontours of ω 0,i for m = 1 and m = -1 (using Π 3 -symmetry) are generally close. For low eccentricity, RH1 is always more unstable than LH1, but for high eccentricity, the converse is true. For a moderate eccentricity of e ≈ 0.5 (not shown here), the dominant mode depends on the specific values of Re Ω and Re z .

For m = 1, 2, some isocontours have discontinuous slopes, as a consequence of a change of dominant saddle point. In figure 10(3a), the 'island' of instability previously presented in figure 5 is shown to be connected to a 'continent', for higher values of Re Ω , via a change of critical saddle point. Indeed, in their analysis of m = 2 in the axisymmetric case, [START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF] restricted their analysis to the reference saddle point defined in §3.1, even though the authors mentioned the existence of other saddle points. Figure 10(3a) gives the complete AI boundary for this case. For higher eccentricities, the 'continent' of instability absorbs the 'island', and for e = 0.7, the saddle point associated to the 'island' is always sub-dominant. When considering the saddle point associated with the 'continent', LH2 is generally more unstable than RH2. On the other hand, for low e and very low Re z , RH2 can be more unstable than LH2 because of the 'island' of instability. For moderate eccentricities, the ordering depends on the specific values of Re Ω and Re z .

Discussion

In this section, we study the critical modes and the production of perturbation kinetic energy. We discuss the results and the instability mechanism in the light of these elements.

Critical modes

In figure 11(1,2), we examine the spatial distribution of the mode m = 0 at critical conditions for Re z = 60 and three representative eccentricities: e = 0.2, 0.4, 0.7. The three-dimensional distribution of axial velocity is shown, together with the corresponding time-averaged distribution of perturbation kinetic energy E = 1 2 ( ũ⊥ 2 + | w| 2 ) (the exp[-k i z] dependence of the mode amplitude is omitted).

Consider the polar angle θ of the maximum of energy, with respect to the inner cylinder center, indicated in figure 2. For low eccentricities or low Re z , the mode is localised in the wide gap, at positive θ. However, for larger e or Re z , this maximum is shifted upstream to the region -90 • θ 0 • . For e = 0.7, the mode concentrates in the vicinity of the inner cylinder, on the wide gap side, but has a radial extent scaling with the small gap. When e increases for Re z = 60, the ratio between the maximum, and the average value increases from less than 2 for e = 0 to almost 14 for e = 0.7, indicating confinement of the mode into a smaller region of the annulus. Finally, the only contribution of ũ⊥ to E accounts for 65% to 85% of the total perturbation energy when e increases from 0 to 0.7.

Similar behaviour is found for LH1 at Re z = 60. For RH1 however, the energy is less tightly concentrated. The maximum of E still occurs at negative angles, but does not go beyond -45 • in this case. It is concluded that the appareance of the peak of perturbation energy at large negative angles is associated with the change of saddle point occurring for both TV and LH1.

Surprisingly, the critical modes of AI and CI peak at completely different regions of the annulus for moderate eccentricities and high Re z . In comparison (cf. [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]), the maximum energy of the critical mode of CI for e = 0.5 is always localised at large positive angles. No direct comparison should be made with AI because critical modes of CI are obtained for different threshold values of Re Ω . It is nonetheless interesting to observe that modes can be localised at either positive or negative polar angles depending on the situation.

Finally, we enumerate characteristics common to all m. First, in-plane motion accounts for the larger contribution to the total perturbation kinetic energy of all m for Re z = 60. Also, as e approaches 0.7 for Re z = 60, all m tend to have similar distributions of energy, with strong localisation close to the inner cylinder, over a radial extent scaling with the small gap. In addition, we recall that the critical wavenumber of all modes takes on large values k 0,r ∼ 6.5 -8 when e = 0.7 and Re z = 60. These observations indicate that for high eccentricities, the critical modes for all m scale with the small gap.

Production of perturbation kinetic energy

To further investigate the instability mechanism, the production of perturbation kinetic energy is calculated. In the well-known Reynolds-Orr equation, the local rate of production of E is given by -u ′ • (u ′ • ∇U). Averaging in time and separating velocities into in-plane and axial components, one can define two contributions (the spatial growth of the mode amplitude along z is ignored again)

P ⊥ = - 1 2 Re {ũ ⋆ ⊥ • (ũ ⊥ • ∇ ⊥ U ⊥ )} and P z = - 1 2 Re { w⋆ (ũ ⊥ • ∇ ⊥ W )} , (4.1)
corresponding to the work of the Reynolds stresses against the in-plane and axial shear respectively. Expressions for the nonlinear terms in the modified bipolar coordinate system are given in [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]. Because the basic flow is axially invariant, these two terms are the only contributions to the production of E.

Figure 12 represents the distribution of P ⊥ and P z for m = 0 and Re z = 60, at e = 0.2, 0.4, 0.7, as in figure 11. Only positive contributions are shown in grey shades, as negative contributions inhibit temporal growth. Distributions of P ⊥ and E look very similar: P ⊥ is maximum in the wide gap for low e, then at negative polar angles for moderate e, then close to the inner cylinder on the wide gap side for e close to 1. The dominant contribution to P ⊥ comes from the Reynolds stress term involving the pseudoradial derivative of the azimuthal velocity V . The same calculation of P ⊥ and P z has been performed for the critical mode of the CI at Re z = 60 and e = 0.5 (m = 2). P ⊥ was also found to account for most of the kinetic energy production (81%), but the peak was located at a positive polar angle of 92 • .

Distributions of P z show larger contributions near the walls, where ∇ ⊥ W is larger. As eccentricity increases, production of kinetic energy close to the outer cylinder decreases, as the mode concentrates in the vicinity of the inner cylinder. For e = 0.7, P z peaks almost in the same region as P ⊥ and E, namely close to the inner cylinder in the wide gap. Integration of P ⊥ and P z over the annular domain for Re z = 60 indicates that inplane shear dominates the production of kinetic energy, P ⊥ always accounting for more than 85% of the total amount, and even more than 97% for e = 0.7. The dominant contribution comes from the term involving pseudo-radial variations of pseudo-azimuthal velocity, characteristic of a centrifugal instability.

Instability mechanism

The striking destabilisation occuring at moderate eccentricities remains difficult to explain even after examining the basic flow, the critical mode and the production of kinetic energy. The transition from a stabilising to a destabilising effect of eccentricity occurs around 0.2 e 0.4 for both TV and LH1, through a change of critical saddle point when Re z is large enough. The recirculation region appears for e ≈ 0.2 for Re Ω ≈ 500 so it seems that the change of behaviour occurs after the recirculation zone has reached a critical size. The basic flow is then significantly different from a classical Couette flow. Indeed, while energy of the mode mostly localises in the vicinity of the inner cylinder, it also partially spans over the recirculation region for moderate eccentricities (cf. figure 11(b)). Surprisingly, a region of the flow located at negative polar angles seems to drive AI for moderate eccentricities, whereas kinetic energy production always peaks at positive angles for CI at criticality. Arguments based on local stability of the flow, however tempting in a quest of explanation, should be avoided here because of strong non-parallelism of the basic flow in the pseudo-azimuthal direction. Indeed, assuming the flow locally parallel in φ leads to completely wrong predictions of instability thresholds of eccentric Taylor-Couette flow (DiPrima 1963;[START_REF] Ritchie | On the stability of viscous flow between eccentric rotating cylinders[END_REF]). The most temporally unstable velocity profile, theoretically located at φ = θ = 0 • , does not coincide with the location of maximum vortex activity found in the experiments [START_REF] Vohr | An experimental study of Taylor vortices and turbulence in flow between eccentric rotating cylinders[END_REF]). On the other hand, global analyses yield good results in this geometry [START_REF] Diprima | Non-local effects in the stability of flow between eccentric rotating cylinders[END_REF][START_REF] Diprima | The nonlinear calculation of Taylor vortex flow between eccentric rotating cylinders[END_REF][START_REF] Eagles | The effects of eccentricity on torque and load in Taylor-vortex flow[END_REF], showing the limits of the local approach. Therefore, localisation of the modes and kinetic energy production at moderate eccentricities should be regarded as a global property of the entire flow field.

The consecutive switchovers of dominant saddle point for e ≈ 0.3 and Re z = 60 occur very near collisions of the distinct pinches into third-order saddle points [START_REF] Davies | Conditions for absolute instability in the cyclotron resonance maser[END_REF] or 'super branch points' [START_REF] Healey | On the relation between the viscous and inviscid absolute instabilities of the rotating-disk boundary layer[END_REF], satisfying simultaneously ∂ω/∂k = 0, ∂ 2 ω/∂k 2 = 0 and ω i = 0. Indeed, s 1 and s 2 collide while being neutral for e = 0.2839, Re z = 50.115, Re Ω = 403.21 while s 2 and s 3 coalesce with ω i = 0 for e = 0.3032 and Re z = 43.188, Re Ω = 353.60. As a result, the saddle points swap dominance but the absolute wavenumber, frequency, and spatial growth rate vary almost continuously through the exchange (see figures 9a,b,c). Consequently, the spatial distribution of the critical mode is little changed and the physical reason for the switchover remains unclear.

As e approaches 1, all m tend to behave in a similar way. After reaching a minimum value, the critical Re Ω increases again as e becomes larger. Instability thresholds of all m 0, and even their respective absolute growth rate for any Re z -Re Ω combination, become close at high e and less sensistive to variations of Re z . More similarities are found by inspecting the critical modes. For Re z = 60, all m have absolute wavenumbers k 0,r in the range 6.5-8, indicating a similar length scale. Indeed, the kinetic energy of all these modes is localised around the inner cylinder, on a radial extent of the order of the smallest gap d(1-e), consistent with a small wavelength. The localisation is so strong that the difference between toroidal or helical structure of the modes is partially 'blurred'. These similarities between m at high e is reminiscent of small gap Taylor-Couette-Poiseuille flow, where critical thresholds associated with different m are very close [START_REF] Ng | On the linear stability of spiral flow between rotating cylinders[END_REF][START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]. Indeed, as eccentricity increases, the Couette-like flow associated with rotation of the inner cylinder scales as d(1 -e), curvature effects become less important, and m behaves more and more like a continuous real wavenumber, as when η → 1. Matching between localisation of the modes and the Couette-like part of the basic flow may also explain the weaker effect of axial advection on the instability thresholds. Indeed, for large eccentricities, most of the axial volume flux passes through the wide gap, and the maximum value of W is located in the recirculation region, far from the inner cylinder. Hence, the region where perturbations are most amplified is spatially separated from the region where they are most rapidly 'blown away'. This observation could explain why critical Re Ω are less sensitive to Re z for e close to one.

Convective versus absolute instability

We conclude this section by comparing the thresholds of CI [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF] and AI in the eccentric Taylor-Couette-Poiseuille flow. Figure 13 represents the critical Re Ω for CI (dotted lines), and for AI (solid lines), as a function of Re z for e = 0, 0.1, . . . , 0.7. For Re z = 0, the two thresholds for m = 0 coincide, as expected for a steady bifurcation. The critical Re Ω increases with Re z for both CI and AI, but with a much larger rate for AI. This was expected as axial advection tends to carry the perturbations away from the source, so a larger driving is required to reach AI. The critical mode is always m = 0 for AI, corresponding to closed pseudo-toroidal Taylor vortices, propagating when Re z = 0. For CI, LH modes of increasing m become critically unstable as Re z increases. Open circles indicate a change of critical m on the CI thresholds. The effect of eccentricity is clearly stabilising for CI, but the effect on AI is more complex. For high enough Re z , as eccentricity increases, the critical Re Ω of AI slightly increases for low e, but then decreases before reaching a minimum and increases again for larger values of e.

Conclusions

In this paper, an absolute instability analysis has been performed for the flow between eccentric cylinders, with rotation of the inner one and a superimposed pressure-driven axial flow. The ratio between cylinder radii was fixed at η = 0.5, a value representative of an oil-well-drilling configuration. Five modes of instability have been considered, m = -2, -1, 0, 1, 2, thus including the most temporally unstable modes up to Re z = 50 found in [START_REF] Leclercq | Temporal stability of eccentric Taylor-Couette-Poiseuille flow[END_REF]. Throughout the range of Re z and e considered in the present study, the mode with the largest absolute growth rate is always the pseudo-toroidal vortex flow corresponding to m = 0. Unlike the temporal growth rate, the absolute growth rate of left-handed pseudo-helical modes (m > 0) is not always larger than that of right-handed ones (m < 0).

Increasing Re z tends to hinder absolute instability because axial flow sweeps perturbations downstream. As a rule of thumb, the rotational velocity of the inner cylinder needs to be approximately one order of magnitude larger than the mean axial velocity to trigger absolute instability.

The effect of eccentricity is more complex and increasing e can result in destabilisation for large enough Re z and moderate eccentricities 0.3 e 0.6. In this case, the critical mode has a complex structure, and the production of kinetic energy peaks at a welldefined region of the annulus, located in the diverging gap region. Outside this range of eccentricities, increasing e has a stabilising effect, increasingly so as the limit of touching cylinders is approached.

The instability mechanism is purely centrifugal in nature and the critical-mode axial wavelength and radial extent scale as the smallest gap d(1-e). For large eccentricities, all the modes localise in the vicinity of the inner cylinder, and their spatial distributions become more and more similar. The effect of Re z on absolute instability thresholds becomes weaker, and the distance between thresholds associated to different m diminishes.

Overall, many valid saddle points were found for the different instability modes, and the critical pinch point switches upon varying the flow parameters. The physical interpretation of the switchovers of saddle point occuring near e ≈ 0.3 for m = 0 and Re z = 60 is unclear, because the saddle points are associated with modes having almost the same spatial distribution and spatio-temporal properties. Indeed, when they swap, the saddle points almost collide into third-order saddle points, where they would be impossible to distinguish.

We believe that the most crucial outlook of this work is additional experiments to confirm our findings on convective and absolute instabilities, since the last measurements made on such a configuration are older than the introduction of absolute instability theory to fluid mechanics [START_REF] Huerre | Absolute and convective instabilities in free shear layers[END_REF]! Even in the convectively unstable régime, the literature is very lean, and only one apparatus seems to have ever existed, whereas the domain of application is vast. On a theoretical viewpoint, it would be particularly interesting to investigate the properties of nonlinear global modes, partly based on local absolute instability properties (Pier et al. 2001), when eccentricity varies slowly along the axis. Indeed, bending of the long drillstring results in axially varying eccentricity. Moreover, such weakly non-parallel open flows are believed to be good candidates to confirm the potential existence of hat modes [START_REF] Pier | Fully nonlinear global modes in spatially developing media[END_REF], 2001), theoretically predicted for model equations, but yet to be identified in a real configuration.

Figure 1 .

 1 Figure 1. Eccentric annulus of radius ratio η = a/b = 0.5 and basic flow U = U ⊥ + W ez.

Figure 2 .

 2 Figure 2. Modified bipolar coordinate system (ξ, φ) fitting the eccentric annular domain. Polar coordinates (r, θ) are centered on the inner cylinder, with θ = φ = 0 along the line joining the cylinder axes.

Figure 3 .

 3 Figure 3. Basic flows for ReΩ = 500: (1) weak eccentricity e = 0.2, (2) high eccentricity e = 0.7. (a) Contours of equispaced in-plane streamfunction with superimposed U ⊥ profiles at θ = 0, π/2, π, 3π/2 (polar angle with respect to the inner cylinder). (b) Equispaced contours of W .

Figure 4 .

 4 Figure 4. Structure (isovalues of axial velocity) of the five modes of interest m = -2, . . . , 2, from left to right. For this example, e = 0, ReΩ = 68.19, Rez = 0 and k = 3.16, corresponding to the critical conditions for temporal instability of the axisymmetric Taylor-Couette flow (see DiPrima & Swinney (1985) for a review).

Figure 6 .

 6 Figure 6. Contours of temporal growth rate at criticality for m = 0, e = 0.3075 and Rez = 60 (ReΩ = 472.24). Saddles are indicated by white dots. The thick line indicates the steepest descent path associated with the pinch points s1 and s2. The white cross indicates a branch point ∂k/∂ω = 0 from which originates a branch cut (vertical dotted line). The temporal growth rate of s2 is larger than that of s1: the pinch s2 dominates the impulse response.

Figure 7 .

 7 Figure 7. Zoom of the bottom-left corner of figure 6 for the same value of Rez, three eccentricities around e ≈ 0.3 and their associated critical ReΩ; (a) e = 0.285, ReΩ = 471.65, s1 is the dominating pinch; (b) zoom of previous figure, s2 dominates; (c) e = 0.33, ReΩ = 468.00, s3 dominates.

Figure 8 .

 8 Figure 8. Critical azimuthal Reynolds number ReΩ versus eccentricity e for (a) m = 0, (b) |m| = 1 and (c) |m| = 2. Curves are drawn for Rez = 0, 10, . . . , 60. Solid lines are for m 0, and dotted lines for m < 0. Filled/open dots indicate a change of saddle point.

Figure 9 .

 9 Figure 9. (a) Absolute wavenumber k0,r, (b) spatial growth rate -k0,i and (c) frequency ω0,r for m = 0 at critical conditions, versus eccentricity e. Curves are given for Rez = 0, 10, . . . , 60.

  Figure10. Absolute growth rate ω0,i maps for e = 0, 0.2, 0.7, m = 0, 1, 2.

Figure 11 .

 11 Figure 11. Critical mode m = 0 for Rez = 60, (a) e = 0.2, ReΩ = 465.82, (b) e = 0.4, ReΩ = 432.74, (c) e = 0.7, ReΩ = 314.16. (1) Distribution of axial perturbation velocity Re( w). Dark (resp. light) grey is for positive (resp. negative) values. (2) Distribution of perturbation kinetic energy E = 1 2 ( ũ⊥ 2 + | w| 2 ). Dark grey indicates high values, and isocontours are equispaced.

Figure 12 .

 12 Figure 12. Distribution of kinetic energy production for m = 0 and Rez = 60, at critical conditions. (1) P ⊥ and (2) Pz. (a) e = 0.2, ReΩ = 465.82, (b) e = 0.4, ReΩ = 432.74, (c) e = 0.7, ReΩ = 314.16. Dark grey indicates large contributions and white corresponds to zero or negative contribution. Isocontours are equispaced.

Figure 13 .

 13 Figure 13. Absolute (solid lines) and convective (dotted lines) instability thresholds: critical ReΩ versus Rez, for e = 0, 0.1, . . . , 0.7. Open circles indicate a change of critical m for CI. For AI, the critical mode is always m = 0.

Table 1 .

 1 Critical azimuthal Reynolds number ReΩ for Rez = 50 and different resolutions N

	± * ) = V ±

AI Figure 5. 'Island' of AI in the axisymmetric case, for mode LH2. Line: present calculation; dots: figure 8(b) in Altmeyer et al. (2011). ξ × K φ .

† Email address for correspondence: colin.leclercq@ec-lyon.fr


Dragan Slavkov Hajdukovic
email: dragan.hajdukovic@cern.ch

Testing the gravitational properties of the quantum vacuum within the Solar System

The existence of the quantum vacuum is well established in the Standard Model of Particles and Fields but completely neglected in contemporary Astrophysics and Cosmology. Independently of any theory it is a major and absolutely urgent task of astronomical observations to reveal if and how quantum vacuum contributes to the gravitational field of baryonic matter in the Universe. We point out that a first signature of the gravitational impact of the quantum vacuum might be seen in the perihelion precession of orbits of satellites of minor planets in the outer part of the Solar System. As an example we consider the minor planet 2002 UX25 and its satellite.

Introduction

So far we had two scientific revolutions in our understanding of gravitation: Newton's law and Al e t Ei stei 's General Relativity. Whatever happens in the future, these two revolutions will remain among the greatest achievements of theoretical physics and the human mind.

We know today that both theories have in common a wrong assumption. The wrong assumption is that the matter of the Universe exists in classical, non-quantum vacuum.

While it is systematically neglected in Astrophysics and Cosmology, quantum vacuum is an inherent part (see for instance [START_REF] Aitchison | Nothi g's ple ty-The vacuum in modern quantum field theory[END_REF][START_REF]Measurement of the running of the fine-structure constant[END_REF][START_REF] Wilson | Observation of the dynamical Casimir effect in a superconducting circuit[END_REF]) of the Standard Model of Particles and Fields which is in perfect agreement with experimental findings. However, without quantum vacuum taken into account, quantum field theory would be in perfect disagreement with experimental findings.

Einstein said (and I agree): Imagination is more important than knowledge. So, I invite you to imagine that, you can switch off and switch on, the quantum vacuum in our Universe. As we live in the Universe with quantum vacuum switched on, you firstly must switch it off. What would happen? In fact I must warn you not to do it; after switching off the quantum vacuum, you will not stay alive to switch it on again! This is not a speculation but prediction based on our best knowledge. For instance, sophisticated experiments [START_REF] Perez | The quark and gluon structure of the proton[END_REF] have revealed that the proton is not an elementary particle but a very complex system that in addition to three valence quarks contains virtual (or see) quark-antiquark pairs and gluons (see Fig. 1).

Fig. 1 Inner structure of a proton revealed at HERA Black spirals represent gluons while purple-green particles denote virtual quark-antiquark pairs (up to 100 of these quark/anti-quark pairs a e isi le at any instant!). Note that there are three more quarks (two up, one down) than anti-quarks. These are the three valence quarks we would normally refer to when speaking of the proton. (Source: DESY in Hamburg)

In simple words, quantum vacuum significantly contributes to the structure of protons (and neutrons as well). If this contribution is switched off, protons would become quite different particles. A radical change of constituents of atoms would perturb everything; the Universe without quantum vacuum would be a completely different place (and certainly without us). Hence, quantum vacuum is not only a strange state of matter in quantum field theory, but also the root of our existence.

Before we continue let us give a simplified (but basically true) description of the quantum vacuum. Quantum vacuum should be considered as a state of matter, completely different from familiar states gas, li uid, solid, plas a… ut as real as they are [START_REF] Aitchison | Nothi g's ple ty-The vacuum in modern quantum field theory[END_REF][START_REF]Measurement of the running of the fine-structure constant[END_REF][START_REF] Wilson | Observation of the dynamical Casimir effect in a superconducting circuit[END_REF]. Popularly speaking, quantum vacuum is a sea of sho t li i g virtual particleantiparticle pairs (like quark-antiquark, neutrinoantineutrino and electron-positron pairs). According to our best knowledge: (1) quantum vacuum is a state with perfect symmetry between matter and antimatter; a particle always appears in pair with its antiparticle, which is totally different from mysterious matter-antimatter asymmetry, i.e. the fact that everything on the Earth (and apparently in the Universe) is made from matter, with only traces of antimatter; [START_REF]Measurement of the running of the fine-structure constant[END_REF] contrary to all other states of matter which are composed from the long living particles (electrons and protons in stars and flowers, have existed before them and will exist after them), the quantum vacuum is a state composed from extremely short living virtual particles and antiparticles (for instance, the lifetime of a virtual electron-positron pair is only about 22 10  seconds). Popularly speaking, as fishes, we live in an o ea ; ou o ea is the quantum vacuum.

Our knowledge about gravitational properties of the quantum vacuum is zero.

There are two possibilities. The first possibility is that quantum vacuum has no impact (or at least has no significant impact) on the gravitational field in the Universe. In fact, contemporary Astrophysics and Cosmology are based on this assumption.

The second possibility is that quantum vacuum contributes greatly to the gravitational field in the Universe. If so, any theory that neglects the existence of the quantum vacuum is blind to some crucial gravitational phenomena, and, as a compensation for the lost phenomena, must invoke some artificial stuff.

The inevitable question is if dark matter and dark energy are such artificial stuff which (in our incomplete theory) mimics well the phenomena which are in fact caused by the quantum vacuum.

Contrary to the other candidates for new physics (supersymmetries, dark matter, dark e e gy… ua tu a uu is ot a spe ulatio but a key feature of Quantum Electrodynamics, Quantum Chromodynamics and Electroweak theory, i.e. a key feature of the Standard Model of Particles and Fields, the most successful and the best tested theory of all times. It is a major and absolutely urgent task of astronomical observations to reveal if and how quantum vacuum contributes to the gravitational field of baryonic matter in the Universe. In the present article we point out that a first signature of the gravitational impact of the quantum vacuum might be seen in the perihelion precession of orbits of satellites of minor planets in outer part of the Solar System. As an example we consider the minor planet 2002 UX25 and its satellite.

Some reflections on quantum vacuum and gravity

We do not know if quantum vacuum contributes or not to the gravitational field in the Universe. The answer must come from astronomical observations and experiments in our laboratories. However it seems plausible to assume that quantum vacuum contributes to the gravitational field.

The first question after this assumption is if the contribution of quantum vacuum is independent or dependent on the quantity and distribution of the immersed matter. Once again the answer must come from observations and experiments. From a theoretical point of view, if quantum vacuum produces a non-zero gravitational field independent of the immersed matter, it means that we should observe a non-homogeneous and non-isotropic background, which is apparently not the case. Hence, it seems plausible that the gravitational contribution of the quantum vacuum depends on the quantity and distribution of the immersed matter. In other words, if we attribute to the quantum vacuum, an effective gravitational charge density, it has a non-zero value only in the presence of an external gravitational field produced by the immersed matter. Let me note that I have used an unusual te g a itatio al ha ge instead of gravitational mass or simply mass; of course, according to the weak equivalence principle all these terms are synonyms for ordinary matter, but it is preferable to keep the term gravitational charge and to stay open for possibility of violation of the weak equivalence principle for non-ordinary matter.

Our theory of gravitation describes the gravitational field of a spherical body (like a star) immersed in the gravitationally featureless classical vacuum. Hence, the second basic question is how quantum vacuum modifies classical gravitational field around a spherical body. The gravitational charge (mass) of the considered spherical body is a constant. Contrary to it, the quantum vacuum around the body should have certain effective gravitational charge density; hence the effective gravitational charge within a sphere of radius r must increase with r . In other words an observer at a larger distance r from a star would observe a larger effective gravitational mass of the star. It is obvious from observations that the effective gravitational charge of the quantum vacuum within the Solar System must be a tiny fraction of the mass of the Sun. Only at very large distances, this fraction can be significant.

It is interesting to note that in Quantum Electrodynamics, quantum vacuum has properties analogous to the above suppositions.

The simplest gravitational signature of new physics

The simplest problem in celestial mechanics is to determine orbit of a point-like body in a central gravitational field. The orbit is an ellipse fixed with respect to the centre of gravity if and only if the central gravitational field has perfect spherical symmetry and the gravitational force strictly follows the Newton inverse square law. Any departure from spherical symmetry and/or the inverse square law of gravity, leads to the precession of the perihelion (see for instance [START_REF] Dermott | Solar system dynamics[END_REF][START_REF] Fitzpatrick | An Introduction to Celestial Mechanics[END_REF]).

General Relativity (describing a spherically symmetric central gravitational field by Schwarzschild metric) is more accurate than Newtonian theory and it predicts a tiny precession even in the case of spherical symmetry. The general relativistic precession is well approximated by It is worth noting that the perihelion precession of planets predicted by classical theory (Newtonian mechanics together with the inverse square law for gravity) is close to the observed values; the largest discrepancy occurs for the Me u y Me u y's o it p e esses at a rate that is about 8% greater than the predicted one). The discrepancy has been explained by general relativistic correction (1) which must be added to the Newtonian result.

In brief, in the case of a central gravitational field with the exact spherical symmetry, Newtonian precession of the perihelion is zero, while the general relativistic result (1) is too small to be detected (with the current accuracy of measurements) for satellites of small central bodies like minor planets in the outer part of the Solar System.

Of course, a minor planet (which orbits around the Sun with a period Sun T ) and its satellite (which orbits around the planet with a period P T ) are not an isolated system but subject to an external gravitational field dominated by the Sun. This external gravitational field produces a Newtonian perihelion shift in the orbit of the satellite: the shift per orbit is well approximated [START_REF] Dermott | Solar system dynamics[END_REF][START_REF] Urbassek | [END_REF] with

2 2 3           Sun P N T T   (2)
Any significant difference between observed precession and the expected value (2) must be considered as a signature of new physics. Typical value of the shift (2) for satellites of trans-Neptunian minor planets is a few tens of arc seconds per century. As we will argue, the shift caused by quantum vacuum might be greater by more than one order of magnitude.

However, independently of any theoretical argument it is of major importance to observe and see if there is any anomalous perihelion shift in the orbits of satellites of trans-Neptunian minor planets. With the existing infrastructure of satellites and telescopes it would be less expensive but potentially not less important than LHC experiments at CERN or detectors devoted to search for dark matter.

An illustration of hypothetical gravitational effects of quantum vacuum

Let us assume that a satellite (with mass m ) orbiting about a minor planet (with mass M) is subject to both Newtonian gravitational force and a tiny additional radial acceleration   4), ( 5), ( 6) and [START_REF] Urbassek | [END_REF] In the simplest but illuminating case,   r A qv has a constant value denoted by qv A . In this case the perihelion shift per orbit can be written [START_REF] Dermott | Solar system dynamics[END_REF]8] in the following way

  dE

  m M G a e A qv qv      2 2 1 2  (9)
Hence, the internal precession (9) caused by the quantum vacuum is mixed with the precession (2) induced by the external Newtonian gravitational field dominated by the Sun. Let us note that according to equations ( 2) and ( 9) quantum vacuum might be dominant only for small systems far from the Sun. Now, as an example let us consider the minor planet 2002 UX25; the needed data taken from Reference [9] are given in Table 1 According to Comparison of ( 10) and (11) shows that in this example the effect of quantum vacuum is strongly dominant (it is larger nearly two orders of magnitude). Hence, there is potential for a clear signal of new physics.

Comments

Of course in addition to 2002 LX25 there are other trans-Neptunian objects with satellites which might be appropriate for the proposed measurement of the perihelion precession. Logos-Zoe, Quaoar-Weywot, (66652)-Borasisi are just a few examples of interesting trans-Neptunian planet-satellite pairs. It seems that otherwise the most attractive candidate, system Eris-Dysnomia [8] is perhaps excluded because of very small eccentricity of the orbit of Dysnomia.

Apparently the best strategy is to focus firstly on the measurement of the perihelion precession of a single satellite (for instance the satellite of 2002 LX25). If the result reveals new physics, in the second step it would be necessary to measure perihelion precession of many more satellites in order to derive insight on how the acceleration   In principle, the hypothetical effect of the quantum vacuum can be revealed by the study of the orbits of inner planets, the satellites of inner planets or even artificial satellites. However in these systems the hypothetical effect is only a tiny fraction of the total perihelion precession.

While it is not directly related to the proposed observations let me finish with one interesting comment. According to the Standard Model of Cosmology we live in a Universe with three spatial dimensions and non-Euclidian geometry; the Big-Bang is considered not as an explosion in space but rather expansion of space (see for instance Reference [12]). However, it is also known that the Big-Bang can be considered as an explosion in a Euclidian space with four spatial dimensions [12] which is presumably not reality. The interesting and never asked question is if the quantum vacuum and matter immersed in it have the same number (three) of spatial dimensions. One interesting possibility is that quantum vacuum has four spatial dimensions; if so the Big Bang might be considered as both, expansion of three dimensional space and explosion in a four dimensional space. I am not claiming that it is so, but just underlining how rich the spectrum of possibilities is, and how important it is to keep open mind.

Appendix: Quantum vacuum and virtual gravitational dipoles

As we have strongly stressed in this paper, the proposed test is of high importance independently of any theory. The study of orbits of satellites of minor planets in the outer part of the Solar system is a sensitive test of the eventual existence of tiny non-Newtonian component of gravity.

In this Appendix we present a far reaching theoretical possibility: the quantum vacuum containing virtual gravitational dipoles might be source of a weak gravitational component that does 't follo the Ne to 's inverse square law.

Let us start with the hypothesis that quantum vacuum contains virtual gravitational dipoles, i.e. the gravitational charges of the opposite sign on a distance smaller than the Compton wavelength of the corresponding virtual particles. Apparently, the most natural is to attribute the hypothetical positive and negative gravitational charge respectively to virtual particles and antiparticles [13][14][15]. However, some caution is needed; we still have to learn a lot about the content of the quantum vacuum and it is possible that the hypothesis of the existence of gravitational dipoles is more robust than the identification of dipoles with virtual particle-antiparticle pairs. Whatever is the nature of gravitational dipoles, if they exist, a gravitational polarization density The simplest possible case of the gravitational polarization of the quantum vacuum is saturation i.e. the case when the external gravitational field is sufficiently strong to align all dipoles along the field. If all dipoles are aligned in the same direction, the gravitational polarization density has the maximal magnitude   max where   denotes the Compton wavelength of a pion (basically a quark-antiquark pair).Additionally, for a central body of mass M the region of saturation is a sphere with radius S R which can be approximated [11, 13, 14, ] with

   m M R S  (15)
In the case of UX25, S R calculated from the equation ( 15) is about three orders of ag itude g eate tha the satellite's se imajor axis. Hence, the satellite of UX25 is in the region of saturation and its orbit should be perturbed with the constant acceleration (13).

  rotation per orbit in radians, a the semi-major axis of the orbit, e the eccentricity of the ellipse and S R the Schwarzschild radius of the central body.

2 

 2 the classical celestial mechanics (see for instance the book of Murray and Dermott, 1999, page 55) in this problem characterized with spherical symmetry, the time evolution of the argument of pericentre  is determined by the equation f cos and dt as functions of the eccentric anomaly E :de otes a e age a gula velocity (or the mean motion) and T is the orbital period. Equations (

  e. the gravitational dipole moment per unit volume) can be attributed to the quantum vacuum.As well known, in a dielectric medium the spatial variation of the electric polarization generates a charge density charge density. In an analogous way, the gravitational polarization of the quantum vacuum should result in a gravitational bound charge density of the vacuum

  leads to the conclusion that in the region of saturation, quantum vacuum acts as the source of a constant acceleration

Table 1

 1

												and Equation (2)
	 	N		09 . 3		10  8	rad		0064 . 0 arc	sec	(10)
	while	 	qv	depends on the choice of qv A in
	Equation (9). As an interesting choice suggested
	recently [10-11], let us take
	A qv		673 . 6		10	11 	m	/	s	2	. It leads to
		qv 	23 . 0 arc	sec	(11)
												Table 1
	Parameters for 2002 UX25 and its Satellite
	UX25 Mass				1.25x10	20 kg
	UX25 Semimajor axis	42.869 AU
	UX25 Orbital Period	280.69 years
	Satellite Semimajor axis	4770 km
	Satellite Orbital Period	8.3094 days
	Satellite Eccentricity	0.17


Marcel Bawindsom Kébré

Francois Ouedraogo

Naon Betaboalé

François Zougmoré

Fabien Cherblanc

François Ouédraogo

Soil-Water Characteristic Curve Modeling at Low Water Content: Empirical and Semi-Empirical Approaches

Keywords: soil, water, arid, modelling, empirical, semiempirical

published or not. The documents may come

L'archive ouverte

INTRODUCTION

The water content in the unsaturated soil is a function of soil water chemical potential. This constitutive relationship is commonly referred to as the soil-water retention curve in soil sciences, and well known as the soil-water characteristic curve (SWCC) in geotechnical engineering. Soil-water characteristic curve appears to be a key hydrodynamic parameter in the fields of soil science and civil engineering.

Using the SWCC, one can compute or determines the hydraulic conductivity function [1], the volumetric variation of the soil [2], the shear strength in the unsaturated soil [3], and the attribution of the water content in subsurface [4]. While knowledge of the two extreme states of water in a soil (saturated and very dry conditions) seems well advanced, the unsaturated state at medium to low water content deserves further investigation. This condition is usually encountered in soil surface layers where strong interactions occur between soil and roots, microorganisms or building foundations. In arid and semi-arid zones, best knowledge of physical flow processes at low water contents would allow development of more effective tools for water resources management. Despite the relatively small amounts of water retained, accurate representation of soil water characteristic curves at the dry end is important for modeling biological processes including plant water uptake and microbial activity in arid environments [5]. Accurate measurement over the full range of water contents is a delicate task. Most or all methods for measuring the water potential have limitations, especially at intermediate and low values of the water content (pendular and hygroscopic states). Frequently the curve is completed by extrapolation, which can lead to instabilities in numerical codes [6]. Experimental methods (infiltration and/or drainage, evaporation and condensation) associated with inverse problems also suffer from shortcomings, including in the instrumental devices, and are mostly robust only at the higher water contents (i.e., the funicular state) where the liquid phase is continuous. It arises from this the question of the validity at low water contents of the matrix headwater content relationship, h -, usually a key function to describe the physics of unsaturated flow. Indeed, this range of the water content comprises the pendular state where the liquid and vapor phases presumably are discontinuous, and the hygroscopic state where only the gas phase is continuous, while water is being strongly adsorbed onto the solid phase. This is the reason why one prefers the use of the water chemical potential to describe the movement of water in all of its states. One of the fundamental properties of the chemical potential stipulates that at local thermodynamic equilibrium, the potential of a component present in two phases should be the same [7]. Experimentally well-defined methods (e.g., [7]) should apply to all important parts of the soil water characteristic curve, from saturation to oven dryness. The main objective of this study is to evaluate the models representating the SWCC in the range of low water contents for an arid soil of Burkina Faso. These results will be used in the modeling of water transfer and management of water resources in arid areas. A comparison is made between the two modeling approaches in low water contents given the constraints of time and cost. This with a view to retain a relatively simple and less expensive method for arid soil hydrodynamic properties determination.

MATERIAL AND METHODOLOGY

TYPICAL SOILS FROM BURKINA-FASO

Morphological Characteristics

Two kind of natural soil from Nasso (Burkina-Faso) are under investigation. The first kind, denoted NH1, was taken from the top layer of ground (0-30 cm depth) and the second one, denoted NH2 comes from deeper layer (30-70 cm depth). Morphological properties are summarized in Tab. 1. Physical analyses have been realized by dry sieving after washing (standard NF P94-056) for particles larger than 80 μm and sedimentation analysis (standard NF P94-057) for smaller particles. Soils granulometric curves are ploted in Figure 1. According to ISSS classification, the textural typologies of these soils are sand for NH1 and loamy sand for NH2 what corresponds globally to sandy soils. Bulk density, s , of each sample have been measured using the Core Sampler Method on undisturbed soil (standard NF X31-501). Assuming the soil particle density to be d = 2650 kgm -3 leads to determine the saturated mass water content w sat :

w sat = r w 1 r s - 1 r d ae è ç ö ø ÷ (1)
where ρ w is the density of water, taken as 1000 kg.m -3 .

Measurement of soil water retention curve

The soil water retention curve was determined by combining measures from tensiometry and sorption isotherms.

Following the ideas put forth by [7], the water content is plotted as function of the water chemical potential, denoted μ in J.kg -1 , rather than pressure head or suction. Indeed, a variety of micro-scale mechanisms governs soil/water interactions and liquid pressure cannot be defined for low water content. From its energetic definition, the chemical potential can consistently describe the thermodynamic state of water over the whole range of water content. It can be seen as a unifying concept well defined from a thermodynamic framework. Classical experimental devices were used: (i) pressure plate apparatus at constant temperature (30 o C) for high water contents corresponding to the funicular state; (ii) saturated saline solutions for low water contents in the hygroscopic state. For both methods, at equilibrium, mass water content was calculated by differential weighting after oven-drying at 105 °C for 48 hours.

Figure 1 : Soil granulometric curves

THEORETICAL MODELLING

Semi-empirical approach

It consists to find an analytical expression that takes into account the variables that influence the flow of the water in the soil. The resulting models are currently based on the distribution of the pore size, distribution of particle size or relationships between measurable characteristics of soil. In formulations underwritten, μ is the absolute value of the water chemical potential.

(i) VG80 [1], where fitting parameters are: w r , α, n and m:

w = w r -w sat -w r ( ) 1+ am ( ) n é ë ù û -m
(2)

(ii) FS95 [10], with w a , α, n and m as fitting parameters and μ 0 = 10 6 Jkg -1 :

w = w a 1- ln m ln m 0 ae è ç ö ø ÷ + w sat -w a 1- ln m ln m 0 ae è ç ö ø ÷ é ë ê ù û ú 1+ am ( ) n é ë ù û -m (3) (iii) FX94 [12]
, where fitting parameters are w r , n, m, μ r and μ i :

w = 1- ln 1+ m m r ae è ç ö ø ÷ ln 1+ m 0 m r ae è ç ö ø ÷ é ë ê ê ê ê ù û ú ú ú ú w sat ln e+ m m i ae è ç ö ø ÷ n ae è ç ö ø ÷ é ë ê ê ù û ú ú m (4)
Numerical minimization procedure based on Levenberg-Marquardt algorithm was applied to objective function (Eq. 5), an ordinary weighted least-squares method, to identify optimal parameters values for each semi-empirical model.

O b ( ) = W i w i -ŵi b ( ) é ë ù û 2 i =1 N å (5)
where b is a vector whose components are all the unknown parameters of the model to identify, N is the number of data pairs for the retention function, W i is the class weight of the individual water content and w i , ŵi are the measured and model predicted values, respectively.

Empirical approach

Several functions are available in the literature. We chose to evaluate an empirical method for determining the soil water curve from basic geotechnical properties proposed by [16]: modified Kovács (MK) model. This model has been successfully evaluated by the authors on the basis of several types of plastic and non-plastic soils, especially in consideration of the thermodynamic states of water in soil: capillary water and adsorbed water.

The MK model is represented by a series of equations that have been developed considering that water is retained in granular materials by capillary forces, which are responsible for capillary saturation (S c ), and adhesion forces causing saturation for adhesion (S a ) [16]. Thus, according to the MK model, the degree of saturation S e in materials depends on two components, S c and S a . These components are originaly express in terme of matric head h (cm). One can show that the suction or matric potential ψ (in Pa) and the matric head h (in m), both widely used in soil hydrology, are closely related to the chemical potential μ (J.kg -1 ) [7]. The relationship between matric potential, chemical potential and matric head is given by μ=ψ/ρ w =gh, where ρ w is the density of water and g the acceleration of gravity (m.s -2 ).

S e = r s r e n w = 1-1-S c ( ) 1-S a (6) S c = 1- h co h ae è ç ö ø ÷ 2 + 1 é ë ê ê ù û ú ú m exp -m h co h ae è ç ö ø ÷ 2 é ë ê ê ù û ú ú (7) S a = a c 1- ln 1+ h h r ae è ç ö ø ÷ ln 1+ h 0 h r ae è ç ö ø ÷ ì í ï ï î ï ï ü ý ï ï þ ï ï h co h n ae è ç ö ø ÷ 2 3 e 1 3 h h n ae è ç ö ø ÷ 1 6 (8)
In equation ( 6), . represents the Macauley brackets ( x = 0.5 x + x ( ) ) and w (kg.kg -1 ), n (-) are the gravimetric water content and the porosity respectively. In equations (7-8), h co (cm) is an equivalent capillary height, a reference parameter in MK model. It is related to an equivalent pore diameter, with depend to the solid surface area. For granular materials, h co is obtained from the following relationship:

h co = 0, 75 eD 10 1+1,17 log C U ( ) é ë ù û (9)
where D 10 (cm) is the diameter corresponding to 10% passing on the cumulative grain-size distribution ; C U = D 60 / D 10 , the uniformity coefficient and e, the void ratio.

The parameter m in equation ( 7) is a pore size coefficient that controls capillary saturation and can be related to C U (m = 1/C U ). In equation ( 8), a c is an adhesion coefficient ; h n , a normalization parameter introduced for unit consistency ; h 0 , the matric head which induces water content of zero (h 0 = 10 7 cm for thermodynamic considerations) and h r (cm), the residual matric head. h r can be estimed with h co : h r = 0,86h co 1,2

(10) The MK model is easy to use. Just have the basic geotechnical properties namely bulk and soil particle densities to compute n and e; the particle size distribution to extract D 10 and D 60 .

RESULTATS AND DISCUSSION

SEMI-EMPIRICAL APPROACH

Saturated water content of each horizon of soil is evaluated through relation ( 1) with parameters values from Table 1. Tables 23-4 summarize the models evaluated for the water retention curve, the values of descriptive parameters identified for soils horizons NH1 and NH2. The values of statistical parameters determining the performance of the models in representing the experimental values as the coefficient of determination, R 2 , and standard deviation, RMSE (Root Mean Square Error), are recorded in this table. Based on these statistical values, the model FX94 is the most representative with R 2 mean = 99.84%, followed by the model FS95 with R 2 mean = 99.74% and finally the most popular model, VG80 where R 2 mean = 99.41%. However, in low water contents, differences are remarkable for the various models tested (Fig. 4) and (Fig. 5). For the sandy material NH1, FX94 is the most efficient while for sandy loam material NH2, the results are mixed. VG80 model, sometimes, fails to reproduce experimental values at low water content [8,10]. Forcing moisture content value to zero for a chemical potential of 10 6 J.kg -1 in FX94 and FS95 could be the cause for a soil with more fine particles.

In addition, the extensions models from FX94 and FS95 in low water contents support the hypothesis by [9]: a linear variation in semi-log between the water content and the chemical potential in the lower margin of the water retention curve. This assumption has restriction for materials containing more fine particles (clay soils). Several studies show that water retention properties in the dry end of the soil are strongly correlated with the specific surface area which increases with the clay content of the soil [5,16]. A proportional relationship between water content and logarithm of the chemical potential is often used to represent the water retention curve in this oven-dry condition [13,17]. Thus, the main challenge is to find a mathematical model, doubly differentiable, to connect the two extreme soil conditions (high saturation and very dry end) where available models seem effective in reproduction of experimental values of the water retention curve.

Another serious concern with regard to models of water retention curves, is the ease of use in a numerical problem of water transfers. When the relation μ = f (w) is used, it is very difficult or impossible to obtain an analytical expression derived from the water retention curve models from FX94 and FS95 [8]. This constraint greatly restricts the use of these models in simulation codes of water transfers, which is not the case for VG80. The model describes well the curve at the low water contents. At high saturation state, results are less representative. However, this model is a good basis for characterization of the water retention curve from data which can be acquired in a limited time and cost. His performance in low water contents are a major asset in the characterization of hydrodynamic properties of sandy soils in arid and semi-arid regions for optimization and management of water resources.

CONCLUSION

Two modeling approaches are used to determine the soilwater characteristic curves of two soil horizons from Burkina Faso. They represent typically the major part of soils in arid and semi-arid regions, regularly subjected to water stress increased by rare rainfall and extreme temperatures. Thus, particular emphasis is focused on the low water contents part, intermediate state between high water saturations and oven-dry states.

It appears from the study that Modified Kovacs model, an empirical approach model, is more appropriate to estimate the curve at low water contents compared to the semiempirical approach models evaluated. Combining terms from capillary state and adsorbed state of soil water gives a physical basis for the Modified Kovacs model. The use of semi-empirical models require enough data from the curve over the entire range of water content while the empirical model MK requires only basic geotechnical properties of soil and a few experimental points to adjust the parameters m and a c of the model. It would be useful to refine correlation functions between these parameters and the physical properties of soil from a large database covering the entire range of soil textural. Additional investigation may also be oriented on estimation of unsaturated hydraulic conductivity, the other soil hydrodynamic function, from the basic geotechnical properties.

Figure 2 :

 2 Figure 2 : Experimental and fitted theoretical SWCC: horizon NH1.

Figure 3 :

 3 Figure 3 : Experimental and fitted theoretical SWCC: horizon NH2.

Figure 4 :

 4 Figure 4 : Experimental and fitted theoretical SWCC: zoom at low water content for horizon NH1.

  Figure 5 : Experimental and fitted theoretical SWCC: zoom at low water content for horizon NH2.3.2. EMPIRICAL APPROACHThe Modified Kovàcs model was tested for the two nonplastic granular soil horizons NH1 and NH2. The values of basic geotechnical parameters (n, e, ρ s , D 10 , D 60 ), derived from the characterization of physical properties are used to calculate the characteristic curves from the set of equations (6) to (10). The important parameters for the model are the coefficient of adhesion a c acting in the of adsorption component S a and m, the pore size coefficient in the capillary state.[16] suggested the values a c = 0.01 and m = 1/Cu for granular materials (no plastic, no cohesive). Using these values for the horizon NH1 gave mixed results shown in Figure (6).

Figure 6 :

 6 Figure 6 : Experimental and fitted theoretical SWCC for horizon NH1 with empirical values of m and a c proposed by [16]. Using our experimental data, it is possible to adjust the parameters a c and m for each material tested to better represent the experimental points in any part of the curve, this is shown in Figure (7) for NH1 and Figure (8) for NH2. Four (04) experimental data regularly distributed on the curve are sufficient to adjust these parameters. The adjusted values are: a c = 0. 018, m = 0.031 for NH1 and a c = 0.0115, m =0.002 for NH2.

Figure 7 :

 7 Figure 7 : Experimental and fitted theoretical SWCC for horizon NH1 with empirical values of m and a c ajusted with experimental data.

Figure 8 :

 8 Figure 8 : Experimental and fitted theoretical SWCC for horizon NH2 with empirical values of m and a c ajusted with experimental data.

Table 1 :

 1 Physical properties of soils used

		Particle size distribution % Bulk
					Density
	Soil type Sand	Silt	Clay	kg.m -3
	NH1	90	7	3	1688
	NH2	85	6.5	8.5	1652

Table 2 :

 2

	VG80 SWCC parameters values
	w sat kg.kg -1	Model parameters n α kg.kg -1 kg.J -1 w r	-	m	-	R 2 %	Statistics RMSE kg.kg -1
	NH1 0.215	0	0.912	1.392 0.282	99.3 0.0048

ACKNOWLEGDEMENT

The authors would like to acknowledge Fondation 2iE (Ouagadougou -Burkina Faso, www.2ie-edu.org) for physical characterisation of materials in his laboratory (Laboratoire Eco-Matériaux de Construction -LEMC). We especially wish to François TSOBNANG, Adamah MESSAN and Koffi KOKOLE.