Joël Tchoufag

Jacques Magnaudet

David Fabre

Linear instability of the path of a freely rising spheroidal bubble

Keywords: drops and bubbles, instability, wakes/jets

Path and wake instabilities of buoyancy-driven oblate spheroidal bubbles with a prescribed shape rising freely in a viscous fluid otherwise at rest are studied using global stability analysis, following the technique recently developed for a coupled fluid + body system by Tchoufag, Fabre & Magnaudet (J. Fluid Mech. vol. 740, 2014, pp. 278-311). The essential role of the wake on the path instability is evidenced by comparing the shape of the global stability diagram with that obtained in the case of a fixed bubble. However, dramatic differences are also found, since the critical curve of the coupled system mostly involves low-and high-frequency oscillating modes, whereas that of a fixed bubble only involves stationary modes. Comparison of the present predictions with results obtained through direct numerical simulation is achieved in several regimes, confirming the predictions of the linear approach but also highlighting some of its limitations when the system successively encounters several unstable modes.

Introduction

Path instability of isolated rising bubbles has been a subject of curiosity for ages, as testified by Leonardo's drawings [START_REF] Prosperetti | Leonardo's paradox (Appendix to 'The added mass of an expanding bubble' by C. D. Ohl, A. Tijink and A. Prosperetti)[END_REF]Prosperetti 2004). Various hypotheses on the origin of this puzzling behaviour were suggested during the second half of the last century (see the review by [START_REF] Magnaudet | The motion of high-Reynolds-number bubbles in inhomogeneous flows[END_REF], including the effects of surfactants and the influence of shape oscillations. The picture was significantly clarified by the direct numerical simulations (DNS) of [START_REF] Mougin | Path instability of a rising bubble[END_REF] and [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF] with freely rising and fixed bubbles, respectively, which showed that path instability is closely related to wake instability. Both studies considered perfectly clean bubbles (i.e. the outer fluid obeyed a shear-free condition at the bubble surface) with a frozen oblate shape and revealed that the wake becomes unstable when the oblateness exceeds some critical threshold, the value of which depends on the Reynolds number based on the rising speed.

Independently, recent theoretical and computational studies devoted to two-and three-dimensional freely rising/falling rigid bodies such as plates [START_REF] Assemat | The onset of unsteadiness of two-dimensional bodies falling or rising freely in a viscous fluid: a linear study[END_REF], infinitely thin disks [START_REF] Auguste | Falling styles of disks[END_REF][START_REF] Chrust | Numerical simulation of the dynamics of freely falling discs[END_REF] or disks of finite thickness [START_REF] Chrust | Effect of solid body degrees of freedom on the path instabilities of freely falling or rising flat cylinders[END_REF][START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF] have revealed that the translational and rotational body degrees of freedom generally play a key role in the path instability mechanism. In particular, they can cause the thresholds and frequencies of the observed paths to differ drastically from those associated with the wake instability of the fixed body. This is why assessment of the possible relation and differences between the primary destabilization of the wake and path of a freely moving bubble and the wake instability of a fixed bubble experiencing a uniform upstream flow is in order. As discussed by [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF], this goal is best achieved by considering a global linear stability analysis (LSA) of the coupled fluid + body system. This is the approach followed here, the outline of which is given in § 2 before the main results are discussed in § 3. This section also provides some quantitative comparisons with the DNS results of [START_REF] Mougin | Path instability of a rising bubble[END_REF] (hereafter referred to as MM02) in various regimes.

Parameterization and formulation of the global stability problem

Following MM02, we assume the bubble to be a non-deformable spheroid characterized by the aspect ratio χ = b/a between the major and minor semiaxes. As shown in several experiments, e.g. [START_REF] Ellingsen | On the rise of an ellipsoidal bubble in water: oscillatory paths and liquid-induced velocity[END_REF], this assumption is approximately valid as far as surface tension effects are of the same order as inertial effects. The problem also depends on the bubble-to-fluid density ratio ρ, and on the Reynolds number Re = U 0 d/ν, where d = 2(b 2 a) 1/3 is the bubble equivalent diameter, U 0 is the equilibrium velocity of its steady vertical rise (or base state) and ν is the fluid kinematic viscosity. Nevertheless, comparisons with experiments and DNS are made easier by replacing U 0 by the gravitational velocity U g = ((1ρ)gd) 1/2 , which yields the Archimedes number Ar = U g d/ν, where g denotes gravity. The conversion from Re to Ar is achieved by determining the drag coefficient C D (Re) in the base state and equating the drag and buoyancy forces, which yields the relation C D (Re)Re 2 = 4Ar 2 /3. As is well known, Re and χ are not independent in real gas-liquid systems. Rather, given the fluid physical properties and bubble volume, they are uniquely determined by the Archimedes number and the Bond number Bo = ρ(1ρ)gd 2 /σ , where σ and ρ stand for the surface tension and liquid density, respectively. Hence, arbitrary variation of Re at a given χ or vice versa as we do here could only be achieved by independent variation of ρ, ν and σ in a laboratory experiment, which is generally not feasible.

Similarly to the formulation employed by [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF], we introduce a system of axes (x, y, z) translating and rotating with the bubble, with x aligned with the spheroid axis and y and z lying in its diametrical plane. The bubble orientation is defined by the vector Ξ whose components are the roll/pitch/yaw angles corresponding respectively to the angular positions of the x/y/z axes with respect to a fixed system of axes (x 0 , y 0 , z 0 ). The Navier-Stokes equations for the absolute velocity V and pressure P are written in the (x, y, z) reference frame; V falls to zero far from the bubble and obeys the zero-penetration condition V • n = 0 and the shear-free condition n × (T • n) = 0 at the bubble surface, where n denotes the unit normal to this surface and T = -PI + ν(∇V + t ∇V) is the stress tensor (I is the Kronecker tensor). The Newton equations expressing the overall force and torque balances on the bubble govern the evolution of its translational velocity U and rotation rate Ω. Finally, for small deviations from vertical, the evolution of the bubble orientation is governed by the equation dΞ /dt = Ω.

The corresponding LSA problem is solved thanks to the finite element solver FreeFem++ already used by [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF] for freely moving disks and by [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF] for fixed oblate bubbles. The reader is referred to these two references, especially the first one, for all technical details. As shown in [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF], dealing with the shear-free condition in the framework of a variational formulation requires the normal stress, Σ = n • T • n, to be considered as an additional unknown. The formulation described in [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF] is modified accordingly and the vector state corresponding to the fluid unknowns is

Q f = [V, P, Σ]
, where all variables depend on space and time. The complete state vector for the bubble

+ fluid system is then Q = [Q f , Q b ], where Q b = [U, Ω, Ξ ]
gathers the time-dependent bubble kinematic degrees of freedom.

In order to perform the global stability analysis, the state vector is split into a base state Q 0 and a linear perturbation q. The base state corresponds to the steady vertical rise of the bubble with zero inclination, i.e. the axis of the oblate spheroid is aligned with gravity, and is computed on a triangular grid (the influence of the grid characteristics is discussed by [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]) via a Newton algorithm [START_REF] Sipp | Global stability of base and mean-flows: a general approach and its applications to cylinder and open cavity flows[END_REF]. By introducing the cylindrical coordinate system (r, φ, x) such that y = r cos φ, z = r sin φ, the perturbation is expanded in normal mode form as q = [q f (r, x)e imϕ , qb ]e λt . By linearizing the governing equations around the base state, the whole problem reduces to a generalized eigenproblem of the form A m q = λB m q (with q = [q f , qb ]) for each value of the azimuthal wavenumber m. Finally, by making use of the SLEPc library, the eigenvalues λ = λ r + iλ i and the associated global modes q are obtained for various sets of the control parameters χ and Re, the bubble-to-fluid density ratio being set to ρ = 10 -3 , which is typical of the air-water system under standard conditions. The results discussed below have been made dimensionless by normalizing lengths with d, velocities with U 0 , time with d/U 0 and gravity with U 2 0 /d.

Results and discussion

3.1. Growth rates and frequencies as a function of Re and χ The most unstable modes all correspond to |m| = 1, as expected from the various available studies with axisymmetric bodies, i.e. fixed spheres and disks [START_REF] Natarajan | The instability of the steady flow past spheres and disks[END_REF], fixed bubbles [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF] or freely moving disks [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF]. These modes, with growth rate λ r and reduced frequency St = λ i /(2π), have been computed over the range χ ∈ [2.1, 2.5] and Re 3 × 10 3 . The upper value of χ was chosen after an attempt with χ = 2.7 showed that it was difficult to accurately determine the near-zero growth rate at the highest values of Re involved (≈5 × 10 3 in that case), due to difficulties in properly capturing the structure of the very thin boundary layer.

The curves λ r (Re) and St(Re) obtained by continuously tracking the most significant eigenvalues as a function of Re for various values of χ are gathered in figure 1(a-b) and (c-d), respectively. The growth rate curves evidence a destabilization-restabilization behaviour of the system, with λ r reaching a maximum at a value of Re that slightly increases with χ and ranges from 400 to 550 in the range of χ explored. Therefore, contrary to the case of buoyancy-driven rigid The dashed lines are associated with modes oscillating at a high quasi-constant frequency (HF modes); the thick solid lines correspond to stationary (S) or low-frequency (LF) modes or to combinations of both. The thin solid lines in (a) and (b) represent a second type of steady mode which at some point merges with the main S mode and gives birth to an LF mode, or sometimes (for χ = 2.25 -2.3 and 200 < Re < 300) to an HF mode.

bodies, the path of a spheroidal bubble is unstable only within a finite range of Re whose span increases with χ. The fact that the vorticity resulting from the shear-free condition reaches a finite value when Re goes to infinity instead of increasing as Re 1/2 when the fluid obeys a no-slip condition is at the root of this specific behaviour, as discussed in detail by [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF]. Figure 1(a) indicates that the growth rate becomes negative whatever the value of Re for χ χ C ≈ 2.15, a critical oblateness slightly less than the aspect ratio χ W ≈ 2.21 corresponding to the threshold of wake instability in the fixed-bubble configuration [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF].

Analysis of the reduced frequencies in figure 1(c-d) allows us to identify three main kinds of unstable modes.

(i) First, the most prominent feature is the existence of an 'S' family corresponding to stationary modes which are found to be unstable as soon as χ 2.215. These modes have the largest range of existence for most values of χ ; their amplification rate reaches a maximum in the range 400 Re 600.

(ii) Second, we can identify an 'HF' family of modes oscillating at a quasi-constant high frequency around the typical value St ∼ 0.1 (dashed lines). For χ 2.55, these modes are unstable in two distinct intervals, the first one being located around Re ≈ 100 and the second being centred at approximately Re ≈ 400. For χ 2.55 (not shown), these two intervals merge into a single one. It should be noted that for χ > 2.25, these modes are the first that become unstable as the Reynolds number increases. This is strikingly different from the picture observed in the fixed-bubble configuration, where the first bifurcation is always stationary [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF][START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF]. Figure 1(a,b) shows that the HF modes always restabilize at high Re.

(iii) Third, we observe an 'LF' family of modes oscillating at a quasi-constant but very low frequency, St ∼ 0.01. These modes are mostly encountered in a range of Reynolds number located just above that corresponding to the existence of the S modes, and their amplification rate is in most cases much smaller than that of the S and HF modes. It should be noted that, compared with the S modes, their range of existence also extends to lower values of χ (see e.g. the case χ = 2.15 in figure 1a).

A closer look at the results shows that, in addition to the aforementioned three main families, a second subdominant stationary mode (thin lines in figure 1a,b) is also present for χ 2.25 in a narrow interval of Re just below the upper threshold Re f of the S mode. This mode helps in understanding the onset of the LF mode, which can be interpreted as a result of the merging of two real eigenvalues into a pair of complex conjugate eigenvalues (see e.g. the curve corresponding to χ = 2.4). A similar process occurs close to the lower threshold of the S mode for χ 2.3, albeit in a much narrower range of χ (see e.g. the case χ = 2.217 in figure 1a). A more detailed view which allows one to fully appreciate how this merging/splitting process occurs is provided in appendix A where the unstable eigenvalue corresponding to χ = 2.217 is tracked in the complex plane.

Combination of the information contained in figure 1(a,b) reveals that all curves involving a stationary mode in some range of Re are actually made of a succession of the above families of modes. Hence, for a given aspect ratio, sudden frequency variations or 'jumps' are observed when the system switches from one mode to another. For instance, as Re increases, St goes from 0.1 to 0 and then to 0.01, corresponding to an HF-S-LF sequence at χ = 2.3, or from approximately 0.01 to 0 and then to 0.01 again for the LF-S-LF sequence at χ = 2.217. As we shall show below, these frequency changes are associated with dramatic changes in the spatial structure of the corresponding global modes.

Stability diagrams

By tracking the Re values such that λ r (Re) = 0 and the corresponding frequencies, we obtain the phase diagrams in the planes (χ, Re) and (χ, St) as drawn in figure 2. It is immediately seen that the HF, LF and S modes are unstable only for aspect ratios larger than the threshold values χ C HF ≃ 2.25, χ C LF ≃ 2.15 and χ S ≃ 2.215, respectively. Hence, path instability of spheroidal bubbles with a frozen shape is found to exist only for χ χ C = 2.15.

The grey lines in figure 2(a) correspond to the neutral stability curves of the wake alone, as obtained by [START_REF] Tchoufag | Linear stability and sensitivity of the flow past a fixed oblate spheroidal bubble[END_REF] for a fixed spheroidal bubble. Comparison of these curves with those of the fluid + bubble system reveals that the destabilization or restabilization thresholds of the wake are quite close to those of the coupled system, underlining the central role of the wake in the path instability mechanism. However, crucial differences may be noticed. First, the primary instability as Re increases is found to be of HF type for χ > 2.25, at odds with the fixed-bubble analysis which predicts the first bifurcation to be stationary for all χ 2.21. Second, instability is found to occur in the form of an LF mode for 2.15 χ 2.215, a range where the wake of a fixed bubble is stable whatever the value of Re. Hence, it is only within the narrow range χ ∈ [2.215, 2.25] that the first instability of the coupled system is well predicted by examining the stability of the sole wake. Clearly, the (χ, Re) and (St, Re) phase diagrams are much more complex for the frequency modes while the thick red dashed line without symbols corresponds to stationary modes. The thin red dashed line shows the location where the merges of (or splits into) stationary modes occur. The green crosses mark the locations of points P 1 -P 4 discussed in § 3.3.

coupled system. In figure 2(a), the critical curve corresponding to the stationary mode and that corresponding to the HF (respectively LF) mode are observed to cross at (χ, Re) ≃ (2.53, 151) (respectively ≃ (2.23, 290)), yielding codimension-two points of steady-Hopf/HF and steady-Hopf/LF type, respectively. Moving along the marginal curve of the primary instability, these codimension-two points are loci where a sudden frequency jump identified by an arrow in figure 2(b) takes place.

Global modes and comparison with DNS results

Figure 3 displays the spatial structure of a few critical global modes corresponding to points P 1 to P 4 in figure 2(a), their magnitude being normalized in such a way that the bubble inclination is unity. The above points were selected to allow comparisons with the DNS results of MM02.

Figure 3(a,b) reveals the wake structure of the LF mode that takes place at P 4 (χ = 2.25, Re = 10 3 ). The corresponding axial vorticity keeps a constant sign over a long distance (typically ten bubble diameters), which makes the near wake look similar to that associated with a strictly stationary mode (compare with panels e and f ). However, [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF] showed that the real and imaginary parts of an unsteady mode actually correspond to two different instants of time in its dynamics. Hence, comparison of panels (a) and (b) proves that the axial velocity and vorticity disturbances change sign as time evolves, although over a long time period. Such dynamics are reminiscent of the DNS observations of MM02 and [START_REF] Mougin | Wake-induced forces and torques on a zigzagging/spiralling bubble[END_REF], who found that in the zigzagging (ZZ) regime, high-Reynolds-number bubbles are followed by a pair of counter-rotating trailing vortices that change sign every half-period of the zigzag (see also the experimental observations in ultrapure water reported by [START_REF] De | Notes on the path and wake of a gas bubble rising in pure water[END_REF][START_REF] Van Wijngaarden | Bubble velocities induced by trailing vortices behind neighbours[END_REF][START_REF] Veldhuis | Shape oscillations on bubbles rising in clean and in tap water[END_REF]. To go beyond this qualitative resemblance, we may compare the frequency of the present mode with that found in the DNS run of MM02 corresponding to χ = 2.25, Ar = 390. In this run, the straight vertical path eventually switches to a ZZ path (figure 4a).

Starting from rest, the Reynolds number Re(t) (based on the norm of U) first increases rapidly until it reaches a maximum Re(t) = Re max close to but slightly less than that corresponding to the equilibrium value in the base state, Re 0 ≈ 1275. Then, Re(t) somewhat decreases when the instability develops (owing to the supplementary drag resulting from the occurrence of the trailing vortices) and eventually oscillates during the ZZ path. Since the bifurcation and the ZZ path both take place at an O(10 3 ) Reynolds number, the comparison with the linear mode at P 4 is relevant. However, it must be kept in mind that this mode was obtained assuming a strictly steady axisymmetric base flow, whereas the flow in the DNS evolves in time at previous stages and is three-dimensional far before the conditions corresponding to P 4 are reached, since figure 2(a) indicates that the stationary mode becomes unstable at Re ≈ 270 for this aspect ratio. Defining the reduced frequency normalized by the gravitational velocity, St = St U 0 /U g , we find it to be St ≃ 0.042 for the linear LF mode at its restabilization threshold (Re ≈ 1180), which compares reasonably well with the DNS finding, St DNS ≃ 0.053, corresponding to figure 4(a). This agreement, although not perfect for the reasons outlined above, supports the idea that, in this range of χ, the zigzagging motion of bubbles with O(10 3 ) Reynolds numbers is associated with an LF mode, as figure 2(a) predicts.

Let now turn to the HF mode at an O(10 2 ) Reynolds number shown in figure 3(c,d). Comparison of the two panels makes it clear that the spatial arrangement of the axial velocity evolves periodically from a 'sign-preserving' type (all 'rolls' in the wake keep the same sign) to a 'sign-alternating' type (two successive rolls have opposite signs), a switch that [START_REF] Tchoufag | Global linear stability analysis of the wake and path of buoyancy-driven disks and thin cylinders[END_REF] showed to be indicative of a strong coupling between the body and its wake at the corresponding bifurcation. It is then no surprise that the fixed-bubble approach could not predict the bifurcation that takes place at P 1 . As figure 2(a) indicates, point P 1 corresponds to the first bifurcation (with increasing Re) for the aspect ratio χ = 2.5. This bifurcation is of Hopf type, and the threshold, once expressed using the gravitational velocity, is found to be Ar ≃ 90.5. These findings are in good agreement with those of MM02, whose DNS predicted the threshold of the zigzagging path to take place at Ar DNS ≃ 86 when χ = 2.5. Moreover, the frequency St is approximately 0.095 for the HF mode, which also compares well with the DNS prediction St DNS ≃ 0.092.

As shown in figure 3(e, f ), the wake associated with the stationary modes at P 2 (Re = 157.3) and P 3 (Re = 1050), both for χ = 2.5, is dominated by the presence of two counter-rotating trailing vortices, a wake structure known to result in a steady non-zero lift force. At P 2 , the near wake also involves two small secondary vortices of opposite sign directly connected to the bubble, a feature that does not exist at P 3 . This difference may be understood by noting that, unlike the case of solid bodies, the standing eddy at the rear of a bubble shrinks as Re becomes very large [START_REF] Blanco | The structure of the axisymmetric high-Reynolds number flow around an ellipsoidal bubble of fixed shape[END_REF][START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF], so that the streamline distribution in the base flow gets closer to that corresponding to the potential flow past an oblate spheroid. This change is also responsible for the narrowing of the wake as Re increases, which tends to 'glue' the trailing vortices to the symmetry axis, as figure 3(f ) reveals. Again, a comparison with the DNS predictions can be achieved in the latter case. With χ = 2.5 and Ar = 390, MM02 found the bubble Reynolds number to reach a maximum Re max ≈ 1035 before the path departed from vertical. Then, this path was observed to experience a marked lateral drift until a periodic ZZ motion eventually set in (see figure 4b). As the Reynolds number at P 3 is close to Re max , it is again reasonable to consider that the linear mode at P 3 represents a good approximation of the mode that first becomes unstable in the DNS. Since we found the former to be stationary, LSA allows us to conclude that the first bifurcation experienced by the system when Re(t) becomes of O(Re max ) yields a steady oblique (SO) path, a prediction in agreement with figure 4(b).

Summary and final comments

In this work, we investigated the global linear stability of the path of an oblate spheroidal bubble rising in a Newtonian fluid otherwise at rest. The computed eigenvalues and modes of the coupled fluid + bubble system revealed a much more complex behaviour than in the case where the stability of the sole wake is considered. In particular, we found that the first instability of the system as the Reynolds number increases may arise through a Hopf bifurcation associated with a low-or high-frequency mode, depending on the aspect ratio, while it always occurs through a stationary bifurcation when the bubble is set fixed. Comparisons with DNS results performed for several sets of χ and Ar showed that the nature and frequency (if any) of the mode that first becomes unstable are well captured by the present LSA approach. However, several open issues remain.

First, as we have already mentioned, comparison of LSA results with experimental and DNS observations is not straightforward because bubbles observed in 'real life' are generally released from rest, which makes their rise Reynolds number highly time-dependent, especially during the initial transient. This raises a question about the type of instability that actually develops when the aspect ratio is such that the physical system successively explores several unstable regions as Re(t) evolves. This is, for instance, the case for the bubble with χ = 2.25, Ar = 390 we discussed in the previous section, which according to figure 2(a) crosses the region where the S mode is unstable when it accelerates and its Reynolds number is in the range 270 Re(t) 630, before it enters the region where the LF mode is unstable when it gets closer to its steady state. In that case, the comparison we performed with DNS data revealed that the final non-vertical path exhibits characteristics consistent with those of the linear LF mode; the only discernible trace of the S mode is probably the slight difference we observed in the values of St, which may presumably be attributed to the three-dimensionality of the flow that has already set in before the LF mode becomes unstable. However, the influence of modes that are unstable during some stage of the initial transient is probably more severe in certain cases, such as for instance with bubbles such that 2.3 χ 2.35 and Ar values yielding final Reynolds numbers slightly beyond the destabilization threshold of the LF mode (the dotted line with blue open circles in figure 2a). Here, the system successively covers ranges of Re(t) where the HF, S and LF modes are unstable. In particular, it spends a significant amount of time in the intermediate region where the stationary mode is unstable. Since the growth rate of this mode is much larger than that of the other two according to figure 1(a), one may guess that the late evolution of the real system carries some footprint of this S mode. Unfortunately, MM02 did not perform any run with such (χ, Ar) sets and we cannot check this conjecture here.

A similar difficulty is encountered after the bubble has reached its maximum Reynolds number and may be illustrated with the example of figure 4(b). Here, Re(t) reaches its maximum value, Re max ≃ 1035, before the SO path sets in. During this second stage, Re(t) sharply decreases due on the one hand to the supplementary drag induced by the trailing vortices and on the other hand to the inclination of the path which lowers the driving component of the buoyancy force. At the end of the SO stage, Re(t) is approximately 800 which, according to figure 2(a), is close to the upper threshold of the HF mode. Therefore, looking at the path displayed in figure 4(b), it is tempting to conclude that the ZZ stage that succeeds the oblique drift corresponds to this HF mode. Actually this is not so, since the frequency of the latter yields St ≃ 0.175 according to figure 2(b), while the frequency determined in the DNS is St DNS ≃ 0.064. A weakly nonlinear approach is required to address the above issues and determine how the primary modes that are successively encountered as Re(t) evolves interact and what the essential characteristics of the final stage of the system are. We are currently developing such an approach.

Comparison of the present predictions with known results for real bubbles raises another issue. According to the experimental findings of [START_REF] Zenit | Path instability of rising spheroidal air bubbles: a shape-controlled process[END_REF] obtained with various silicone oils, the path may become unstable for aspect ratios typically 10 % smaller than the critical value predicted here. Moreover, for O(10 2 ) Reynolds numbers, [START_REF] Zenit | Measurements of the streamwise vorticity in the wake of an oscillating bubble[END_REF] showed that the wake structure is very similar to that corresponding to the LF mode displayed in figure 3(a,b) (with, in addition, the pair of secondary vortices found in the same range of Re in the S mode of figure 3c). We would guess that both features result from the slight fore-aft asymmetry of real nearly spheroidal bubbles: not only might this asymmetry shift the marginal curve toward smaller χ, as suggested by the numerical study of [START_REF] Cano-Lozano | Wake instability of a fixed axisymmetric bubble of realistic shape[END_REF], but it also seems to deform the stability diagram in such a way that the threshold corresponding to the destabilization of the LF mode is shifted toward much smaller Reynolds numbers than for perfectly spheroidal bubbles and may, over a significant range of χ, be the first that becomes unstable as Re increases. We are currently working on this aspect by introducing more realistic bubble shapes in our LSA code.

Finally, the existence of a crescent-shaped region in figure 2(a) within which the S mode is unstable provided the aspect ratio is larger than 2.215 raises an interesting puzzle. Indeed, this prediction suggests that bubbles whose final Reynolds number falls in this region should eventually follow an SO path. This is especially clear within the narrow range 2.23 χ 2.25 where the S mode is associated with the primary bifurcation, i.e. no previous bifurcation can alter the final behaviour of the system. The point is that we are not aware of any experimental or computational observation that supports the existence of such an inclined path, except during short transients as in figure 4(b) (again MM02 did not consider any (χ, Ar) set whose final state should have fallen in this region). Only zigzagging or spiralling paths have been observed to subsist over long periods of time, suggesting that the SO path is not a solution of the full nonlinear problem. This is at odds with the case of solid spheres and disks for which such a path has been consistently reported, in both experiments and computations, and was predicted to be a stable non-trivial solution of the fluid + body problem though a weakly nonlinear analysis by [START_REF] Fabre | The steady oblique path of buoyancy-driven disks and spheres[END_REF]. Crucial in this analysis is the fact that the torque experienced by the body vanishes for a small non-zero inclination of the path. It might be that the shear-free rather than no-slip condition that holds at a bubble surface changes the nature of the bifurcation that yields the SO path: it was found to be supercritical for spheres and disks (except when the latter are extremely thin) and could be subcritical for bubbles. We are currently repeating the analysis developed by [START_REF] Fabre | The steady oblique path of buoyancy-driven disks and spheres[END_REF] with a bubble to check this possibility.

FIGURE 1 .

 1 FIGURE 1. Variations of the growth rate (a,b) and reduced frequency (c,d) as a function of Re for several bubble aspect ratios: (a) and (c) χ = 2.15-2.25; (b) and (d) χ = 2.3-2.5.The dashed lines are associated with modes oscillating at a high quasi-constant frequency (HF modes); the thick solid lines correspond to stationary (S) or low-frequency (LF) modes or to combinations of both. The thin solid lines in (a) and (b) represent a second type of steady mode which at some point merges with the main S mode and gives birth to an LF mode, or sometimes (for χ = 2.25 -2.3 and 200 < Re < 300) to an HF mode.

FIGURE 3 .

 3 FIGURE 3. Several global modes (m = 1) of the coupled system at points P 1 -P 4 in figure 2(a). The upper (lower) half of each subfigure displays the isovalues of the axial velocity (vorticity) of the mode; gravity acts from left to right. (a,b) The LF mode at P 4 (χ = 2.25, Re = 10 3 ). (c,d) The HF mode at P 1 (χ = 2.5, Re = 90.6) (in both cases, (a,c,e) and (b,d,f ) display the real and imaginary parts, respectively). (e,f ) The stationary modes at P 2 (χ = 2.5, Re = 157.3) and P 3 (χ = 2.5, Re = 1050), respectively; for these two modes, the axial velocity (vorticity) is purely real (imaginary).

FIGURE 4 .

 4 FIGURE 4. Path of two oblate bubbles corresponding to Ar = 390 with aspect ratio (a) χ = 2.25 and (b) χ = 2.5. In both cases, the whole system is initially at rest (reprinted from Mougin 2002).

FIGURE 5 .

 5 FIGURE 5. Complex path of the unstable eigenvalue found with a bubble of aspect ratio χ = 2.217. The left grey part of the plane indicates the region where the corresponding eigenmode is stable.

† Email address for correspondence: magnau@imft.fr

Acknowledgement

This work was supported by the Agence Nationale de la Recherche under grant ANR-09-BLAN-0132 OBLIC.

Appendix. Merging/splitting of eigenvalues in the complex plane

To facilitate the understanding of the merging/splitting behaviour exhibited by several eigenvalues in figure 1, figure 5 shows how the eigenvalue with the largest real part found for χ = 2.217 evolves in the complex plane. We selected this example as a prototype of eigenvalues whose imaginary part vanishes in some range of Re, then yielding a stationary mode. The arrows in figure 5 indicate the direction of increasing Re, from 1 to 4. Following arrows 1, we first observe a pair of complex conjugate eigenvalues with a negative real part. This real part then changes sign, leading to an instability corresponding to an LF mode, until the two eigenvalues merge at point (0.007, 0). There they give birth to a pair of real eigenvalues. Following the arrow 2 pointing to the right, one of them (which corresponds to the main stationary mode in figure 1) is seen to increase until it reaches the point (0.048, 0). The second real eigenvalue (corresponding to the subdominant stationary mode in figure 1) decreases along the second arrow 2 and changes sign (i.e. the corresponding mode restabilizes) until it reaches the point (-0.025, 0). Then, following arrows 3, these two eigenvalues go through their previous paths in the reverse direction. They eventually merge at point (0.014, 0) where they split into a new pair of complex conjugate eigenvalues (corresponding to a second LF mode). Following arrows 4, the magnitude of their imaginary part sharply increases without much change in the real part, until it saturates at a value of ±0.017. Then, the real part sharply decreases and eventually becomes negative, leading to restabilization.


Olivier Baldellon

Jean-Charles Fabre

Matthieu Roy

Distributed Monitoring of Temporal System Properties using Petri Nets

Keywords: Distributed Monitoring, Online Verification, Petri nets

Supervising a system in operation allows to detect a violation of system specification or temporal properties, and is the first step required by any reconfiguration mechanism.

In this work, we focus on run-time verification of temporal system properties in distributed and realtime systems. Based on a description of a property that includes events and temporal constraints, expressed as an arc timed Petri net, we automatically derive a monitoring system responsible for checking this property. The proposed approach enables the distributed verification of system properties.

Our contribution is twofold. On the theoretical side, we introduce a slight modification of the semantics of Petri nets to be able to execute it in partial executions and noisy observation environments. On the practical side, we show how to use this formal framework to provide a distributed and efficient monitoring system, and describe its current implementation.

Introduction 1

In the different means to provide dependability guarantees in systems, supervising [MR10], [JRR94], [ZSLL09],

[CJ05], or monitoring systems and applications states is a requisite to detect a possible violation of system specification and envisage a recovery action . In a first section, we introduce a slight modification of the semantics of Petri nets to be able to execute it in partial executions and noisy observation environments. Then, we show how to use this formal framework to provide a distributed and efficient monitoring system, and describe its current implementation.

Monitor Alarm

(e1, τ1) (e3, τ3) Figure 1. Conceptual view
Our monitor tool is a distributed service that executes a model of the system expressing temporal and behavioral properties on reception of timed events from the system.

t 1 [0, ∞] t 2 [3, 6] t 3 [0, 5] Figure 2. A simple Petri net I. Our Approach A. Petri Nets For System Properties An arc-timed Petri net is a tuple (P , T , • (•), (•) • , M 0 , I)
, where P is a finite set of places, T is a finite set of transitions, M 0 is the initial marking (a function that associates to each place a set of tokens), • (•) and (•) • are the backward and forward incidence functions that associate to each transition a set of places, and I a function that associates to each arc between a place and a transition a time interval. We say that there is an arc between a place p and a transition t if p ∈ • t.

We assume there is a bijection between transitions and events generated by the system. Thus a timed event (e i , τ i ) can unambiguously be written (t i , τ i ) where t i is the transition associated to e i .

B. Petri Net Execution

Properties on system events that need to be monitored are expressed using the Petri net formalism. Our approach consists in executing a Petri net on-the-fly to detect failures. The monitoring tool takes as an input a sequence of events, i.e., a sequence of couples (t, τ ) where t represents the transition associated to an event and τ the date of the event. The monitoring tool executes the model using such events sequence, and possibly raises an alarm when it detects an incorrect behavior.

A simple execution: Let us take as an example the simple Petri net depicted in Figure 2: if the monitor receives, the sequence of events (t 1 , 10); (t 2 , 15); (t 3 , 21) in this order, it will first fire t 1 , storing that this firing was done at time 10. When the event corresponding to t 2 is received, the transition t 2 is fired because 15 -10 ∈ [3,6]. The reception of (t 3 , 21) will raise an error because the event occurred too late; the token stayed in the place of • t 3 six units of time (τ 3 -τ 2 = 21 -15 = 6), which is out of the specified interval [0, 5]. However, the transition t 3 will still be fired to enable the execution of the Petri net model to continue, and an error will be reported.

In this scenario, notice that the minimal amount of information to detect a failure is (t 2 , 15); (t 3 , 21): whatever the occurrence time of t 1 , t 3 happened too late with respect to t 2 . If the monitor only receives (t 2 , 15); (t 3 , 21) but (t 1 , 10) is not yet received, the usual semantics of Petri nets forbids to fire transitions t 2 and t 3 ; in other words, to fire t 2 and consequently t 3 , the transition t 1 needs to have been already fired. This crucial issue can be solved using negative tokens.

An execution with negative tokens: Our approach consists in firing a transition t as soon as it is received and anticipating the removal of tokens in places of • t by negative tokens in these places. Figure 3 shows the result of firing t 2 , with a negative token represented as a black circle, a positive token as a black disk.

p 0 t 1 [0, ∞] p 1 t 2 [3, 6] p 2 t 3 [0, 5] p 3 Figure 3. After the firing of t 2
The firing of t 3 will add one negative token in p 2 and one positive token in p 3 . The removal of a token always depends on the presence of its negative counterpart. To know how long the positive token stayed in the place p 2 , we need to compare the date of the event that created the negative token with the date of the event that created the positive one. If the difference is not in the time interval I(p 2 , t 3 ) = [0, 5], then an error is raised. In all cases, both tokens are removed.

In summary, the main difference between our approach and classical semantics is that, in the classical one, a transition is allowed to be fired only if all required tokens are present. In our approach, transitions are always fired speculatively, and the fireability property is checked a posteriori.

II. Protocol & Implementation

A. The protocol

Due to lack of space, the formal description of the protocol is omitted. It implements a simulation of the Petri net by a set of thread, one for each place and one for each transition. The protocol is then a direct transcription of the execution semantics described in I-B with negative tokens. The protocol assigns a monitoring thread to each place and each transition of the Petri net. The communication graph of the Petri net of Figure 2 is presented in Figure 4: when the system generates an event, this event is sent to the thread associated to the corresponding transition. In return, transition threads send tokens to threads associated to places when they receive events.

When a place thread receives a positive or negative token, it starts a timer. We make the assumption that the network is synchronous (bounded delay) and thus when a token is received, the place can compute the maximal delay to wait for reception of the negative version of this already received token.

(t 1 , τ 1 ) (t 2 , τ 2 ) (t 3 , τ 3 ) Figure 4. Communications Graph

B. Implementation

We have developed a proof-of-concept of our approach in the Erlang language. The full code of our implementation, called Minotor, is available on the author website2 . The implementation is very scalable, since we were able to run series of tests on Petri nets of size up to 2 20 (more than one million) transitions and places.

Conclusion

Differently of static analysis, a run-time verification system receives system events on the fly, and possibly out-of-order. To cope with this problem, we introduce the notion of virtual, or signed token, i.e. we execute every transition associated to an event as soon as this event is caught by the monitoring system, no matter the current state of the Petri net. The actual verification of timing constraints and events ordering is performed a posteriori by checking respective dates of signed tokens.

The decoupling of transitions firing and timing constraints verification allows us to completely distribute the verification, as we show in the article. Our strategy is to associate to each transition, and each place in the Petri net a conceptual thread, in charge of executing an atomic subpart of the Petri net and verifying locally that timing constraints are valid. References [JRR94] F. Jahanian, R. Rajkumar, and S.C.V. Raju. Runtime monitoring of timing constraints in distributed real-time systems. Real-Time Systems, 7(3):247-273, 1994.

[MR10] P. Meredith and G. Roşu. Runtime verification with the rv system. In Proceedings of the First international conference on Runtime verification, pages 136-152. Springer-Verlag, 2010.

[ZSLL09] W.

  and C. Jard. Time supervision of concurrent systems using symbolic unfoldings of time petri nets. Formal Modeling and Analysis of TimedSystems, pages 196-210, 2005.

This work has been partially supported by ANR, the French Science Foundation, under contract number ANR-BLAN-SIMI10-LS-100618-6-01.

http://www.olivier.baldellon.eu/documents/minotor-srds.tar.gz


Félicie Korostelev
email: felicie.korostelev@upmc.fr

Clémence Basuyau

Sylvie Leroy

Christel Tiberi

Abdulhakim Ahmed

Graham W Stuart

Derek Keir

Frédérique Rolandone

Ismail Al Ganad

Khaled Khanbari

Lapo Boschi

Crustal and upper mantle structure beneath south-western margin of the Arabian Peninsula from teleseismic tomography

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The Afar triple junction, where the Red Sea, East African and Gulf of Aden rifts intersect (Fig. 1), is a key region to understand how continental breakup occurred under the influence of a plume and abnormally hot mantle [e.g. Bastow et al., 2 0 1 1 The YOCMAL project (YOung Conjugate MArgins Laboratory) deployed 23 broadband stations in a network running from the Red Sea margin to Aden city, passing through the Yemeni highlands and Sana'a city (Fig. 2). Using classical teleseismic tomography [Aki et al.,1 9 7 7 ]o nr e c o r d i n g sf r o mt h e s es t a t i o n st o g e t h e rw i t hap e r m a n e n tG e o f o ns t a t i o n (DAMY), we image the relative velocity variations of P-waves in the crust and upper mantle down to 300 km depth. We thus: [START_REF] Jordan | A new method for the 3-D joint inversion of teleseismic delaytimes and Bouguer gravity data with application to the French Massif Central[END_REF] characterize the lithospheric structure of the rifted margins of the Gulf of Aden and Red Sea system of western Yemen; (2) locate

DRAFT

A p r i l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT X-5

the presence of asthenospheric upwellings in the region and their interaction with the lithosphere.

Geodynamic setting

The Red Sea and Gulf of Aden rifts are connected to the East African rift at the Afar [START_REF] Leroy | Recent off-axis volcanism in the eastern Gulf of Aden: Implications for plumeridge interaction[END_REF] (Fig. 1). The increased magmatism caused by extension above a plume was first thought to stop at the Shukra al Sheik fracture zone [e.g. [START_REF] Hébert | Lithospheric structure of the nascent west Gulf of Aden spreading center inferred from gravity data[END_REF], but these recent results indicate that the limit could be further east (Fig. 1), as proposed by Leroy et al. [2010a].

The Gulf of Aden is characterized by two stretched continental margins systems: non volcanic margins in the east and volcanic margins in the west near the Afar triple junction.

The volcanic margins are associated with syn-rift SDRs up to 5 km thick [START_REF] Tard | The volcanic passive margin in the vicinity of Aden, Yemen[END_REF][START_REF] Leroy | From rifting to oceanic spreading in the Gulf of Aden: a synthesis[END_REF]. SDRsarealsosparselyfoundintheeastespeciallyattheoceancontinent transition [Autin et al., through the Yemeni highlands (Fig. 2). The sensors deployed were Guralp 40T (sampling rate 50 sps, 30 s natural period), 6TD (sampling rate 40 sps, 30 s natural period) and ESPD (sampling rate 40 sps, 60 s natural period). This network configuration allows imaging structures with a high resolution, and down to 300 km depth (the surface extent of our network). We selected 200 teleseismic events with clear first P-wave arrivals from earthquakes with magnitude ≥5.5. Among them, 142 events arrived as P, 35 events as PP and 23 events as PKP (Fig. 3), and were picked on consistent peaks in the first cycle.

We used 2456 delay times calculated relative to the IASP91 Earth model [Kennett and Engdahl ,1 9 9 1 ] . F o re a c hr e s i d u a l ,ap i c k i n ge r r o rw a sa s s i g n e dw i t h i nt h er a n g e±0.05

to ±0.2 s. As the events registered come mostly from east / northeast and have a range of epicentral distances, the rays are crossing, and we expect our best resolution in this direction.

Method

Our delay time residuals were inverted using the ACH teleseismic tomography method first used by [START_REF] Aki | Three-dimensional seismic velocity anomalies in the crust and upper-mantle under the USGS, California seismic array[END_REF], then developed by [START_REF] Weiland | Nonlinear teleseismic tomography at Long Valley caldera, using three-dimensional minimum travel time ray tracing[END_REF], [START_REF] Zeyen | Joint Inversion of Teleseismic Delay Times and Gravity Anomaly Data for Regional Structures. Theory and Synthetic Examples[END_REF] and [START_REF] Jordan | A new method for the 3-D joint inversion of teleseismic delaytimes and Bouguer gravity data with application to the French Massif Central[END_REF] [START_REF] Ahmed | Crustal structure of the rifted volcanic margins and uplifted plateau of Western Yemen from receiver function analysis[END_REF] to compute a time correction for each station (Fig. 6). This correction accounts for lateral variations in the crustal thickness, which are difficult to constrain by our teleseismic data alone; seismic velocity within the crust is treated as a free parameter in our inversion. The correction for each station is computed on the absolute residuals and is in the range -0.611 s (for the station located on the thinnest crust) to 0.160 s (for the station located above the thicker crust). This allows a reduction in the propagation of crustal structure into the deeper layers of the velocity model.

The smoothing factor, which limits the short wavelength velocity variations, and the ini- l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT X-9

DRAFT A p r i
tial standard deviation associated with each node of the initial model for each of the nine iterations were chosen after a series of tests and are respectively 0.001 and 0.007 km.s -1 .

The smoothing value chosen is the same as [START_REF] Tiberi | Asthenospheric imprints on the lithosphere in Central Mongolia and Southern Siberia from a joint inversion of gravity and seismology (MOBAL experiment)[END_REF][START_REF] Leroy | Recent off-axis volcanism in the eastern Gulf of Aden: Implications for plumeridge interaction[END_REF].

Tests showed that the results are not significantly different for a standard deviation b etween 0.005 and 0.01 km.s -1 .W e t e s t e d t h e s e i n v e r s i o n p a r a m e t e r s i n o r d e r t o g e t a stable model, choosing parameters that decrease the root mean square of residuals (RMS)

through the inversion's nine iterations. The overall decrease of the RMS is 55% in the final model.

Results

Checkerboard test

We use a checkerb oard test to assess the resolution of our inversion. The synthetic checkerboard model consisted of rectangular velocity anomalies of +5% and -5% velocity variations at depths of 70, 200 and 300km depth (Fig. 7a). The size of the anomalies increases around the edges as the method requires the nodes to be further apart at the edges of the model. We use the same inversion parameters (smoothing, standard deviation) as in our actual inversion.

The results of the checkerboard test (Fig. 7b) show that at 70 km depth, due to the concentration of crossing rays under the stations, we have the best resolution, with ∼ 20 % recovered amplitude. At 200 km depth, the crossing rays cover a larger area, so the anomalies are well retrieved. At 300 km depth, the eastern retrieved anomalies are much better constrained east of 44 • E, due to the concentration of ray paths coming from the east. The Fig. 7b shows the piercing points under our network at 300 km depth, at our l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT maximum investigation depth. Below 300 km depth, the rays are too dispersed to give good resolution.

DRAFT A p r i

Crustal-scale structures

The thickness of the crust in our study area is estimated from receiver function analysis to be between 14 km at the coast and 35 km inland [Ahmed et al., 2 0 1 3 ] . T h ec r u s t is represented by the first two layers of our model, and is characterized by the highest anomaly contrasts in P-wave velocity, which have a range of ±4% (Fig. 4 and5). These anomalies are related to geological structures observed at the surface. At 30 km depth, the strong low velocity anomaly (-4%) beneath the center of our network (stations MAWI to DALA) correlates with the high topography of the plateau (>2000 m above sea level). We interpret this pattern as due to the thicker crust (>30 km) beneath the Yemeni highlands The high velocity anomalies under the stations ANID and SUGH could be due to a thin crust as they are in the Red Sea coastal area, or to Tertiary granitic intrusions of Jabal Hufash and Jabal Bura respectively [Geoffroy et al., 2002], see Fig. 2. The high velocity anomaly corresponding to Jabal Hufash is imaged down to 45 km depth, which is slightly

DRAFT A p r i l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT X -11
deeper than that of Jabal Bura (30 km depth). This can be explained by significant smearing due to the higher amplitude of the Jabal Hufash anomaly. There are lower velocities under the stations UAYA and ZUWA, on the Tihama plain (Fig. 5b), probably due to ∼4000 m of low velocity sediments [El-Anbaawy et al.,1 9 9 2 ;Davison et al.,1 9 9 4 ] .

Upper mantle structures

The resulting upper-mantle P-wave anomalies are in the range of ±2%. The most striking pattern is a low velocity anomaly located under the Yemeni highlands at 45 km depth, apparently dipping northeastward down to 300 km. It reaches its maximum amplitude at 70 km depth (east of the DAMY station) beneath the volcanic field of Dhamar (Fig. 2),

where there are two active volcanoes [START_REF] Manetti | Magmatism of the eastern Red Sea margin in the northern part of Yemen from Oligocene to present[END_REF]. The northern part of this anomaly is located, at 70 km depth, under the volcanic field of Sanaa, which is also still active [START_REF] Manetti | Magmatism of the eastern Red Sea margin in the northern part of Yemen from Oligocene to present[END_REF], see Fig. 5.

There is a second low-velocity anomaly located under the southwestern corner of Yemen and the stations MOKA and HAKI (Fig. 5c andd). This low velocity anomaly is nearly vertical and is recognized from the surface to 300 km. It is located just beneath the volcanic area of Jabal An Nar (Fig. 5d), which was active during late Miocene, around 10 Ma, [Manetti et al.,1 9 9 1 ] .

Even if we corrected the residuals from the crustal portion, our results are likely to include effects from Moho variations. This is because the corrections are based on receiver functions which display a strong trade-off between Moho depth and crustal velocity model.

In addition, we took a 1D velocity model to be consistent with the receiver functions study of [Ahmed et al., corrections than we have done. To estimate the resolution of our models in both the crust and the mantle, we proceed to synthetic tests. We test whether the low velocity anomaly beneath the high plateau is dipping towards the northeast because of smearing along ray paths (Fig. 3) in the next section. We investigate by means of synthetic tests whether:

(1) the velocity variations observed in our resulting model are smearing downwards into the mantle, (2) the low velocity anomalies under the southwestern corner of Yemen and under the high plateaus are artefacts, and whether we can determine at what depth they are located, (3) the dipping low-velocity anomaly is related to the presence of partial melt or not.

6. Synthetic tests and presence of melt

Propagation of crustal signal

We test the smearing with depth of a crustal anomaly (0-30 km) without a crustal correction by introducing a -5% anomaly around stations in the Yemeni Highlands and a +5% anomaly under the two continental margins (Fig. 8a). The rest of the input model is laterally uniform. The results of the inversion of this synthetic model (Fig. 8b andc) show that the velocity variations are well recovered in location but with 40% lower amplitude.

Between 45 and 100 km depth, we can discern very low amplitude anomalies (<0.5%), corresponding to a small amount of vertical smearing. At 100 km depth, we have a -1% anomaly under the Yemeni highlands. Approximately half of this signal could be due to a smearing of the crust signal. No significant anomaly can be seen beneath 100 km depth. We conclude that crustal velocity anomalies do not propagate to deeper layers of our model and that there is an authentic low velocity anomaly in the upper mantle. This

DRAFT A p r i l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT X -13
test shows that the use of a crustal correction is relevant and is necessary to limit the extent of the propagation of crustal velocity structure into deeper layers of the model.

Low velocity anomalies beneath Southwestern Yemen

Using a synthetic input model, we simulate the resulting geometry of our final P wave model (see auxiliary material). We computed a series of tests (available in auxiliary material), and here we present the most relevant example. We first place a -5% anomaly under the three Yemeni volcanic fields of Sana'a, Dhamar and Marib from the base of the crust down to 200 km depth. The other -5% anomaly is placed from the base of the crust to the base of our model, under the southwestern corner of Yemen and southern Red Sea. Figure 9 shows the retrieved inversion image for a SW-NE profile (compare with the observed results in Fig. 9c). Although the anomaly amplitudes within the crust are not retrieved, the dipping anomaly is quite well retrieved in the synthetic output model, as well as the low velocity anomaly beneath Jabal An Nar volcanic field (Fig. 9). This tests shows that the dipping low velocity anomaly under the Yemeni volcanic fields can be explained by a 220x260 km mantle upwelling between the base of the crust and 200 km depth. Moreover, the low velocity anomaly beneath Jabal An Nar can be explained by a large zone of hotter mantle.

Presence of partial melt

Traveltime tomography gives the present state of the upp er mantle in terms of velo city variations but it precludes any direct interpretation concerning their origin. Indeed, strates that a purely thermal origin leads to a linear relationship between P and S residuals with a slope of 2.9. Gao et al. [2004] proposed that a slope higher than 2.9 between P and S residuals for the same events highlights the presence of partial melting and several authors used the relationship to explain the presence of negative velocity anomalies in We thus selected the events coming from northeast (Russia, Japan, China), and picked the S arrival on the transverse component for stations located above the Yemeni Highlands, so that the rays chosen are passing through the northeastward dipping low-velocity anomaly.

In Fig. 10, following Gao et al. At shallower depths of 0-30 km in our model, seismic anomalies are directly related to the geological units observed at the surface. The granitic intrusions of Jabal Hufash and Jabal Bura are associated with high velocity anomalies (up to +4%, see Fig. 2). The depth extent of the anomaly below the granitic intrusions of Jabal Hufash may indicate a deeper root, as hypothesized by Denèle et al. [2012]. It could also be explained by a stronger smearing effect due to the high amplitude of the Jabal Hufash anomaly.

Deep structure of the margins

Our synthetic tests show that deep anomalies cannot be explained by the smearing in depth of crustal anomalies. We interpret the low velocity anomaly (-2%) between ∼35 km and 300 km depth under the highest topography as abnormally hot mantle upwelling. This low velocity hotter mantle is located below the stations DAMY, RUSA and YSLE (Fig. 5), which are located on the thick Oligocene flood basalts and three more recent volcanic fields (15 Ma to present) volcanic areas, e.g. Dhamar and Sana'a [Davison et al.,1994;[START_REF] Pik | Timing of East African Rift development in southern Ethiopia: Implication for mantle plume activity and evolution of topography[END_REF][START_REF] Leroy | Contrasted styles of rifting in the eastern Gulf of Aden: A combined wide-angle, multichannel seismic, and heat flow survey[END_REF], see Fig. 2. Additionally, we infer presence of partial melt in the crust or uppermost mantle to be responsible for this low velocity anomaly (Fig. 10). A similar pattern of upper mantle off-axis upwelling has also been found in the southern Red Sea rift of Afar and explained by small diapiric upwellings l2 1 ,2 0 1 4 ,3 : 2 6 et al.,2003]. Moreo v er, [START_REF] White | Magmatism at rift zones: the generation of volcanic continental margins and flood basalts[END_REF] explained that an increase of about 150 • Cinthema n tlele a dstoadyna micuplif t,buttha tthea dditio n of dense igneous material under a thinned crust produces an immediate subsidence of DRAFT A p r i l2 1 ,2 0 1 4 ,3 : 2 6 p m DRAFT more than 2 km in order to maintain isostatic equilibrium. We observe a similar pattern on the Red Sea and Gulf of Aden margins, with high Yemeni plateau dynamically supported by a hot mantle (Fig. 5a, b andd), and a subsiding area, for example the Tihama Plain, underlained by seaward dipping reflectors, and ultra-mafic bodies accreted under the crust. This hypothesis should however be tested for Yemen by gravity and isostatic modeling.

DRAFT A p r i
The weak low velocity anomaly imaged under the southwestern corner of Yemen, beneath the MOKA station, is underneath the Miocene volcanic area of Jabal An Nar. Our synthetic test (Fig. 9) shows that this low velocity anomaly may be explained by a large zone of hotter mantle. This feature could be due to the Afar plume signal, located only ∼300 km away, in the Afar depression.

Conclusions

We p erformed an inversion of P-wave teleseismic data to image lithospheric structure beneath the SW of the Yemen, the southern Red Sea and western Gulf of Aden margins.

The crustal part of the model is dominated by a possible ultra-mafic underplating beneath the Red Sea and Gulf of Aden margins, a sudden thinning of the crust for this volcanic margin, and Tertiary granitic intrusions (Jabal Hufash and Jabal Bura). In the mantle, we image a low velocity anomaly in which we infer partial melting just below the highest topography, the thick Oligocene flood basalts and other off-axis volcanic regions (from 15 Ma to present). This low velocity anomaly could correspond to an abnormally hot mantle and could be responsible for dynamic topography and recent magmatism in western Yemen. Some hot material has also been inferred under the southwestern corner of Yemen and may be due to the Afar plume signal. l2 1 ,2 0 1 4 ,3 : 2 6 Benoit, M., A. Nyblade, J. VanDecar, and H. Gurrola (2003), Upper mantle P wave velocity structure and transition zone thickness beneath the Arabian Shield, Geophys.

DRAFT A p r i
Res. Lett, 30 (10), 1531.

Benoit, M., A. Nyblade, and J. VanDecar (2006), Upper mantle P-wave speed variations beneath Ethiopia and the origin of the Afar hotspot, Geology, 34 (5), 329-332. Bertrand, H., G. Chazot, J. Blichert-Toft, and S. Thoral (2003), Implications of widespread high-µ volcanism on the Arabian Plate for Afar mantle plume and litho- l2 1 ,2 0 1 4 ,3 : 2 and events coming from Northeast (Japan, China, Russia). The solid red line is the least square fit to our data (with the slope of the line), and the dashed line is a slope of 2.9 (thermal effect only). Standard deviation is 0.31 for the best-fit line and 0.36 for the 2.9 gradient.

DRAFT A p r i

  ] . M a n y seismological studies have been carried out in north-east Africa and Arabia to determine the depths of the Moho [e.g. Ahmed et al.,2 0 1 3 ;Mechie et al.,2 0 1 3 ]a n dt oi m a g eu p pe r mantle structure and understand regional geodynamics [e.g. Bastow et al.,2 0 0 5 ;Benoit et al.,2 0 0 3 ,2 0 0 6 ;Montagner et al.,2 0 0 7 ;Zhao,2 0 0 7 ;Sicilia et al.,2 0 0 8 ;Chang and Van der Lee,2 0 1 1 ;Koulakov , 2011] and how it is connected with global mantle flow [e.g. Montelli et al.,2 0 0 6 ;Boschi et al.,2 0 0 7 ,2 0 0 8 ;Forte et al.,2 0 1 0 ;Moucha and Forte, 2011]. However no previous study has the resolution in continental Yemen on the Gulf of Aden margin due to the lack of seismic stations. Surface wave studies [e.g. Debayle et al.,2 0 0 1 ;Sebai et al.,2 0 0 6 ;Montagner et al.,2 0 0 7 ;Sicilia et al.,2 0 0 8 ;Chang et al., 2011; Chang and Van der Lee,2 0 1 1 ]l a c ks u ffi c i e n tl a t e r a lr e s o l u t i o nt oi m a g et h ed e t a i l of lithospheric and uppermost mantle structures.

  triple junction, in the 'Horn of Africa'. In the Afar triple junction region, the presence of flood basalts [e.g. Baker et al.,1 9 9 6 ;Hofmann et al.,1 9 9 7 ;George et al.,1 9 9 8 ]a n da n abnormally low velocity upper mantle anomaly [e.g. Debayle et al.,2 0 0 1 ;Bastow et al., 2005] could be due to the presence of a mantle plume (Fig.1). Around 35 Ma ago, the current rift system began to open, under the influence of the Afar plume [e.g. Leroy et al., 2012]. The flood basalts of Ethiopia and Yemen are the signatures of voluminous eruptions produced during the Paleogene [Ebinger and Sleep,1 9 9 8 ]w i t hh i g h e s te r u p t i o nr a t e sa t 31 Ma [Baker et al.,1 9 9 6 ;Hofmann et al.,1 9 9 7 ;George et al.,1 9 9 8 ] . Ar e n e w e dp h a s e of volcanism took place 24 Ma ago that corresponds with the synchronous appearance of basaltic dikes and gabbros along the Red Sea, from Afar and Yemen to northern Egypt. From 25 to 16 Ma, a series of basaltic, trachytic and rhyolitic dikes were emplaced along the Red Sea margins [Zumbo et al.,1 9 9 5 ] ,a tt h es a m et i m ea se m p l a c e m e n to fl a r g e granitic batholiths (Jabal Hufash and Jabal Bura, see Fig.2). These granites are oriented north-south [George et al.,1 9 9 8 ]a n dl o c a t e da tt h eb o r d e ro ft h eG r e a tE s c a r p m e n t ,a 1000 km long sudden change in altitude, from 200 m on the west in the Tihama Plain, to more than 1000 m in the Yemeni highlands. In addition, syn-rift (30 to 16 Ma) basaltic flows dipping towards the sea (seaward dipping reflectors -SDRs) have been imaged under the Tihama Plain [Davison et al.,1994]. IntheRedSea,asecondphaseofopeningbegan 14 Ma ago synchronous with the formation of the Dead Sea transform fault to the north. .: TELESEISMIC TOMOGRAPHY IN SW ARABIA This phase is accompanied by increased extension rates and rift-flank uplift [Courtillot et al.,1 9 9 9 ] . From 12 Ma to present, magmatic provinces developed within a 2500 km radius of Afar and southwestern Yemen as far away as western Saudi Arabia, Jordan and northern Syria, [Zumbo et al.,1 9 9 5 ;Bertrand et al.,2 0 0 3 ;Coulié et al., 2003], see Fig.1. Recent magmatism has occured offshore in the ocean-continent transition of the Eastern Gulf of Aden margin [Lucazeau et al.,2 0 0 9 ;Autin et al.,2 0 1 0 ;Watremez et al.,2 0 1 1 ] ,o ff -a x i s Sheba ridge [d'Acremont et al.,2 0 1 0 ]a n db e l o wt h es o u t h e r nO m a nc o n t i n e n t a lm a r g i n

  2 0 1 0 ;[START_REF] Leroy | Contrasted styles of rifting in the eastern Gulf of Aden: A combined wide-angle, multichannel seismic, and heat flow survey[END_REF] but no syn-rift volcanism has been found east of the longitude 46 • E, see Fig.1 [Leroy et al., 2012]. The study region is near the edge of the African superplume, as shown by large scale global tomographic models [e.g. Debayle et al.,2 0 0 1 ;Sebai et al.,2 0 0 6 ;Boschi et al.,2 0 0 7 ,2 0 0 8 ]a n db y petrologically derived temperature estimates of the mantle [Rooney et al.collected from 23 temporary broadband stations deployed from March 2009 to March 2010 (YOCMAL project), and from one permanent station (DAMY, Geofon). The network extends from the Red Sea margin to the Gulf of Aden margin, passing

[

  Ahmed et al.,2 0 1 3 ] .Under stations located near the Red Sea (north-westernmost part of the network) and Gulf of Aden (southernmost part of the network) margins there are high velocity anomalies above 30 km due to the thinner crust (<30 km).[START_REF] Ahmed | Crustal structure of the rifted volcanic margins and uplifted plateau of Western Yemen from receiver function analysis[END_REF] estimated the thickness of the crust to be ∼22 km in the coastal areas, and less than 14 km for the Red Sea margin. These high velocity anomalies are located beneath SDRs, which is consistent with the emplacement of sub-aerial volcanic material during rifting[Tard et al.,1 9 9 1 ; Davison et al.,1 9 9 4 ;Bastow and Keir ,2 0 1 1 ;Leroy et al.,2 0 1 2 ] .

  several factors, such as temperature, chemical composition or anisotropy can affect the velocity of seismic waves [e.g. Karato,1 9 9 3 ;Sobolev et al.,1 9 9 6 ] .[START_REF] Karato | Importance of anelasticity in the interpretation of seismic tomography[END_REF] demon-

  the upper mantle [e.g. Bastow et al.,2 0 0 5 ;Basuyau et al.,2 0 1 0 ] . T h i sm e t h o dc o n s i d e r s relative delay times, so that problems associated with amplitude recovery (e.g. due to differing numbers of traveltime observations and regularization levels) and other artefacts associated with the inversion procedure (i.e. parameterization and ray-path accuracy) do not complicate the comparison of the data.

  [2004],Bastow et al. [2005] and Basuyau et al.[2010], we plot S versus P relative travel-time residuals, and find out that our data are consistent with a slope >2.9, thus implying the presence of partial melt along the rays coming from northeast. The presence of Holocene volcanoes (in Sana'a, Dhamar and Marib volcanic fields) on the surface helps support the idea that there is partial melt being created in the asthenosphere and intruding the lithosphere. Moreover, the isotopic signatures of the melts from the three volcanic fields suggest a strong asthenospheric component[Manetti et al.,1 9 9 1 ] ,w h i c hi sc o n s i s t e n tw i t ho u rr e s u l t s .We produce a relative velocity model for the propagation of P waves down to 300 km beneath western Yemen. The low velocity anomaly (∼ -4%) located below the Yemeni highlands between 0 and ∼35 km corresponds to a region of 30-35 km thick crust[START_REF] Ahmed | Crustal structure of the rifted volcanic margins and uplifted plateau of Western Yemen from receiver function analysis[END_REF]. At 30 km depth, we observe a dramatic transition from very low to high velocities (-4% to +4%) under the coast of the Red Sea and Gulf of Aden. We interpret this sudden short length scale (<40 km for the Red Sea margin, <100 km for the Gulf of Aden margin) variation as a lateral transition between a thick crust and a thinned, intruded and stretched crust. In addition, our new seismic images showing lower mantle velocities under the southwestern corner of Yemen is consistent with ongoing rifting above a thermal anomaly in the underlying mantle, due to the Afar plume. Under these conditions, melt generated by adiabatic decompression of the asthenosphere beneath thinned and stretched lithosphere migrates upwards to intrude or underplate continental crust and extrude as mostly basalt flows (Oligocene flood basalts) [White and McKenzie,1 9 8 9 ] . O u rp o s i t i v e Vp anomaly near the base of the crust under the Red Sea and Gulf of Aden margin are consistent with melt produced from an abnormally hot mantle which enriches the melt in MgO [White and McKenzie,1 9 8 9 ] ,c a p a b l eo fp r o d u c i n gi n t r u s i o n s / u n d e r p l a t eo fu pt o 7.2 km/s rather than 6.8 km/s from melting normal mantle. These high velocity anomalies are focused mainly into narrow zones of denser material, away from the most stretched areas. That is because the mantle temperature was highest at the start of continental break-up [White et al.,2 0 0 8 ] . S u c hl o w e rc r u s t a li n t r u s i o n s / u n d e r p l a t i n ga r eac o m m o n feature of volcanic margins such as the north Atlantic [Geoffroy,2005;Mjelde et al.,2005; the timing of the underplating is difficult. If emplaced before the eruption of SDRs and continental break-up as has been proposed for the north Atlantic (at least 10 to 15 Ma off Norway, Gernigon et al. [2006]) then it would have an influence on the structural development of the margin and partly consist of high pressure granulite/eclogite lower crustal rocks [Gernigon et al.,2 0 0 6 ] .

  )[Hammond et al.,2 0 1 3 ] . W es u r m i s et h a tas i m i l a rm e c h a n i s mi st h eb e s t explanation for our observations. This idea is supported by the large asthenospheric component in the magma from the three Yemeni volcanic fields inferred from trace element and isotope geochemistry[Manetti et al.,1 9 9 1] . O b s e v a t i o n so fr e c e n td i k ei n t r u s i o n sa t Harrat Lunayirr in Saudi Arabia [Pallister et al.,2 0 1 0 ]s h o wr i f t e dm a r g i nm a g m a t i s m may be important in accomodating extension after breakup along the full length of the Red Sea margin [Ebinger and Belachew ,2 0 1 0 ] . Our new relative P-wave velocity model beneath western Yemen suggests the young rift flank volcanoes on the margin of the Red Sea rift are caused by focused small-scale diapiric upwelling from a broad region of hot mantle beneath the area. Our work shows that relatively hot mantle, along with partial melting of the mantle, can persist beneath rifted margins after breakup has occurred. These findings have important implications for interpreting the thermal history and deformation of volcanic rifted margins after breakup is achieved since most breakup models assume rift margin volcanism ceases after seafloor spreading starts. Buoyant hot mantle may contribute to a dynamic topography of the Yemeni high plateau. Almost all the topography in this region could indeed have a dynamic origin, because the rift-flank uplift from flexure [Daradich et al.,2 0 0 3 ]i sn o ts u ffi c i e n tt op r o d u c eh i g h topographies over such a large area. Numerical modeling of the plume/lithosphere interaction predict an uplift of about 500 to 1500 m in less than 0.3 to 0.5 Ma after the plume initiation [d'Acremont

  p m DRAFT Oman), Geophysical Journal International, 180 (2), 501-519. Baker, J., L. Snee, and M. Menzies (1996), A brief Oligocene period of flood volcanism in Yemen: implications for the duration and rate of continental flood volcanism at the Afro-Arabian triple junction, Earth and Planetary Science Letters, 138 (1), 39-55. Bastow, I., and D. Keir (2011), The protracted development of the continent-ocean transition in Afar, Nature Geoscience, 4 (4), 248-250. Bastow, I., G. Stuart, J. Kendall, C. Ebinger, et al. (2005), Upper-mantle seismic structure in a region of incipient continental breakup: Northern Ethiopian Rift, Geophysical Journal International, 162 (2), 479-493. Bastow, I., D. Keir, and E. Daly (2011), The Ethiopia Afar Geoscientific Lithospheric Experiment (EAGLE): Probing the transition from continental rifting to incipient seafloor spreading, Volcanism and Evolution of the African Lithosphere, 478,5 1 -7 6 . Basuyau, C., C. Tiberi, S. Leroy, G. Stuart, A. Al-Lazki, K. Al-Toubi, and C. Ebinger (2010), Evidence of partial melting beneath a continental margin: case of Dhofar, in the Northeast Gulf of Aden (Sultanate of Oman), Geophysical Journal International, 180 (2), 520-534.

Figure 2 .

 2 Figure 2. (a) Topographic map of Southwestern Yemen. The Yemeni highlands, above 1000 m

Figure 3 .

 3 Figure 3. Azimuthal distribution of the events used on our study.

Figure 4 .

 4 Figure 4. Final P-wave velocity model obtained from inversion. Seismic stations are located

Figure 5 .

 5 Figure 5. Cross sections in the final P-wave velocity model obtained from inversion. (a) Cross

Figure 6 .Figure 7 .

 67 Figure 6. Crustal correction applied for the stations. These corrections were computed from

Figure 8 .

 8 Figure 8. Synthetic test for the propagation of crustal signal. (a) Synthetic output model, (b)

Figure 9 .

 9 Figure 9. Synthetic test for two low velocity anomalies under the Southwestern Yemen. (a)

Figure 10 .

 10 Figure 10. Plot of P-wave versus S-wave relative arrival-time residuals for all the stations

Figure 11 .

 11 Figure 11. Auxiliary material. Series of synthetic tests for constraining the low velocity

  Earth and Planetary Science Letters, 206 (3-4), 379-396. d'Acremont, E., S. Leroy, M. Maia, P. Gente, and J. Autin (2010), Volcanism, jump and propagation on the Sheba ridge, eastern Gulf of Aden: segmentation evolution and Figure 1. Geodynamic map of Arabia. Yellow triangles are for YOCMAL Network stations in southwestern Yemen. The magmatism older than 20 Ma is represented in purple, whereas the younger magmatism is in pink (modified from Davison et al. [1994]). J.: Jordan, UAE: United Arab Emirates, SS.: Shukra el Sheik Fracture Zone.

			X -21	X -29
	sphere composition, Chemical Geology, 198 (1), 47-61.
	Boschi, L., T. Becker, and B. Steinberger (2007), Mantle plumes: Dynamic models and
	seismic images, Geochemistry, Geophysics, Geosystems, 8 (10).
	Boschi, L., T. W. Becker, and B. Steinberger (2008), On the statistical significance of implications for oceanic accretion processes, Geophysical Journal International, 180 (2),
	correlations between synthetic mantle plumes and tomographic models, Physics of the 535-551.
	Earth and Planetary Interiors, 167 (3), 230-238. Daradich, A., J. Mitrovica, R. Pysklywec, S. Willett, and A. Forte (2003), Mantle flow,
	Calò, M., C. Dorbath, D. Luzio, S. Rotolo, and G. d'Anna (2008), WAM tomography dynamic topography, and rift-flank uplift of Arabia, Geology, 31 (10), 901-904.
	in the southern Tyrrhenian region. Petrological inferences and hypotheses on the fluid Davison, I., M. Al-Kadasi, S. Al-Khirbash, S., A. Al-Subbary, J. Baker, S. Blakey,
	circulation in the subducting Ionian slab and adjoining mantle domain, Bollettino di D. Bosence, C. Dart, R. Heaton, K. McClay, et al. (1994), Geological evolution of the
	Geofisica Teorica ed Applicata, 42 (2), 136-141. southeastern Red Sea Rift margin, Republic of Yemen, Geological Society of America
	Chang, S., and S. Van der Lee (2011), Mantle plumes and associated flow beneath Arabia Bulletin, 106 (11), 1474-1493.
	and East Africa, Earth and Planetary Science Letters, 302 (3-4), 448-454.
	Chang, S., M. Merino, S. Van der Lee, S. Stein, and C. Stein (2011), Mantle flow beneath
	Arabia offset from the opening Red Sea, Geophysical Research Letters, 38 (4), L04,301. Earth and Planetary Science Letters, 193 (3), 423-436.
	Coulié, E., X. Quidelleur, P. Gillot, V. Courtillot, J. Lefèvre, and S. Chiesa (2003), Com-Denèle, Y., R. Pik, and S. Leroy (2012), Thermal and tectonic histories of the Jabel-Bura
	parative K-Ar and Ar/Ar dating of Ethiopian and Yemenite Oligocene volcanism: im-per-alkaline granite on the south-eastern margin of the Red Sea (Yemen), Magmatic
	plications for timing and duration of the Ethiopian traps, Earth and Planetary Science Rifting and Active Volcanism Conference, Addis Abeba, Ethiopia.,p p .1 0 7 -1 0 8 .
	Letters, 206 (3), 477-492. Ebinger, C., and M. Belachew (2010), Geodynamics: Active passive margins, Nature
	Courtillot, V., C. Jaupart, I. Manighetti, P. Tapponnier, and J. Besse (1999), On causal Geoscience, 3,6 7 0 -6 7 1 .
	links between flood basalts and continental breakup, Earth and Planetary Science Let-Ebinger, C., and N. Sleep (1998), Cenozoic magmatism throughout east Africa resulting
	ters, 166 (3), 177-195. from impact of a single plume, Nature, 395 (6704), 788-791.
	d'Acremont, E., S. Leroy, and E. Burov (2003), Numerical modelling of a mantle plume: El-Anbaawy, M. I. H., M. A. H. Al-Aawah, K. A. Al-Thour, and M. Tucker (1992),
	Miocene evaporites of the Red Sea Rift, Yemen Republic: Sedimentology of the Salif
	DRAFT	6 p m A p r i l2 1 ,2 0 1 4 ,3 : 2 6 p m	DRAFT DRAFT

the plume head-lithosphere interaction in the formation of an oceanic large igneous

DRAFT A p r i

l2 1 ,2 0 1 4 ,3 : 2 6

p m DRAFT province, Debayle, E., J. Lévêque, and M. Cara (2001), Seismic evidence for a deeply rooted lowvelocity anomaly in the upper mantle beneath the northeastern Afro/Arabian continent,

Acknowledgments. This project was funded by the ANR YOCMAL (Agence Nationale de la Recherche), CNRS-INSU-PICS Yemen, GSMRB Yemen and is in the framework of the Actions Marges program. Seismic equipment from SEIS-UK is funded by NERC under loan number 873. We thank David Hawthorn, Alex Brisbourne and Victoria Lane for their efforts during the deployement and servicing of network, the French Embassy in Yemen (J. G. Sarkis, J. Dechezlepretre and C. Bousquet), local governers and the people of the Yemen governerates for their help during the field work.

X -23

halite: Sedimentary Geology, Sedimentary Geology, 81, 6 1 -7 1 . Evans, J., and U. Achauer (1993), Teleseismic velocity tomography using the ACH method: theory and application to continental-scale studies, in Seismic Tomography Theory and Practice, pp. 319-360, chapman and Hall, London.