François Lehmann

High-altitude free fall and parameter estimation for undergraduate numerical techniques laboratory

In this article, we present the development and treatment of an inverse problem applied to data from the Felix Baumgartner stratospheric jump. This jump is well documented with a lot of data. Therefore, it makes a particularly well-suited example for teaching in a numerical techniques laboratory. The major aim of the article is to give guidelines in order to construct a simple, but not simplistic, inverse problem with real data for junior undergraduate students. Students should master classical mechanics and have some skills in numerical modelling. We use the programming language Python and various libraries in order to build a model and solve the entire problem. This programming language is increasingly used and understood by students, which allows them to focus on the physical and numerical aspects of the involved problem. The fairly new strategy presented in this article is an attempt to estimate the angle of attack from acceleration measurements, and to give an uncertainty estimation of Baumgartner's free-fall speed.

Introduction

"Oh no, not the free-fall again!" Most first-ever physics courses start with classical mechanics, and especially with the study of the motion of solids. On 14 October 2012, the Red Bull Stratos team [1] shared with the world an extraordinary experiment, a leap from a capsule suspended 40km above Earth, over Roswell, New Mexico. During this jump, the skydiver Felix Baumgartner was the first human to break the sound barrier without any thrust force. He maintained for a total of 30 seconds a speed greater than the speed of sound in air.

Even many years later, the extraordinary experimental conditions of this free fall still make it a veritable 'playground' for any physicist who wants to teach classical mechanics and numerical modelling. Like any physics problem under extreme conditions, this experiment is more complex than it seems. It is a problem involving several scientific fields: mechanics, fluid dynamics, thermodynamics, numerical method and many others, such as pathophysiology or data acquisition. Anyone can understand it very quickly by reading the Red Bull report [START_REF] Bull | Full Scientific Data Report, Findings of the Red Bull Stratos Scientific Summit[END_REF]. This is probably the reason why it makes a particularly well-suited example for teaching. It allows students to make connections between different scientific fields, to break the disciplinary boundaries and to learn to work as a team on a small project. This example could be realized by undergraduate students after a set of laboratory classes, such as that proposed by Samsonau [START_REF] Samsonau | Computer simulations combined with experiments for a calculus-based physics laboratory course[END_REF]. To carry out this work in their third year, sophomore college students should have followed an introductory course in Python language, with an emphasis on applications in the physical sciences and engineering, on basic problem solving, programming techniques, and fundamental algorithms. Students who have experience in programming with Python of about 7 hours per week (1:30 lectures, 4h computer laboratory time and 1:30 personal work) for 6 to 7 weeks can easily tackle this problem.

Modelling such an experiment is difficult given the non-linearity of the processes in an environment where air pressure, air temperature and air density are interrelated, and change rapidly with altitude.

The mechanistic approach, which consists of describing phenomena in the context of the laws of conservation, is the most commonly used approach. It constitutes a powerful tool for understanding and seems best suited for predictive simulations. However, we should also not forget that the key to the modelling approach is to build a model that is simple, realistic and feasible. The most commonly used approach for modelling such a high-altitude free fall is to use Newton's second law to consider gravitational force, and an aerodynamic force that opposes the skydiver's motion through the air. This model has been used by many authors [4 -8] in order to estimate the position and velocity of the skydiver. An excellent introduction for students to this kind of problem can be found in the book written by Barger and Olsson [9]. This paper can be seen as a continuation of the earlier works of Colino and Barbero [START_REF] Colino | Quantitative model of record stratospheric freefall[END_REF] and Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF].

In order to build the model, we need to completely define the forces included in the equation of motion, consisting of all the parametric functions and the initial conditions. Most of the time, not all the parameters are well known. Sometimes, we don't know the exact value of a specific parameter. Some of these parameters cannot be directly measured, so in such conditions we can use parameter estimation techniques in order to estimate the parameter values and also, importantly, the parameter uncertainty and ultimately, calculating the associated uncertainty of the computed position and velocity of the skydiver, whilst taking into account the available measures. This paper follows a three-fold structure: first, building the direct problem in which the physical system of the free fall is modeled, with all parameters and with the initial conditions. Second, the inverse problem is merged with the direct problem and the data. Finally, the last section illustrates the global model application to determine an estimation of the parameters and their associated uncertainties.

II. The direct problem

In order to establish the equation of motion, Newton's second law is applied to a body of mass m subject to two forces, gravitational force 𝑃 ⃗ and drag force 𝐹 𝐷 ⃗⃗⃗⃗ :

𝑚𝑎 = 𝑃 ⃗ + 𝐹 𝐷 ⃗⃗⃗⃗ (1)
𝑎 is the body acceleration in m s -2 . If the free fall is strictly vertical and the positive z direction is chosen to be up, the equation of motion, as described in detail by Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF], is given by:

𝑚 𝑑𝑣 𝑧 𝑑𝑡 = -𝑚𝑔(𝑧) + 1 2 𝐶 𝐷 (𝑀𝑎)𝐴 ⊥ (𝛽)𝜌(𝑧, 𝑇)𝑣 𝑧 2 (2)
where 𝐶 𝐷 (𝑀𝑎) is the drag coefficient, which depends on the Mach number Ma, 𝐴 ⊥ is the projected area depending on the angle of attack 𝛽 (the angle between the z axis reversed and a longitudinal reference line on the body, head first), 𝜌(𝑧, 𝑇) the air density depending on altitude z and air temperature T, 𝑔(𝑧) the acceleration due to gravity depending on z. For this equation of motion, the state variables are the altitude z and the speed 𝑣 𝑧 . The aim of the direct problem is to compute the state variables, which are dependent variables, against time (the independent variable). This can be done by rewriting the second-order differential equation (2) as a system of two first-order equations:

{ 𝑑𝑧 𝑑𝑡 = 𝑣 𝑧 𝑑𝑣 𝑧 𝑑𝑡 = -𝑔(𝑧) + 1 2𝑚 𝐶 𝐷 (𝑀𝑎)𝐴 ⊥ (𝛽)𝜌(𝑧, 𝑇)𝑣 𝑧 2 (3)
To determine the evolution of the body, it is necessary to define its initial state, that is the initial position and the initial velocity. At this step, it is also necessary to define all the parametric functions which complete the model. The same relations as those used by Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF] are included within the model.

The parametric functions and the data

Before defining explicitly all the parametric functions, it is important to fix the geographical coordinate system in which the trajectory is studied. The initial GPS position (𝜑:latitude, 𝜆:longitude, z:altitude),

when Baumgartner left the capsule on 14 October 2012, at 18:06:32.0 UTC, was (𝜑 = 33.3408417°, 𝜆 = -103.7679067°, 𝑧 = 38969.4 m). To study the trajectory and the speed of Baumgartner during his fall, we choose a local tangent plane coordinates system with an origin given by (𝜑 = 33.3408417°, 𝜆 = -103.7679067°, 𝑧 = 0.0 m); by convention the east axis is labeled 𝑥 𝐸𝑎𝑠𝑡 , the north 𝑦 𝑁𝑜𝑟𝑡ℎ and the up 𝑧 𝑈𝑝 . In this local ENU (East, North, Up) coordinate system, the trajectory versus time is depicted in figure 1. This representation has not been used by other authors who studied the free fall, but it is very interesting, more intuitive and practical. In figure (1), we can clearly identify two stages of the leap's evolution. The first stage, from exit time t=0 to 75 s, is a quasi-perfect vertical one dimensional movement. The second stage from t=75 to t=260 s (main chute deployment, at 18:10:52.0 UTC) looks more like a classic free fall of a skydiver who does not try to increase his speed but instead tries to fall slowly and steadily in a nice regular "belly-to-earth" position.

In most papers [4,5 and 8], the acceleration due to gravity is taken as constant, but the change of g can be easily considered by using the world geodetic system ellipsoidal gravity formula (WGS 1984) with the first-order free-air correction factor:

𝑔(𝑧, 𝜑) = 9.7803 1+0.00193•𝑠𝑖𝑛 2 (𝜑) √1-0.00669•𝑠𝑖𝑛 2 (𝜑) -3.086 × 10 -6 • 𝑧 (4)
at the origin of local ENU system 𝜑 = 33.3408417° 𝑔(𝑧, 𝜑) = 9.7959 -3.086 × 10 -6 • 𝑧

For z varying from the initial position to the main chute deployment, 2566.8 ≤ 𝑧 ≤ 38969.4 m, the acceleration due to gravity is 9.6756 ≤ 𝑔 ≤ 9.7879 𝑚 𝑠 2 ⁄ , that is a relative error of 1% if g is taken as a constant. The total mass of Baumgartner taken for the simulations is 𝑚 = 121.2 𝑘𝑔 [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF].

For this experience, the Red Bull Stratos team did everything to ensure the safety of Baumgartner's jump [1; 10]. Various measures have been taken, analyzed and processed by the team to validate several world records, and now they are available to us. Baumgartner was equipped with several GPS apparatuses, and in particular a Garmin 18X-5 WAAS GPS with a sampling frequency of 5 position measurements per second. The data files made available also contain temperature, air pressure, speed of sound and triaxial acceleration data. Baumgartner was equipped with a triaxial accelerometer positioned at chest level in order to understand and analyze the environmental stressors experienced [START_REF] Garbino | Physiological monitoring and analysis of a manned stratospheric balloon test program[END_REF].

Atmospheric density and speed of sound (𝑀𝑎𝑐ℎ 1 ) profiles are estimated based on the pressure 𝑃(𝑧)

and temperature 𝑇(𝑧) measurements with:

𝜌(𝑧) = 𝑃(𝑧)
𝑅 𝑠 𝑇(𝑧) and 𝑀𝑎𝑐ℎ 1 (𝑧) = √1.40 • 𝑅 𝑠 𝑇(𝑧) (6 and 7)

with 𝑅 𝑠 = 287.0 𝐽 (𝐾 𝑘𝑔) ⁄ . The profiles of temperature, pressure, density and speed of sound are depicted in figure 2 (red curve) and are compared to two classical models, the US Standard atmosphere, 1976 [12] and NRLMSISE-00 Atmosphere Model [START_REF] Picone | NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues[END_REF][START_REF] Papitashvili | Web interface for NRLMSISE-00 Atmosphere Model, The Community Coordinated Modeling Center[END_REF], and to the sounding operated by the Red Bull Stratos team at the Santa Teresa observatory and Roswell airport. The most important data for studying the free fall is for 20 ≤ 𝑧 ≤ 40 km. In this part, we have good correlation between the data and the NRLMSISE-00 Atmosphere Model. In the later simulations we will use the data depicted on the red curve. To our knowledge, the data was obtained by the Red Bull Stratos team with a numerical weather forecast model.

For the drag coefficient and the projected area, we use the relations proposed by Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF]:

𝐶 𝐷 = 𝐴 for 𝑀𝑎 = 𝑣 𝑧 𝑀𝑎𝑐ℎ 1 ≤ 0.6 (8)
𝐶 𝐷 = 𝐴 + 𝐵(𝑀𝑎 -0.6) 2 for 0.6 ≤ 𝑀𝑎 ≤ 1.1

𝐶 𝐷 = 𝐴 + 0.25 • 𝐵 -𝐶(𝑀𝑎 -1.1) for 1.1 ≤ 𝑀𝑎 [START_REF] Blue | Overview of medical operations for a manned stratospheric balloon flight[END_REF] and

𝐴 ⊥ (𝛽) = 𝐴 𝑥 sin(𝛽) + 𝐴 𝑧 cos(𝛽) (11)
with 𝐴 𝑥 = 1.19 𝑚 2 the effective front area and 𝐴 𝑧 = 0.525 𝑚 2 the effective top area [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF], A, B and C three parameters that should be estimated in the inverse problem.

The last needed information, and probably the most difficult to obtain is the angle of attack 𝛽 during the leap. Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF], estimate an angle of attack 𝛼 = 𝜋 2 -𝛽 with the full video from Baumgartner's point of view during the leap. In this study, the angle of attack is estimated based on the triaxial acceleration measurements (fig. 3). Accelerometers are sensitive to both linear acceleration and the local gravitational field. In the absence of linear acceleration, the accelerometer output is a measurement of the rotated gravitational field vector and can be used to determine the accelerometer orientation angles. At the start of a fall, the vector sum of acceleration will tend toward 0 g (fig 3)this is called the phenomenon of weightlessness. Until t=25 s, the longitudinal and lateral accelerations do not change, after that time and until t=75 s, these accelerations fluctuate a lot. During this time, Baumgartner was in a dangerous situation of spinning, which has been studied in detail by Garbino et al. [START_REF] Garbino | Physiological monitoring and analysis of a manned stratospheric balloon test program[END_REF]. The vertical acceleration has a different signal: it changes gradually, probably due to the drag force's module and direction acting on the body (fig. 3). Even if we aren't in a perfectly static situation, and we also neglect the lateral acceleration 𝐴 𝑦 , and if we suppose likewise that Baumgartner didn't rotate around his body's z-axis (head first), the angle of attack is estimated with [15, eq. 55 page 19]:

𝛽 = 𝑐𝑜𝑠 -1 ( 𝐴 𝑧 √𝐴 𝑥 2 +𝐴 𝑧 2 ) (12)
with 𝐴 𝑧 and 𝐴 𝑥 representing the vertical and longitudinal acceleration respectively. In order to smooth the signal and remove potential linear acceleration, a rolling mean of 2.5 s is applied to the angle and used later in a linear interpolation lookup table for the simulations (fig. 4).

The numerical implementation

The programming language we use is Python, a language which gives readable code that closely resembles the equations. The equation (3) is a system of two first-order equations that can be solved by numerical integration. There are a lot of numerical tools called "ODE solver" (ordinary differential equation) to solve this kind of equation [START_REF] Linge | Programming for Computations -Python : A Gentle Introduction to Numerical Simulations with Python 3.6[END_REF]. Here, we have chosen to work with the SciPy ecosystem [START_REF] Virtanen | Fundamental Algorithms for Scientific Computing in Python[END_REF], and in particular with the module scipy.integrate.solve_ivp. This module allows us to choose between all the classical numerical integration methods, and trigger events to stop the integration time if needed. By default, we select the Runge-Kutta method. We can also set the maximum allowed time step size during integration, which allows us to properly control any variations in the parametric functions, like the angle of attack. At this stage, we have all the necessary information allowing us to program and to solve the direct problem, apart from the parameters A, B and C appearing in the drag force. A useful solution is to generate uniformly distributed random parameters between bounds and/or to take some values from literature and to build the code and do tests. This is an important step before tackling the inverse problem, since the direct code must be solved in as little CPU time as possible, and be reliable and robust regardless of the parameters used.

III. The inverse problem

The inverse problem covers a very wide range of formulations and fields of application. In this part, we will present an application of a classical non-linear inverse problem. Here, the aim is to find the model parameters we don't know that produce the data we have measured. From a mathematical point of view, we formulate this with the help of a classical cost function or objective function, i.e. the mean sum square error between model output and the measures:

𝑂(𝑝; 𝑦, 𝑦 𝑚𝑒𝑠 ) = 1 𝑛𝑚-𝑛𝑝 ∑ (𝑦 𝑖 (𝑝) -𝑦 𝑖 * ) 2 𝑛𝑛 𝑖=1 (13)
with 𝑝 = (𝐴, 𝐵, 𝐶) the vector of parameters, nm and np the number of measurements and parameters respectively, 𝑦 𝑖 (𝑝) the model output, e.g. the speed and 𝑦 𝑖 * the measurements. The aim is to find the parameters that minimize the objective function, which is done iteratively. Generally, the algorithm starts with a set of initial parameters and calculates the corrections to be made to the parameters, so that the objective function decreases until a criterion of the objective function is reached and/or until the parameter values don't change anymore. There are plenty of algorithms that can be used to minimize a function. The most widely used is probably the Levenberg-Marquardt algorithm. The choice of a good algorithm is problem dependent, on the number of parameters, the bounds on the parameters, the number of available measurements, the correlations between parameters, the structure of the model and so forth. We can distinguish two main classes of algorithms: global and local algorithms. According to the author's own experience, the most efficient strategy when solving a minimization problem is to be able to use several kind of algorithms and analyze the different results.

Here, we have chosen to use LMFIT: Non-Linear Least-Squares Minimization and Curve-Fitting for Python [START_REF] Newville | LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python (Version 1.0[END_REF]. This library allows us to use more than twenty methods of minimization, from the local Levenberg-Marquardt least squares method to the Markov Chain Monte Carlo (MCMC) algorithm, amongst others. A complete discussion of inverse methods is beyond the scope of this document; the aim here is simply to give students access to up-to-date methods and algorithms as "gray-box", and not completely "black-box".

The use of lmfit.minimize() module is straightforward, as long as the user starts with one of the numerous examples available on the GitHub repository ( https://lmfit.github.io/lmfitpy/examples/index.html) and adapts it to their problem. Principally, the user has to build a callable function that corresponds to the direct problem. This function should return the difference between model and measure, i.e. the residual vector:

𝜀 𝑖 = 𝑦 𝑖 (𝑝) -𝑦 𝑖 * (14)
In our example, the state variable of interest is the speed. 𝑦 𝑖 (𝑝) is the solution of the system of equation ( 3), i.e. 𝑣 𝑧 (𝑡 = 𝑡 𝑖 ), with 𝑡 𝑖 the discrete time observations for 𝑖 = 1 ⋯ 𝑛𝑚. In the local ENU coordinate system, we can compute the speed by using the evolution of the altitude z over time with a first order finite difference equation:

𝑦 𝑖 * ≡ 𝑣 𝑖 ≃ 𝑧 𝑖+1 (𝑡 𝑖+1 )-𝑧 𝑖 (𝑡 𝑖 ) 𝑡 𝑖+1 -𝑡 𝑖 (15)
The GPS used has a sampling frequency of 5 position measurements per second; for the whole duration of the leap (from t=0 to 260s), this gives us a total number of available measurements of 1300. In figure

(1), we have seen that the vertical free fall is well approximated during only the first part of the leap, that is for a period from t=0 to 75 s. For the parameter estimation, we only use the data obtained during this first part of the leap and with one measurement per second, that is a total number of measurements 𝑛𝑚 = 75, for a total number of parameters to be estimated 𝑛𝑝 = 3.

IV. Results and discussion

In order to start the inverse problem, it is often necessary to have some prior knowledge of the parameters we want to estimate. Sometimes, these parameters have no physical meaning, but even in this case we should give an initial set of parameters and some lower and upper bound values (Table 1). For this simulation, we use the trust region reflective method to minimize the objective function [START_REF] Newville | LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python (Version 1.0[END_REF]. At the end of the estimation, after 11 iterations, the final set of parameters and the 95% confidence interval on the parameters are: 𝐴 = 0.617 +/-0.035 (5.7%)

𝐵 = 0.840 +/-0.259 (30.8%)

𝐶 = 1.425 +/-0.344 (24.2%)

The results differ from those of Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF] for the parameters B and C. We find slightly larger values for almost identical confidence intervals. It should be noted that the parameters A and B have a correlation of -0.92 and for B and C of 0.89. Since we didn't use exactly the same angle of attack as Guerster and Walter [START_REF] Guerster | Aerodynamics of a highly irregular body at transonic speeds-Analysis of STRATOS flight data[END_REF], it is possible to find different parameters that give almost identical results for the speed. This is a key point for parameter estimation, as there is not always a unique solution to the inverse problem, but instead depends greatly on the model structure and on available data. The drag force in the model depends on the product of the drag coefficient and the projected area, so this creates a perfect correlation between them. Colino and Barbero's [START_REF] Colino | Quantitative model of record stratospheric freefall[END_REF] strategy was to use the product and take different values of this product depending on the stages of the leap. The results of the inverse problem for the altitude, speed and Mach number are depicted in figure [START_REF] Theilmann | Supersonic freefall-a modern adventure as a topic for the physics class[END_REF]. We can observe a fairly good match between the measurements and the modelling. The gray region (fig. 5) represents the 95% confidence interval of the estimated speed. In the second stage, from t=75 to t=260 s, we observe larger differences between the measurements and the model, probably due to a 3D effect on the trajectory. A more detailed graphic concerning the supersonic stage is given in figure [START_REF] Colino | Quantitative model of record stratospheric freefall[END_REF]. We found a total supersonic time of 30s, a maximum speed, reached at t=50 s, 𝑣 = 375.7 ± 10.8 𝑚 𝑠 ⁄ and a

Mach number 𝑀𝑎 = 1.25 ± 0.03. We can also see the impact of variations over time of the angle of attack on the speed, for t=40 s and t=74 s.

Now that the whole model has been built (direct and inverse), it is relatively simple to carry out additional tests if needed, for example to test other sets of initial parameters, or to assume that g is constant, or to propose a different model for the drag coefficient or for the air density equation, or to modify the time parametric function for the angle of attack.

If one wishes to have more reliable information on the parameters and their uncertainties, it is possible with LMFIT to use an MCMC algorithm to determine the posterior distributions for the parameters, and not only a local uncertainty estimation.

V. Conclusion

The good thing about building a full model, i.e. direct and inverse, is that it gives students the opportunity to learn multiple skills by discussing the problem on the conceptual level, together with mathematical formulation, finding numerical methods, appropriate algorithms and writing the computer code, and testing the programs and using them to perform numerical experiments, to analyze the problem, and to compare the results to real data. All of this is now possible, especially with Python if we agree to use some libraries as a "gray-box" and explain the principles of these algorithms to students. Additionally, representing the computational results graphically is always very important, and this is facilitated for students with Python [START_REF] Weber | The benefit of computational modelling in physics teaching: a historical overview[END_REF]. Finally, teaching the building of an inverse problem

  has a lot of similarities with creating an Investigative Science Learning Environment [20, figure 1.3 page 1-7 ].

Figure 1 :Figure 2 :

 12 Figure 1: Skydiver position evolution through the leap, expressed in a local tangent plane coordinates system, as origin (𝜑 = 33.3408417°, 𝜆 = -103.7679067°, 𝑧 = 0.0 m)

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Measured triaxial acceleration

Figure 6 :

 6 Figure 6: Measured and predicted altitude, speed and Mach number with the final set of parameters for the supersonic stage.

Acknowledgments

A special thanks to Jonathan Clark and Alex Garbino, Department of Medicine, Section of Emergency Medicine and Center for Space Medicine, Baylor College of Medicine, Houston -TX, for support and making the flight data available to us. A sincere thank you to our colleague, Rob Simmons, from Institut Universitaire de Technologie Louis Pasteur, for his diligent proofreading of the article.


Seismic evidence for partial melt below tectonic plates

Eric Debayle 1 , Thomas Bodin 1 , Stéphanie Durand 1 and Yanick Ricard 2 1 Univ Lyon, Univ Lyon 1, ENSL, CNRS, LGL-TPE, F-69622, Villeurbanne, France 2 Univ Lyon, ENSL, Univ Lyon 1, CNRS, LGL-TPE, F-69007, Lyon, France

The seismic low velocity zone (LVZ) of the upper mantle is generally associated with a low-viscosity asthenosphere that plays a key role for the dynamics of plate tectonics 1 . However, its origin remains enigmatic, some authors attributing the reduction in seismic velocity to a small amount of partial melt 2,3 , others invoking solid-state mechanisms near the solidus [4][5][6] , or the effect of volatile contents 6 . Observations of shear attenuation provide additional constraints to unravel the origin of the LVZ 7 . Here, we report the discovery of partial melt within the LVZ from the simultaneous interpretation of global 3D shear attenuation and velocity models. We observe that partial melting down to 150-200 km depth beneath mid-ocean ridges, major hotspots and back-arc regions feeds the asthenosphere. A small part of this melt (<0.30%) remains trapped within the oceanic LVZ. The amount of melt is related to plate velocities and increases significantly between 3 and 5 cm yr -1 , similar to previous observations of mantle crystal alignement underneath tectonic plates [START_REF] Debayle | Seismic observations of large-scale deformation at the bottom of fast-moving plates[END_REF] . Our observations suggest that by reducing viscosity [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] , melt facilitates plate motion and large-scale crystal alignment in the asthenosphere. Melt is absent under most of the continents.

Our finding results from the simultaneous analysis of two upper mantle tomographic models of shear wave velocity (Vs) and attenuation (parameterised with Qs the quality factor, a measure of energy dissipation). Until now, most global tomographic studies of the upper mantle and their thermochemical interpretation have focused on shear velocity 4,7 . Recent experiments on olivine suggest that wave speed and attenuation are insensitive to water, which would imply that elevated water contents are not responsible for the LVZ [START_REF] Cline | Redox-influenced seismic properties of uppermantle olivine[END_REF] . However, Vs is sensitive to temperature, composition and melt content, and deciphering the causes of its variations represents a strongly non-unique inverse problem, severely limiting our understanding of the Earth's interior. Shear attenuation has a different sensitivity to these quantities, and therefore provides complementary constraints on the origin of seismic heterogeneities [START_REF] Deschamps | Radial thermo-chemical structure beneath Western and Northern Pacific from seismic waveform inversion[END_REF] . Shear attenuation is negligibly dependent on major element chemistry [START_REF] Shito | Towards Mapping the Three-Dimensional Distribution of Water in the Upper Mantle From Velocity and Attenuation Tomography[END_REF] , and exponentially dependent on temperature [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] . The relation between attenuation and melt is debated, with some authors arguing for a weak dependence based on experiments, models and seismological observations [START_REF] Shito | Towards Mapping the Three-Dimensional Distribution of Water in the Upper Mantle From Velocity and Attenuation Tomography[END_REF] , while others suggest a larger effect 3 . Measuring shear attenuation is nevertheless difficult. For this reason, only a few global Qs models have been published in the last 20 years [START_REF] Romanowicz | 1.25 -Deep Earth Structure: Q of the Earth from Crust to Core[END_REF] , and the only recent joint interpretation of 3D Qs and Vs tomographic models at global scale is based on models built from different datasets and modeling approaches [START_REF] Dalton | Global seismological shear velocity and attenuation: A comparison with experimental observations[END_REF] .

The novelty of our study is to simultaneously interpret two recent global Vs and Qs models that are consistent as obtained from the same Rayleigh wave dataset, at the same resolution and using the same modelling approach. These Vs (DR2020s) and Qs (QsASR17) models are displayed in Figure 1. Details of our tomographic procedure can be found in Methods. Before interpreting these two models in the light of laboratory experiments, a few words are needed to emphasize in what temperature and pressure range our interpretation is pertinent. The attenuation models derived by mineral physicists 4,[START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] are valid for temperatures T larger than 900°C, which correspond to the base of the lithosphere and to the asthenosphere.

They consider thermally activated processes varying exponentially with 1/T that would imply quality factor in excess of 2,000 in the upper 100 km of the lithosphere and reaching several million near the surface (see Methods). However, finite Qs is observed in the crust and in the cold mantle lithosphere where attenuation is most likely related to non-thermal processes such as scattering or fluid-fracture interactions [START_REF] Romanowicz | 1.25 -Deep Earth Structure: Q of the Earth from Crust to Core[END_REF] . Furthermore, the attenuation observed in seismology accumulates along the seismic ray and the observation of Qs is only possible when the amplitude of a wave is measurably smaller than in a pure elastic model. Given the uncertainties on amplitude data, it appears impossible to resolve quality factors larger than ~2,000 with long period Rayleigh waves (at periods of 100s, assuming velocities of 4.5 km s -1 and ray lengths of 10,000 km, the amplitude reduction with respect to a pure elastic model, would be less than 3.5%). Our inversion leads indeed to a Qs model with strong lateral variations, by two orders of magnitudes, but with a maximum Qs of ~1,750. Our attenuation model is therefore mostly adapted to depths greater or equal to ~100 km, where Qs values between fifty and a few hundreds are expected [START_REF] Dalton | Global seismological shear velocity and attenuation: A comparison with experimental observations[END_REF] . Therefore, our interpretation applies to the oceanic asthenosphere and the mantle structure at depths greater or equal to ~100 km, where our tomographic models are acurate and where conditions similar to those used in laboratory experiments prevail. [START_REF] Xu | The effect of bulk composition and temperature on mantle seismic structure[END_REF] , given by two anelasticity models 4,[START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] based on laboratory experiments, appropriate for asthenospheric conditions (see Methods for details). These models only explain a limited part of the velocities and attenuations of our dataset. We consider that the first theoretical curve 4 is compatible with a given Vs-Qs observation if it falls within the typical uncertainties of the Vs and ln(Qs) observations, 1% and 10% respectively (ivory colour). We show in Methods and supplementary Figures S1-10 that reconciling the low values of Vs with Qs cannot be done by invoking radial anisotropy, elastically accomodated grain boundary sliding, or the effect of composition or water. However, this can be done by adding partial melt, thus reducing Vs, and shifting the theoretical curves to the left in Figure 2, since as discussed below, adding partial melt reduces Vs but has little effect on Qs. Warm colours in Figure 2 indicate the amount of melt (<0.7%) needed to reconcile observations with the first theoretical model 4 (dark blue curve), which requires the smallest amount of melt. On the right side of the ivory region, points in grey are those for which Vs is too high compared to the theoretical value. They are associated with the lithospheric depletion of cratonic roots. In this case, the grey intensity quantifies the departure in percent from the theoretical curve assuming a meltless pyrolitic mantle. Results using the second theoretical model [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] are shown in Figures S11-12. They lead essentially to the same conclusions, but require larger melt fractions (up to 1%), which are more difficult to reconcile with the very small melt fraction (~0.1%) suggested by geochemical studies 4 .

The effect of melt on Vs has been estimated to 7.9% reduction per percent of melt based on model calculations [START_REF] Hammond | Upper mantle seismic wave velocity: Effects of realistic partial melt geometries[END_REF] . Recent experiments 3 are in qualitative agreement but require a slightly larger Vs reduction (Figure S13). The effect of melt on attenuation is not well constrained and depends on the mechanisms of attenuation. We show in Methods that melt may have a large effect on Qs at short period (1 s), but not in the period range of surface waves (50-250s). We therefore neglect the effect of melt on Qs and we model its effect on Vs based on recent experiments 3 . Figure 2 shows that the slowest shear velocities require less than 0.7 % of melt.

Figure 3 displays global maps of melt content at different depths, using the same colour-coding as in Figure 2. The associated mantle temperatures, derived by our approach are on average slightly above the solidus in oceanic regions, between 100 and 200 km depth (Figure S14, panel g). The differences between the maps of temperature (Figure S14 panels a-f) and those of melt content are attributed to a variable amount of volatile [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] . A higher amount of volatiles above subduction zones may lower the solidus and favor melting, while a dryer mantle in other regions may impede melting. The heterogeneities in Figure 3 display a strong correlation with surface tectonics. Regions where Vs is too high (in grey) are located beneath continents down to 150 km depth. The discrepancy is likely due to our assumption of a homogeneous and pyrolitic mantle. This assumption is reasonable in the well-mixed convective mantle. However, beneath cratons, compositional heterogeneities and depletion of incompatible elements are contributing to the high seismic velocities. Our observations under cratons, at 150 km depth, are on average 2.4% faster than a pyrolitic mantle which is compatible with compositional effects [START_REF] Bruneton | Layered lithospheric mantle in the central Baltic Shield from surface waves and xenolith analysis[END_REF] .

At 100 km depth, melt is required below mid-ocean ridges, some hotspots near the Atlantic ridge and in the south Pacific Ocean, back-arc basins around the Pacific Ocean including the eastern margin of Asia and some other active tectonic regions (Afar, Tibet, West of North-America and Southwest Pacific, including the North Fidji basin and the Northfolk ridge south of the Vanuatu arc). In these regions, the amount of melt exceeds 0.3% and can reach 0.7%. Melt is not required beneath most of the remaining oceanic and continental Phanerozoic lithosphere where temperature variations alone explain our observations. The depth range 150-200 km corresponds to the oceanic LVZ, where a small amount of melt (<0.3%) is required over broad regions. The largest amount is beneath hotspots, ridges, and volatile-rich back-arc regions and can extend deeper than the LVZ, suggesting that these deep regions feed the asthenosphere with partial melt. For example, partial melt is observed down to 250 km beneath Hawaii and 300 km beneath the Afar and East-African rift, the hotspots located on the western part of North-America, the region of the Balleny Islands in Antarctica, the western Pacific and the Indian ocean near the Ninety-East ridge. A melting anomaly near 300 km on the eastern part of the Tibetan plateau is also observed.

The depth extent and amount of melt beneath oceanic regions are consistent with local studies. The NoMELT experiment [START_REF] Lin | High-resolution seismic constraints on flow dynamics in the oceanic asthenosphere[END_REF] was performed beneath a location relatively far away from Pacific hotspots (blue star in Figure 3) where we too confirm the absence of melt in what could be a volatile-poor mantle. We observe melting beneath the East Pacific Rise (green star in Figure 3) down to 250 km depth, in agreement with the MELT experiment [START_REF] Yang | Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone[END_REF] . The depth extent of melting beneath Hawaii down to 250 km is supported by measurements of water abundance in the deep region that feeds the plume [START_REF] Wallace | Water and partial melting in mantle plumes: Inferences from the dissolved H2O concentrations of Hawaiian basaltic magmas[END_REF] . Our results agree with observations of significant melting beneath the Philippine Sea plate and the Western Pacific 2 . Finally, melting down to at least 300 km beneath the Balleny Islands hotspot is also consistent with previous observations [START_REF] Sieminski | Seismic evidence for deep low-velocity anomalies in the transition zone beneath West Antarctica[END_REF] .

The origin of partial melt within most of the oceanic LVZ is uncertain. Melting at mid ocean-ridges can exceed 1% in the depth range 40-80 km [START_REF] Yang | Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone[END_REF][START_REF] Key | Electrical image of passive mantle upwelling beneath the northern East Pacific Rise[END_REF] , but smaller volume melting may occur down to depths of 150-250 km [START_REF] Key | Electrical image of passive mantle upwelling beneath the northern East Pacific Rise[END_REF] . After melting at midocean ridges, a small amount of melt may remain unextractable from the mantle peridotite [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF][START_REF] Faul | Melt retention and segregation beneath mid-ocean ridges[END_REF][START_REF] Selway | A small, unextractable melt fraction as the cause for the low velocity zone[END_REF] . Our observations suggest that decompression melting with at least 0.3 to 0.7% of melt occurs also beneath some hotspots and back arc basins down to about 200 km. Most of the melt produced in these regions is extracted and incorporated to the oceanic crust, but a small amount remains in the oceanic LVZ as it ages. The quantity of melt, if any, decreases close to continents below 0.1%, where a simple model without melt explains most of our observations. This may be due to the difficulty of melting the depleted continental lithosphere. However under tectonically active regions, beneath the Afar and East-African rifts, Tibet, Western North-America and Transantarctic Mountains, the asthenosphere contains a small amount of melt (Figure 3).

The quantity of melt needed to reconcile Vs and Qs under large swaths of the oceanic LVZ is larger than usual estimates of unextractable melt, which range from very small values 9 to a maximum of 0.1% [START_REF] Faul | Melt retention and segregation beneath mid-ocean ridges[END_REF] . Although melt is likely connected and able to percolate even at very small porosity [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] , surface tension resists the phase separation [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF][START_REF] Hier-Majumder | Role of grain boundaries in magma migration and storage[END_REF] . The ability of the melt to rise depends on various parameters, surface tension, buoyancy, permeability, melt and matrix viscosities, which are all known with large uncertainties [START_REF] Selway | A small, unextractable melt fraction as the cause for the low velocity zone[END_REF] . The melt fractions greater than 0.1% that we obtain are therefore plausible. This melt concentration is in overall agreement with the range of estimates derived from electromagnetic studies [START_REF] Ni | Electrical conductivity of hydrous basaltic melts: Implications for partial melting in the upper mantle[END_REF] of the LVZ, which often propose even larger values (see Methods). A structure made of magma-rich sills embedded in a meltless mantle might also be mapped by tomography as an average medium with moderate melt content. This partially molten layered model has been proposed for the northwest Pacific and Philippine plates 2 , but it could extend more generally to the entire oceanic LVZ. This layering would explain both the radial anisotropy observed within the oceanic LVZ [START_REF] Chang | Joint inversion for global isotropic and radially anisotropic mantle structure including crustal thickness perturbations[END_REF] and the sharp velocity and viscosity contrasts at the lithosphere-asthenosphere boundary 2 .

Finally, the amount of melt exhibits a very peculiar relation with the plate velocity [START_REF] Demets | Effect of recent revisions to the geomagnetic reversal time-scale on estimates of current plate motions[END_REF] expressed in a no-net-rotation reference frame (Figure 4). The melt fraction in the asthenosphere, abruptly increases by a factor close to 2 when the velocity is larger than 4 cm yr -1 . This variation is very similar to the observed variation of azimuthal anisotropy with present-day plate motion [START_REF] Debayle | Seismic observations of large-scale deformation at the bottom of fast-moving plates[END_REF] . Our results suggest that plate-scale crystal alignment beneath fast-moving plates is associated with a greater amount of melt. This requires either that melt facilitates deformation [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] or that deformation favours melt retention in the LVZ 2 , or both. In any case, the small amount of melt observed beneath large swaths of the oceanic LVZ is likely to significantly decrease viscosity (by one to two orders of magnitude [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] , see Methods) and to play a significant role in the decoupling of tectonic plates from the mantle. and light blue curves are the theoretical predictions assuming a pyrolitic composition [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] in the absence of melt using the anelasticity models of Takei 4 and Jackson et al. [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] , respectively. The upper colour scale indicates the amount of melt in percent required to explain our observations using Takei 4 's model. Adding melt enables to lower the predicted velocity without changing Qs. The lower grey scale indicates the misfit in percent between theory and observations, in regions where Vs is too high and cannot be reconciled with model predictions assuming a pyrolitic mantle. (The velocities [START_REF] Demets | Effect of recent revisions to the geomagnetic reversal time-scale on estimates of current plate motions[END_REF] are expressed in a no-net-rotation reference frame. The percentage of melt is averaged for all geographical points with similar plate velocities, using a sliding window of ±2 cm yr -1 . The amount of melt increases significantly in the asthenosphere (100-200 km) for plate velocities between 3 and 5 cm yr -1 . This result links with previous observations that only plates moving faster than 4 cm yr -1 can organize the flow at large scale in the underlying asthenosphere [START_REF] Debayle | Seismic observations of large-scale deformation at the bottom of fast-moving plates[END_REF] , suggesting that melt reduces viscosity [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] and facilitates large-scale crystal alignment.

Methods

Qs and Vs tomographic models. Figure 1 presents maps of DR2020s, a new global Vs model and QsADR17 a recent global Qs model 30 . Both models are built from the same massive Rayleigh wave measurements [START_REF] Debayle | A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements[END_REF] , and are obtained from a similar tomographic procedure. The first step is an automated waveform inversion approach which was applied to approximately 375,000 Rayleigh seismograms [START_REF] Debayle | A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements[END_REF] .

From a single surface wave seismogram, the waveform inversion extracts simultaneously a path-average depth-dependent shear velocity profile, Vs, and quality factor, Qs. By jointly interpreting the amplitude and phase of each waveform, we ensure that the interplay between Vs and Qs is accounted for, and that the shear quality factor and velocity profiles are constrained within the same period range. The waveform analysis is performed in the period range 50-250 s and accounts for the fundamental and up to the 5th higher mode of Rayleigh waves, thus ensuring a good depth resolution for Vs and Qs from 50 km depth down to the transition zone. It is a non-linear iterative process, which also produces frequency-dependent phase velocity and attenuation curves compatible with the recorded waveform. The effect of physical dispersion due to attenuation is accounted for in the modelling.

The second step is a regionalization of the 1D path-average models. DR2020s is obtained from the regionalization at each depth of the path-average Vs models using a continuous regionalization approach [START_REF] Debayle | A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements[END_REF] . This tomographic inversion yields 3D absolute velocities. QADR17, a global model of Rayleigh-wave attenuation was obtained after adapting the same regionalization approach to our dataset of Rayleigh wave attenuation curves, parameterized as ln(Q) [START_REF] Adenis | Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the Pacific Ocean[END_REF] Using other 1D or 3D models of radial anisotropy does not affect our conclusions, as discussed below. QsADR17 was obtained from the inversion at depth of QADR17 30 , using the same vertical smoothing as for DR2020s. DR2020s and QsADR17 are therefore consistent since they are derived from the simultaneous inversion of the same waveforms and inverted using the same regionalization approach with the same Gaussian filtering.

Prediction of Qs and Vs using a temperature-dependent model 4 . The dark blue curve in Figure 2 is obtained by predicting Vs and Qs at each depth for a range of temperatures, using a temperature-dependent anelasticity model 4 and assuming a pyrolitic composition [START_REF] Xu | The effect of bulk composition and temperature on mantle seismic structure[END_REF] . This model parameterizes the relaxation spectrum with a small number of variables determined from experimental data. It accounts for a monotonic background spectrum plus a broad temperature-dependent absorption peak in the seismic frequency band, whose amplitude and width increase below the solidus. This produces a significant enhancement of polycrystal anelasticity before melting. This pre-melting effect induces large Vs reductions under the solidus. We choose the solidus according to Hirschman [START_REF] Hirschmann | Mantle solidus: Experimental constraints and the effects of peridotite composition[END_REF] . We use the model parameters given in Tables 1 and2 of Takei 4 , except for the unrelaxed shear modulus 𝜇 ! . We estimate 𝜇 ! for a pyrolitic mantle using the mineralogic phase diagram computed by the Perple X software [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] . In Figures S15-S16 we show the results obtained with the shear modulus of Takei 4 , which reduces the amount of melt in oceanic regions but fails to explain our ln(Qs)-Vs observations beneath old oceanic basins at 100 km depth. It also significantly increases the misfit beneath continents.

Prediction of Qs and Vs from the experimental results of Jackson et al. [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] . The light blue curve in Figure 2 is obtained by a two-step process. First, we use experimental results on melt-free polycrystalline olivine to predict a theoretical quality factor Qt. Computations are performed in the temperature range 800-1800 K and the pressure range 1.49-12.99 GPa, corresponding to the depth range 50-390 km. We use the following relation [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] :

𝑄 ! !! 𝑇 ! , 𝑇, 𝑃, 𝑑 = 𝐴 𝑇 ! 𝑑 !! exp -(𝐸 + 𝑃𝑉) 𝑅𝑇 ! Eq. A.1
In this formulation, 𝑄 ! !! at period T0 depends on the temperature T, the pressure P and to a lesser extent the grain size d. R is the gas constant, E is the activation energy and V is the activation volume. Table S1 summarizes the values of the different parameters of Eq. A.1. In the range of values compatible with experiments [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] , we choose T0=100 s, the average period of our long period Rayleigh waves and d=0.01m. The theoretical curves corresponding to Eq. A1 are shown in Figure S9, with the effect of changing grain size.

Second, we use Perple X [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] to estimate an isotropic Vs for a pyrolitic model [START_REF] Xu | The effect of bulk composition and temperature on mantle seismic structure[END_REF] . We compute Vs for the same temperatures and pressures as Qt. Perple_X produces unrelaxed S-wave velocities at infinite frequency, above the absorption band of seismic attenuation, while our long period (>50 s) seismic surface waves are likely to see relaxed velocities. We therefore correct for the effect of anelasticity on Vs [START_REF] Karato | Importance of anelasticity in the interpretation of seismic tomography[END_REF] :

𝑉 ! ! = 𝑉 !! 1 - ! ! !! (! ! ,!,!,!) !!"# !" ! Eq. A.2
where 𝑉 !! is the unrelaxed velocity and 𝑉 ! ! is the velocity corrected from the effect of attenuation. According to Eq. A2, the stronger the attenuation, the stronger the velocity reduction is. The light blue curve in Figure 2 is finally obtained by plotting at a given depth the obtained ln(𝑄 ! ) as a function of 𝑉 ! ! for each temperature.

Comparison between our observations and the theoretical relations. At each geographical location and depth, we extract the observed value of Qs in QsADR17

and Vs in DR2020s and we compare our observations with the chosen theoretical relation. We consider that a theoretical curve is compatible with a given Vs-Qs observation if it falls within ±1% and ±10% of the Vs and ln(Qs) observations, respectively (Figure 2). These errors are on the conservative side of recent estimates [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF][START_REF] Resovsky | Error bars for the global seismic Q profile[END_REF] . If the theoretical relation cannot explain our observations, two situations can arise: 1) the observed Vs is lower than the theoretical V ! ! . In this case melt can be added to reduce V ! ! . This is done assuming melt has no effect on Qs; 2) the observed Vs is higher than V ! ! . This is mostly the case under continents and can be explained by the depletion of the pyrolitic mantle increasing the theoretical

V ! ! 18 .
Effect of grain size. Grain sizes in the shallow upper mantle are likely 1-20 mm as observed in upper mantle xenoliths, harzburgite and dunite bodies [START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] . The effect of grain size on Eq. A1 is shown in with a smaller amount of melt. However, in both cases, amounts of melt are larger than for our preferred model using the equations of Takei 4 (Figure 3).

Amounts of melt.

Although the amount of melt present in the LVZ depends on the choice of the anelasticity model and on the effect of melt on Vs, we show in Figures S12, S16, S17, S18 that different choices of parameters yields similar maps and that melt fractions away from source regions always stay below 1%, in agreement with the partially molten asthenosphere model 2 . We note however that our preferred anelasticity model 4 generally explains our Vs and ln(Qs) observations without the need of partial melt above and below the asthenosphere (Figure 3). This is therefore a conservative choice, able to reconcile our Vs and ln(Qs) observations in most of the upper mantle, while minimizing the amount of melt required in the LVZ.

Effect of radial anisotropy. We used the anisotropic parameter ξ=(Vsh/Vsv) 2 of PREM to convert our Sv observations into an isotropic model, DR2020s. We checked that 3D models of ξ obtained from anisotropic models like SEMUM2 38 or S362ANI [START_REF] Kustowski | Anisotropic shear-wave velocity structure of the Earth's mantle: A global model[END_REF] would not affect our main conclusions. As an example, we computed the anisotropy ξ that would be needed to account for the discrepancy between observed Qs and Vs, without invoking mantle depletion or partial melt (Figure S3). Comparison of this figure with Figure S2 shows that neither the patterns nor the amplitude would fit (a 3 to 4 times larger anisotropy than in SEMUM2 would be needed to explain our observations). The anisotropy therefore does not affect our conclusions.

Elastically accommodated grain boundary sliding (EAGBS) hypothesis.

Anelastic relaxation caused by EAGBS can produce a sharp velocity reduction and may explain the LAB and the mid-lithosphere discontinuity observed beneath continents [START_REF] Karato | Mechanisms and geologic significance of the midlithosphere discontinuity in the continents[END_REF] . We tested the EAGBS hypothesis as an alternative explanation to partial melting for reconciling Vs and Qs in the LVZ. EAGBS is characterized by a sharp attenuation peak at a characteristic frequency 𝜔 !"#$% , followed by a diffused absorption band at lower frequencies. We first compute the characteristic frequency for the transition between unrelaxed and relaxed shear moduli. It is given by Karato et al., [START_REF] Karato | Mechanisms and geologic significance of the midlithosphere discontinuity in the continents[END_REF] :

𝜔 !"#$% = 𝐴. 𝑑 !! . 𝑒𝑥𝑝 - ! * !!! * !.! . 1 + ! ! ! !! ! Eq. A3
where Cw is the water content in wt%, A=2.3x10 -14 m s -1 , Cw0=10 -4 wt%, r=1, E * =350 kJ mol -1 , V * =10 -6 m 3 , d=5x10 -3 m. Following Cobden et al. 7 we assume a dry mantle with Cw=10 -6 wt %.

Surface waves analysed at frequencies lower than 𝜔 !"#$% should see relaxed The Ghahremani [START_REF] Ghahremani | Effect of grain boundary sliding on anelasticity of polycrystals[END_REF] equation, which produces the smallest velocity reduction is :

𝑉 ! !"#$% = 𝑉 ! ! . !.!"!!.!"! !!! • !!! !.!"!!.!"! Eq. A.4
where 𝜈 is the Poisson's ratio (~0.3 for upper mantle minerals). Regions where this can reconcile our Vs and Qs observations are displayed in red in Figure S4.

The two other formulations give too large velocity reduction. Figure S4 shows that EAGBS fails at reconciling our long period observations, excepted in a few localized regions (in red). This does not preclude EAGBS to play a significant role in explaining the sharp velocity reduction associated with the LAB, which cannot be resolved by long period surface waves. However, the strong velocity reduction observed within the LVZ requires another mechanism.

Effect of composition.

Compositional heterogeneities exist in the upper mantle, likely in the domain limited by pyrolite, harzburgite and pyroxenite [START_REF] Hirschmann | A possible role for garnet pyroxenite in the origin of the '"garnet signature[END_REF] . Here, we check the effect of changing the composition from pyrolite to these end-members using Perple X [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] to compute the velocities. Pyroxenites have variable bulk compositions [START_REF] Lambart | Markers of the pyroxenite contribution in the majorelement compositions of oceanic basalts: Review of the experimental constraints[END_REF] and we consider both a silica-deficient composition (see, Stixrude and Lithgow-Bertelloni 44 , Table 1) and a silica-excess composition similar to that of a subducted basalt (see, Xu et al., [START_REF] Xu | The effect of bulk composition and temperature on mantle seismic structure[END_REF] , Table 1). A silica-excess pyroxenite would only be faster that pyrolite after eclogitization, below 150 km depth (Figure S5, orange curve). At 100 km depth, a mantle made entirely of silica-excess pyroxenite would be too slow to explain the observed velocities (Figure S6). At larger depth, after eclogitization, it would lead to large melt fractions (Figure S6). The silicadeficient pyroxenite and the harzburgite are characterized at asthenospheric conditions by much faster or comparable velocities than pyrolite (Figure S5) and hence cannot explain the low velocities that we attribute to the presence of melt.

Harzburgite would give melt fractions and distribution similar to pyrolite (Figure S7), while the silicate-deficient pyroxenite would lead to higher melt fractions (up to 1.7%) and a ubiquitous melt presence below 150 km depth (Figure S8).

Pyroxenite, whether enriched or depleted in silica, should only be present in a small proportion (~5%) [START_REF] Hirschmann | A possible role for garnet pyroxenite in the origin of the '"garnet signature[END_REF] that hardly affects the seismic velocities of the peridotitic mantle. At any rate, large-scale chemical heterogeneities (> 1000 km)

are not likely to exist in the LVZ, where viscosity is small (below 10 20 Pa s) and the convective mantle should be well stirred and mixed. Therefore, our observations cannot be explained by reasonable variations of compositions in the LVZ.

Possible effect of water. Although, it has been recently shown that in olivine, wave-speed and attenuation are insensitive to water [START_REF] Cline | Redox-influenced seismic properties of uppermantle olivine[END_REF] , some previous studies had suggested an effect of water quantified following the equation [START_REF] Behn | Implications of grain size evolution on the seismic structure of the oceanic upper mantle[END_REF] :

𝑄 ! !! 𝑇 ! , 𝑇, 𝑃, 𝐶 !" , 𝑑 = 𝐴 𝑇 ! 𝑑 !! ! !" ! !"(!"#) ! exp -(𝐸 + 𝑃𝑉) 𝑅𝑇 ! Eq. A.5
where COH is the water concentration, r is a dimensionless constant and all other parameters are as in Eq. A1. The average H 2 O content of Earth's upper mantle is estimated in the range 50-200 ppm [START_REF] Hirschmann | Water, melting, and the deep Earth H2O cycle[END_REF] . We test a value close to the upper bound increase of water content must be balanced by a decrease of temperature and therefore by a higher expected velocity 𝑉 ! ! . Reconciling this higher 𝑉 ! ! with our observations requires larger amounts of melt (Figure S10). Therefore, if there were an effect of water on Vs and Qs, it would increase the required amount of melt.

Effect of melt on attenuation.

A small amount of melt might have a large effect on attenuation [START_REF] Faul | Shear wave attenuation and dispersion in meltbearing olivine polycrystals: 2. Microstructural interpretation and seismological implications[END_REF] in the case of grain boundary sliding, where a broad dissipation peak is observed in the seismic frequency band. Body wave studies have reported low Qs values at short period (1s) with Qs≤25 beneath the Juan de Fuca and Gorda ridges [START_REF] Eilon | High seismic attenuation at a mid-ocean ridge reveals the distribution of deep melt[END_REF] and 25≤Qs≤80 in the back-arc mantle of Central America, the Marianas and the Lau Basin [START_REF] Abers | Reconciling mantle attenuation-temperature relationships from seismology, petrology, and laboratory measurements[END_REF] . A recent experimental study suggests that 0.2% of melt produces Qs=48 3 . However such a large effect of melt on Qs is difficult to reconcile with long period (>30 s) surface waves observations of Qs=80 beneath the fast spreading southern East pacific Rise [START_REF] Yang | Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone[END_REF] where up to 1-2% of melt is expected [START_REF] Dunn | Imaging the transition between the region of mantle melt generation and the crustal magma chamber beneath the southern East Pacific Rise with short-period Love waves[END_REF] . If the attenuation mechanism is melt squirt, then the dissipation peak may lie outside the seismic frequency band and cause little attenuation. Model calculations [START_REF] Hammond | Upper mantle seismic wave attenuation: Effects of realistic partial melt distribution[END_REF] , longperiod surface waves [START_REF] Yang | Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone[END_REF] and experiments 12 favour such interpretation. This is also supported by our long period (100 s) global surface wave tomographic models, which suggest that strong Vs and Qs reductions are not necessarily correlated 30,[START_REF] Adenis | Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the Pacific Ocean[END_REF] .

For these reasons, we assume that in the period range of surface waves (50-250s), melt does not significantly affects attenuation and we neglect its effect on Qs.

Temperature in the upper mantle. Our Qs and Vs seismic models are interpreted using mineral physics results 4,[START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF] . These experimental results are valid at temperatures larger than 900°C or 1000°C, which are reached at the base of the lithosphere and in the asthenosphere. The predictions of our model can be checked with respect to the wellknown thermal behaviour of the cooling oceanic lithosphere [START_REF] Turcotte | Geodynamics: Applications of Continuum Physics to Geological Problems[END_REF][START_REF] Stein | A model for the global variation in oceanic depth and heat flow with lithospheric age[END_REF] although this is certainly pushing the model outside its applicability range. with geodynamic and petrologic expectations. In agreement with our findings, under oceans, the average 1D temperature profile is above the solidus 4 and the adiabat [START_REF] Katsura | Adiabatic temperature profile in the mantle[END_REF] in the depth range 100-200 km, where we predict partial melting. This overshoot of the temperature, above the adiabat, is indeed found in all numerical simulations of mantle convection [START_REF] Curbelo | Numerical solutions of compressible convection with an infinite Prandtl number: comparison of the anelastic and anelastic liquid models with the exact equations[END_REF] . Under continents, the temperature appears to be below the solidus 4 and the adiabat [START_REF] Katsura | Adiabatic temperature profile in the mantle[END_REF] . It reaches the adiabat around 250 km depth (Figure

S14, panel g).

Compatibility of our results with electrical conductivity. Our results are consistent with the interpretation of a number of recent studies of electrical conductivity. In the depth range 100-150 km, the joint interpretation of electrical and seismic data 56 requires 0.3-2.5% of melt beneath the mid-Atlantic ridge, 1% or less melt beneath Hawaï and less than 5% of melt beneath the East-pacific rise, in good agreement with Figure 3. We also confirm the absence of melt in the region of the NoMELT experiment [START_REF] Lin | High-resolution seismic constraints on flow dynamics in the oceanic asthenosphere[END_REF] , in agreement with conductivity data in this region [START_REF] Selway | A small, unextractable melt fraction as the cause for the low velocity zone[END_REF] .

Previous experimental results suggest that between 0.3 and 2% of hydrous basaltic melt can account for the observed electrical conductivity in the LVZ [START_REF] Ni | Electrical conductivity of hydrous basaltic melts: Implications for partial melting in the upper mantle[END_REF] . A more recent study [START_REF] Freitas | Simultaneous measurements of electrical conductivity and seismicwave velocity of partially molten geological materials: effect of evolving melt texture[END_REF] has refined these results by simultaneously measuring wave velocity and electrical conductivity on a simplified partial melt analogue. They conclude that the low velocity zone away from spreading ridges can be explained by 0.3-0.8% volatile-bearing melt, the upper bound of our observations.

Implications for viscosity. The variation of viscosity as a function of melt content occurs in two steps [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] . The onset of melting brings already a significant effect on viscosity when a connected network of melt tubules is formed. The viscosity η is expected to decrease by one or two orders of magnitude before the melt fraction φ reaches 0.1%. For larger melt content, η decreases further [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF] with dln(η)/dφ = -26 but this effect is minor for the low melt content that we observe. Under oceanic plates where φ is around 0.3%, η should be one or two orders of magnitude lower than under continents. Fig. S3: Dependence of shear velocity on melt fraction 𝜑. The linear Vs reduction of 7.9% per percent of melt 2 is shown in black from a reference velocity Vref. The polynomial expression ( 𝑉 ! = 0.065𝜑 ! -0.5565𝜑 + 𝑉 !"# ) derived from experimental results 3 is shown in red. The anelastic effect expected for seismic waves at high temperature 4 is shown assuming Qs=80, the value of PREM 5 in the asthenosphere, for two values of α (blue curve for α=0.38 and green curve for α=0.26). For small melt fractions (<1%), a stronger Vs reduction is obtained using the polynomial expression and the choice of α has a small effect. Vs (km/s) Fig. S4: same as Fig. 3 but using a linear 7.9% Vs reduction per percent of melt 2 instead of the polynomial expression derived from experimental results 3 . The upper color scale is slightly modified to allow melt content up to 1% (the maximum value at 100 km depth). Fig. S5: same as Figure 2 but the unrelaxed shear modulus 𝜇 ! needed to compute the temperature-dependent model (dark blue curve) is calculated using parameters proposed for of the temperature model of the Pacific 6 , instead of those deduced from Perple X assuming a pyrolitic composition.

Figure 2 displays

 2 Figure 2 displays Qs as a function of Vs for each pixel of the maps, at different depths in the upper mantle. The curves in dark and light blue represent the theoretical relations due to temperature variations for a meltless pyrolitic

  mantle

Fig. 1 :

 1 Fig. 1: Shear velocity and attenuation in the upper mantle. Panels a-c-e-g: perturbations in shear wave velocities from DR2020s shown in percent with respect to a mean value Vref given in km s -1 above the colour scale. Panels b-d-f-h: maps of our Qs tomographic model QsADR17 at different depths in the upper mantle. Qs is plotted with a logarithmic scale. Its geometric average is given above the colour scale. Hotspot locations are shown with black triangles.

Fig. 2 :

 2 Fig. 2: Scatter plot depicting the observed shear attenuation as a function of shear velocity compared with theoretical predictions. The velocities (from DR2020s) and attenuation values (from QsADR17) are plotted at 100 km (panel a), 150 km (panel b) and 200 km depth (panel c). Each dot corresponds to a geographical location. The dark blue

Fig. 3 :

 3 Fig. 3: Melt content at different depths in the upper mantle. This melt content is derived from the joint interpretation of QsADR17 and DR2020s. The colour coding is identical to Figure 2. Melt content in percent is indicated with warm colour from ivory (0% melt) to brown (0.4 to 0.7% melt). The grey scale indicates the misfit in percent between theory and observations, in regions where Vs is too high compared with predictions. Hotspots locations are indicated with black triangles. The blue and green stars indicate the location of the NoMELT and MELT experiments, respectively.

Fig. 4 :

 4 Fig. 4: Percentage of melt at different depths as a function of absolute plate velocity.

  . The regionalisation of the path-average attenuation curves accounts for frequency-dependent effects like focussing-defocussing, which can have important effects on the amplitude of Rayleigh waves. The logarithmic parameterization brings the distribution of the quality factor dataset close to a Gaussian, allows the large variations of Q documented by local seismic studies and guarantees to avoid negative attenuation values in the inverted model. The horizontal smoothing in DR2020s and QADR17 is determined by a Gaussian a priori covariance function controlled by an angular correlation length. Adenis et al. 32 chose a conservative value of 10° (meaning that the Q model is resolved accurately up to spherical harmonic 12). Debayle and Ricard 31 used a shorter correlation length of 3.6° for their Vs model DR2012. We re-inverted their dataset using a correlation length of 10° in order to obtain a Svwave tomographic model at the same horizontal resolution and vertical smoothing than QsADR17. To minimize biases due to un-modelled radial anisotropy, we computed the isotropic Voigt average of our Sv model and the Sh model obtained by adding the radial anisotropy of PREM to our Sv model. The resulting isotropic model DR2020s is plotted in Figure 1.

Figure S9 .

 S9 FigureS12with a smaller amount of melt. However, in both cases, amounts of

Figure S1 is similar to Figure 3

 3 but based on our Sv observations corrected with the anisotropic ξ parameters of SEMUM2 (FigureS2). Only minor changes are observed (for example in the Pacific or Indian oceans at 150 and 200 km), but none that would change the interpretations of this paper. The lateral variations of ξ in SEMUM2 (FigureS2) are too small to affect our results. To confirm this point,

  velocities within the absorption band (see Fig. SI-2 of Karato et al.[START_REF] Karato | Mechanisms and geologic significance of the midlithosphere discontinuity in the continents[END_REF] ), providing depth and temperature are within the activation ranges of EAGBS. EAGBS is activated between 60 and 160 km depth in the temperature range 900-1,350 K. We consider that EAGBS may affect the interpretations of our long period (>40 s) Qs and Vs maps when 𝜔 !"#$% > 0.025 Hz and, from figure4of Karato et al.[START_REF] Karato | Mechanisms and geologic significance of the midlithosphere discontinuity in the continents[END_REF] , in the following temperature ranges for the depths of our tomographic inversion: • 70 km: 930-1,140 K • 90 km: 980-1,180 K • 100 km: 1,010-1,220 K • 125 km: 1,060-1,270 K • 150 km: 1,130-1,330 K Regions where EAGBS can apply are displayed in green in Figure S4. These regions correspond mostly to continental areas. Assuming that EAGBS correction applies, Karato et al. 40 list three alternative formulations for the velocity reduction.

(COH=1000H/10 6

 6 Si, corresponding to 125 ±75 ppm by weight of water) with r=1, a standard value for moderate water content and COH(ref)=50H/10 6 SI 7 . The corresponding theoretical curves are shown in blue on Figure S9. At fixed temperature, water increases attenuation. However to keep Q constant, an

  Figure S19 (panel a),displays the temperature variations beneath the Pacific as a function of sea-floor age predicted from our seismic models using experimental results4 . It demonstrates that even at depths shallower than 100 km, we retrieve the well-known age-dependence of temperature in oceanic regions (we also pick up the cold signal of the West Pacific subductions at old ages). However, our temperature variations between the ridge and an old lithosphere (~100°) are lower than the predictions of the plate-cooling model (250° at 75 km depth, panel b). Panel c displays the quality factors that would be deduced from the plate-cooling model using the experimental results 4 used in this paper. In the oceanic lithosphere, Qs values much greater than 2000 (ln(Qs) > 7.6) would be predicted, that cannot be retrieved by long period seismology. The limitations of experimental data at low temperatures together with the inability of surface wave seismology to quantify precisely negligible attenuations, explain why we do not interpret results at depths shallower than 100 km. Below the oceanic plates, FigureS14, panel g shows that our predicted temperatures beneath oceans are perfectly compatible

FigFig. S7 :

 S7 Fig.S6: same as Figure3but instead to estimate the unrelaxed shear modulus 𝜇 ! using Perple X and a pyrolitic model, we use fitting parameters of the temperature model for the Pacific 6 .

Fig

  Fig. S8: same as Figure S2 but for a grain size of 100 mm.

Table S1 :

 S1 Reference parameters for Eq.A.1 after Jackson et al.1 Fig.S1: same as Figure2but the color scales indicate the departure from the model of Jackson et al.1 (light blue curve), for a grain size of 10 mm. The upper color scale indicates the amount of melt in percent required to explain our observations (ivory color from 0 to 0.01% of melt underlines data for which the model can reconcile our Qs and Vs observations). The lower grey scale indicates the misfit in percent between theory and observations, in regions where Vs is too high and cannot be reconciled with model predictions. Fig. S2: same as Figure 3 but color scales indicate the departure from the model of Jackson et al.1 , for a grain size of 10 mm. Melt content in percent is indicated with warm colors (upper color scale) from ivory (0% melt) to brown (0.4 to 0.7% melt). The lower grey scale indicates the misfit in percent between the theory and observations, in regions where Vs is too high compared with predictions.

	0100 km					0150 km
	0200 km					0250 km
	0300 km	ln(Qs)				0350 km
							0200 km
		4.0	4.2		4.4	4.6		4.8
					Vs (km/s)
		0.00	0.01	0.10	0.15	0.20	0.30	0.40	1.00
	0.00	0.01	0.10		0.15	0.20		0.30	0.40	1.00
		0	1		2	3	4		5
	0	1			2		3		4	5

Acknowledgements. We thank the Iris and Geoscope data centers for providing seismological data. We also thank two anonymous reviewers for their comments, J.P. Perrillat and M. Behn for discussions on mineralogy and attenuation models, and F. Dubuffet for preparing data for sharing as IRIS Data Products. The European Union Horizon 2020 research and innovation program funds T. B. under grant agreement 716542. The LABEX Lyon Institute of Origins (LIO, ANR-10-LABX-0066) of the University of Lyon funded a beowulf cluster hosted and maintained at ENSL, and used in this study. The world map figures were created with open software GMT 4.5.13.

Data availability:

The dataset generated during this study (3D Vs, Qs models and melt fraction models) is available as IRIS data products at https://doi.org/10.17611/dp/emc.2020.dbrdnature.1.

Code availability: Numerical modelling codes related with this paper can be downloaded from https://doi.org/10.17611/dp/emc.2020.dbrdnature.1. Requests about the numerical modelling codes associated with this paper should be sent to