Vladimir Lyakhovsky
email: vladimir.lyakhovsky@gmail.com

Eyal Shalev

Ittai Kurzon

Wenlu Zhu

Laurent Montesi

Nikolai M Shapiro

Effective seismic wave velocities and attenuation in partially molten rocks

Keywords: partial melt, phase transition, wave velocity, attenuation 2

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Melt fraction is largely present in the crust and uppermost mantle of the Earth. The partially molten rocks are characterized by enhanced mobility because of their elevated buoyancy and reduced viscosity and, therefore, play a key role in many geodynamical processes. Partial melting of silicate rocks occur in regions within the crust and in the uppermost mantle where the solidusliquidus transition can be reached by different mechanisms such as: decompression (e.g., [START_REF] Niu | Mantle Melting and Melt Extraction Processes beneath Ocean Ridges: Evidence from Abyssal Peridotites[END_REF], heating by uprising hot material such as mantle plumes (e.g., [START_REF] Campbell | Large Igneous Provinces and the Mantle Plume Hypothesis[END_REF], compositional changes such as the addition of water (e.g., [START_REF] Asimow | The importance of water to oceanic mantle melting regimes[END_REF], or by mechanical processes such as frictional melting (e.g., [START_REF] Toro | Natural and experimental evidence of melt lubrication of faults during earthquakes[END_REF]. The most direct manifestations of the presence of the molten rocks are the volcanic deposits observed at the surface, and the magmatism imprinted in geological records. Initial magmas form at depths and driven by buoyancy, rise toward the surface, through the uppermost mantle and the crust. While volcanic deposits of molten rock reach the surface, most of the magmas and molten rock remain stored at depth. Understanding how the magma and molten rocks are stored at depth is one of the major challenges in volcanology and igneous geology. Traditional views, implying the existence of concentrated "high-melt" volumes, forming well identified intrusions or magma chambers, is challenged today with new concepts suggesting a more distributed storage of melts in a form of low melt-fraction "mush reservoirs" [START_REF] Jackson | Chemical differentiation, cold storage and remobilization of magma in the Earth's crust[END_REF]. The form of such storage at depth strongly affects the evolution of magmas/melts and their mobility, changing the hydraulic connectivity of volcano-magmatic systems [START_REF] Pinel | On the relationship between cycles of eruptive activity and growth of a volcanic edifice[END_REF]. Melt distribution within the shallow layers of the crust is of great importance for estimating the potential volcanic activity. Estimation of the volume of erupted volcanic products from geological records shows that past rates of volcanism were very variable both for individual volcanic systems and for the whole Earth. In particular, the level of volcanic activity that occurred during human history is rather weak compared to some pre-human major events such as explosive super-Plinian eruptions, the emergence of flood basalts and the generation of large igneous provinces (e.g., [START_REF] Bryan | The largest volcanic eruptions on Earth[END_REF]; such strong events are often called upon to explain environmental changes and mass extinctions [START_REF] Sobolev | Linking mantle plumes, large igneous provinces and environmental catastrophes[END_REF]. The possibility of occurrence of new major volcanic episodes with great impact at a global scale cannot be excluded and its understanding requires the estimation of melt distribution within the shallow layers of the crust.

The presence of melt is also of major importance for global mantle dynamics and plate tectonics. Mantle convection in the vicinity of the major convergent and divergent plate boundaries is strongly influenced by the melt fraction, present in the subduction zone mantle wedges (e.g., [START_REF] Gerya | Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts[END_REF] and beneath the mid-oceanic ridges (e.g., [START_REF] Niu | Mantle Melting and Melt Extraction Processes beneath Ocean Ridges: Evidence from Abyssal Peridotites[END_REF]. The partial melt is also suggested to contribute significantly to the mechanical weakness of the asthenosphere (e.g., [START_REF] Chantel | Experimental evidence supports mantle partial melting in the asthenosphere[END_REF] and to facilitate the widespread small-scale sub-lithospheric convection(e.g., [START_REF] Van Hunen | New evidence for dislocation creep from 3-D geodynamic modeling of the Pacific upper mantle structure[END_REF]. Partial melt might be present in the continental crust driving fast deformation of the rocks, and affecting the orogenesis [START_REF] Chen | Melting conditions in the modern Tibetan crust since the Miocene[END_REF]. Also, partial melt can lubricate seismogenic faults [START_REF] Toro | Natural and experimental evidence of melt lubrication of faults during earthquakes[END_REF] and, therefore, affect earthquake dynamic and seismo-tectonic cycles.

Understanding many ongoing key geodynamic processes and related natural hazards requires information on the distribution of partial melt storage systems and reservoirs in the Earth's crust and upper mantle. This information can be obtained only by methods of deep geophysical sounding, searching for the areas with abnormal physical and mechanical properties. One possible approach is finding the presence of partial melt in the subsurface by searching for elevated electrical conductivity [START_REF] Watanabe | The relationship between electrical conductivity and melt fraction in a partially molten simple system: Archie's law behavior[END_REF]. So far, large scale electromagnetic and magneto-telluric surveys were used to study partially molten rocks located at depths beneath mid oceanic ridges (e.g., [START_REF] Baba | Mantle dynamics beneath the East Pacific Rise at 17 S: Insights from the Mantle Electromagnetic and Tomography (MELT) experiment[END_REF] and volcanoes (e.g., [START_REF] Hill | Distribution of melt beneath Mount St Helens and Mount Adams inferred from magnetotelluric data[END_REF].

Another approach is to use the significant reduction in mechanical properties, i.e., elastic moduli and seismic wave velocities, as well as enhanced anelastic attenuation in the areas with partially molten rocks that can be measured with different seismological methods [START_REF] Debayle | Seismic evidence for partial melt below tectonic plates[END_REF]. This approach is much widely used because of the existence of numerous permanent and temporary seismic networks. In highly active volcanic systems the mobility of melt can lead to recordable seismic radiation [START_REF] Chouet | Long-period volcano seismicity: its source and use in eruption forecasting[END_REF][START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF]. However, in most cases, the detection of molten rocks requires using different methods of seismic imaging. Since it is widely accepted that the shear modulus and S-wave velocity are most affected by the presence of melt (or more generally fluid fraction), the S-wave seismic tomography is particularly suitable to image the regions with presence of significant portions of melt such as in mid-oceanic ridges (e.g., [START_REF] Conder | Asthenospheric flow and asymmetry of the East Pacific Rise, MELT area[END_REF] hotspots [START_REF] Allen | Imaging the mantle beneath Iceland using integrated seismological techniques[END_REF]. Strong negative P-wave seismic velocity anomalies can also indicate the presence of melt, especially beneath volcanoes, and in many cases, the Vp/Vs ratio is used to trace the melts (e.g., [START_REF] Koulakov | Mantle and crustal sources of magmatic activity of Klyuchevskoy and surrounding volcanoes in Kamchatka inferred from earthquake tomography[END_REF]. Melt-induced variations of the elastic moduli also affect the reflectivity structure in the crust and in the uppermost mantle. Therefore, the melt-enriched regions can be imaged with the reflection seismic methods or with its passive seismology analog: the receiver functions (e.g., [START_REF] Rychert | Seismic imaging of melt in a displaced Hawaiian plume[END_REF]. Different methods of attenuation measurement using amplitudes of direct waves (P, S) and the coda of seismograms can be used in different geological settings (e.g., [START_REF] Prudencio | 3-D seismic attenuation structure of Long Valley caldera: looking for melt bodies in the shallow crust[END_REF]. The separation of the intrinsic attenuation and scattering is a particularly challenging issue. The latter can be strongly enhanced by the presence of small pockets of partially molten rocks.

As mentioned above [START_REF] Jackson | Chemical differentiation, cold storage and remobilization of magma in the Earth's crust[END_REF], magma and melt may be stored not within large intrusions or reservoirs but rather distributed over a system of small melt pockets or layers with dimensions smaller than seismic wavelengths. In such cases, the preferred orientation of these small partially molten volumes may result in seismic anisotropy. This mechanism can contribute to the well-known radial seismic anisotropy in the upper mantle [START_REF] Ekström | The unique anisotropy of the Pacific upper mantle[END_REF].

More recently, radial seismic anisotropy has been discovered beneath some large volcanic systems (e.g., [START_REF] Jaxybulatov | A large magmatic sill complex beneath the Toba caldera[END_REF] where it is believed to be caused by the presence of large sill complexes, composed of many nearly horizontal layers of magmatic rocks with different degree of partial melting.

Interpretation of the seismological observations requires understanding of the effects of the properties of the partially molten rocks on the melt fraction and topology of melt inclusions. Zhu et al. (submitted) approach this problem by numerical simulations of the effective elastic properties of different samples prepared in laboratory from a powder mixture of natural basalt with wide range of melt fraction values. They demonstrate that in-spite of complex topology of melt distribution the simple geometrical model of the elastic media with spherical fluid inclusions provides a good approximation for the mechanical behavior of the solid-fluid system. However, their static modelling, similar to most of previous studies, ignores the impact of heat and mass transport across the solid-melt interface associated with pressure oscillations. We utilize the simple geometry of spherical melt inclusions and develop analytical solution predicting the impact of the solid-melt phase transition on the bulk elasticity and wave attenuation.

In the next section we briefly review averaging models for effective properties of partially molten rocks. Later we provide key aspects of the existing mechanical models of partially molten rocks. Then we discuss the technique of sample preparation and the applied numerical method.

Section 5 describes results of numerical simulations and their comparison with existing models of effective elasticity of partially molten rocks (see Appendix A for model details). Then we summarize the main relations connecting bulk elasticity and wave attenuation with moving boundary between solid and melt phases. The detailed derivations are presented in Appendix B.

Our findings lead to practical suggestions that may improve the interpretation of seismological observations in areas with partially molten rocks.

Effective properties of Partially Molten Rocks

Bulk elasticity and seismic wave velocities

A partially molten rock is a composite material with solid and liquid (melt) fractions.

Effective elastic properties of composite materials have been studied for decades starting with original works by [START_REF] Voigt | Bestimmung der Elasticitätsconstanten des brasilianischen Turmalines[END_REF] and [START_REF] Reuss | Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle[END_REF] who suggested upper and lower bound values for any given volume fraction of constituents, ignoring the topology of the partially molten rock (see Appendix A for detailed description of these and other models). The Voigt (1890) iso-strain model assumes that all constituents have the same strain and provides the upper bound value, while the [START_REF] Reuss | Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle[END_REF] iso-stress assumption leads to the lower bound value. In the case of the solidmelt system with zero shear modulus of the liquid phase, lower bound value for the shear modulus is zero. Converting elastic moduli to V p and V s values allows estimating their expected reduction as a function of the melt fraction (Fig. 1). The expected reduction of the V p values is relatively minor, smaller than 20% velocity reduction for melt fractions of up to 25% (blue curves in Fig. 1).

However, the possible range of V s values is large (red curves in Fig. 1). The iso-strain, Voigt model predicts only minor reduction, whereas the iso-stress model predicts zero V s values. This wide range is not useful for practical applications and cannot be used to constrain the melt fraction.

More realistic bounds giving the narrowest possible range of bulk elasticity without specifying the geometries of the melt inclusions is the Hashin-Shtrikman-Walpole models. Several models consider sphere or ellipsoid inclusion shape and developed analytical solutions (equations a3-a5, Appendix A) connecting effective elastic properties with melt fraction. Several numerical methods such as self-consistent and differential effective medium methods were proposed to account for inclusion interaction and wide range of possible geometries. These methods as well as more accurate evaluation of the effective elastic properties achieved by characterization of the microstructure geometry are briefly discussed in the Appendix A.

Wave attenuation

Quantitative interpretation of seismic wave attenuation in partially molten rocks addresses two general mechanisms. One mechanism considers viscous dissipation due to melt flow or squirt flow in thin films (e.g., [START_REF] Mavko | The rock physics handbook Second Edition[END_REF][START_REF] O'connell | Viscoelastic properties of fluid-saturated cracked solids[END_REF]. Another mechanism addresses dislocation-diffusion stress relaxation in distributed melt pockets or inclusions (e.g., [START_REF] Karato | Defect microdynamics in minerals and solid-state mechanisms of seismic wave attenuation and velocity dispersion in the mantle[END_REF][START_REF] Morris | Diffusionally assisted grain-boundary sliding and viscoelasticity of polycrystals[END_REF]. Both approaches were examined in theoretical and laboratory studies, using natural rocks and synthetic materials, quantifying the seismic attenuation, and trying to establish a physical understanding of these mechanisms [START_REF] Faul | Transient Creep and Strain Energy Dissipation: An Experimental Perspective[END_REF] [START_REF] Yamauchi | Polycrystal anelasticity at near-solidus temperatures[END_REF]. Each approach has its shortcomings; while viscous dissipation, due to melt flow, requires unrealistically high melt viscosity [START_REF] Faul | Transient Creep and Strain Energy Dissipation: An Experimental Perspective[END_REF], the dislocation-diffusion mechanism does not explain the strong association of attenuation with the onset of melting [START_REF] Fontaine | Temperature dependence of shear wave attenuation in partially molten gabbronorite at seismic frequencies[END_REF].

The relaxation time for a melt-squirt within elliptical inclusions was derived by [START_REF] O'connell | Viscoelastic properties of fluid-saturated cracked solids[END_REF]. In the case when fluid fills narrow crack-like zones, the scaling relations should account for the crack aspect ratio or the ratio between opening and length (c/a). They concluded that the most efficient attenuation of fluid-saturated cracked solids is expected for a frequency range defined by two bounding frequencies:

3 1 ÷ ø ö ç è ae = a c K h w and ÷ ø ö ç è ae = a c h µ w 2 (1)
This frequency range, as well as the general time scale for the transient deformation, is associated with the stress relaxation process characterized by the Maxwell relaxation time defined as the ratio between fluid viscosity and elastic modulus of the rock mass. Viscosity of dacites and rhyolitic melts, that are expected to be present in the magmatic systems of the continental crust, is relatively high. The viscosity of rhyolite may be as high as 10 11 Pa s for dry melts and decrease to about 10 6

Pa s depending on water content. The Maxwell relaxation time for these melts is of the order of 10 -3 -10 +2 s, overlapping the seismological frequency range. Therefore the (O'Connell and Budiansky, 1977) analysis of fluid-filled cracks could well encompass the seismological observations. However, for typical viscosity range of basaltic melts 10 0 -10 2 Pa s and elastic modulus of the order of tens GPa, the Maxwell relaxation time is of the order of 10 -7 -10 -9 s, and the effects associated with visco-elastic relaxation are negligibly small after very short time intervals. Accordingly, we may expect significant attenuation only at very high frequency waves.

More accurate estimation (eq. 1) predicts that even for a small aspect ratio c/a ~ 10 -2 , the typical frequency range is above 1 KHz, which is well above the seismological frequency range. This was already noted by [START_REF] Mavko | Velocity and attenuation in partially molten rocks[END_REF][START_REF] Mavko | The rock physics handbook Second Edition[END_REF], who also discussed squirt flow mechanism.

These estimations are compatible with results discussed by [START_REF] Faul | Shear wave attenuation and dispersion in meltbearing olivine polycrystals: 2. Microstructural interpretation and seismological implications[END_REF] who noted that experimentally observed attenuation in melt-bearing olivine polycrystals reported by Jackson et al. ( 2004) is not compatible with melt squirt flow through an observed small pipe with aspect ratios around c/a ~ 10 -1 . They also noted that this flow mechanism requires unrealistically high melt viscosity (~10 4 Pa s). [START_REF] Faul | Shear wave attenuation and dispersion in meltbearing olivine polycrystals: 2. Microstructural interpretation and seismological implications[END_REF] suggested an alternative stress relaxation mechanism introducing effective viscosity for olivine-olivine grain boundary regions. Adjusting high "grainboundary viscosities" (~10 6 Pa s) allow fitting the observed attenuation; however, the background physical mechanism, responsible for temperature-dependent attenuation, is still not clear.

In addition, most of the recent works discuss different mechanisms for shear attenuation and neglect bulk attenuation. [START_REF] Vaišnys | Propagation of acoustic waves through a system undergoing phase transformations[END_REF] and [START_REF] Mavko | Velocity and attenuation in partially molten rocks[END_REF] discussed the mechanism of bulk attenuation and noted the role of the phase transition between solid and melt phases due to waveinduced pressure oscillations. [START_REF] Li | Effect of phase transitions on compressional-wave velocities in the Earth's mantle[END_REF] reported results of a complicated laboratory measurement demonstrating effects associated with the presence of phase changes. Their results support the velocity decrease and bulk attenuation observed during loading with a period of 1,000 s and the softening observed at low amplitudes and shorter periods. [START_REF] Durand | Constraining the kinetics of mantle phase changes with seismic data[END_REF]. However, this attenuation mechanism is not very efficient for body waves with typical frequency range 0.1 -1 Hz. In the case of the solid-melt phase transition, the characteristic time scale of the of the heat transport generated at the solid-melt interface (latent heat) varies between milliseconds to tens of seconds for melt inclusions of the sizes ranging from tens microns to centimeter scale [START_REF] Mavko | Velocity and attenuation in partially molten rocks[END_REF]. Because latent heat is very efficient in releasing energy, it may cause significant attenuation. In this study we develop the analytical solution predicting the impact of the heat generated at the solid-melt phase transition on the bulk elasticity and wave attenuation, considering infinitely fast kinetics.

Method

3.1 Bulk elasticity, no mass/heat transfer. [START_REF] Zhu | Microtomography of partially molten rocks: three-dimensional melt distribution in mantle peridotite[END_REF] obtained high-resolution 3-D determinations of melt distribution in experimentally produced olivine-basalt aggregates using X-ray synchrotron microtomography.

Olivine-basalt aggregates with melt fractions of 0.02, 0.05, 0.10, and 0.20 were synthesized from powder mixtures of pulverized natural forsterite and MORB. The spatial resolution of the synchrotron microtomographic images is ~0.7 microns. Comparison to the conventional 2-D scanning electron microscopy (SEM) images of the same samples demonstrated that the microtomographic method captured realistic melt distribution in 3-D over a large representative volume [START_REF] Zhu | Microtomography of partially molten rocks: three-dimensional melt distribution in mantle peridotite[END_REF]. More information regarding the experimental charges and the image analysis could be found in [START_REF] Miller | Experimental quantification of permeability of partially molten mantle rock[END_REF]. (schematically shown in Fig. 3) predicts surprisingly well the elastic properties of partially molten rocks, and will be used in the next section to study impact of heat and mass transport across the solid-melt interface.

Analytical solution of solid-melt phase transition, bulk elasticity and wave attenuation.

A small wave-induced pressure change at a solid-melt interface disturbs the thermodynamic equilibrium and leads to certain motion of the interface. We develop a model for cyclic load with a given amplitude and frequency keeping in mind effects associated with propagation of monochromatic waves. Transient effects related to a final duration of the wave package and temporal variations of the oscillation amplitudes are not considered here. Pressure and deviatoric stress in the solid phase and pressure in the melt, are coupled by stress continuity at the solid-melt interface during external wave-induced load and fluid pressure oscillations. For typical wave amplitude 10 -5 -10 -7 m and wavelength between 0.5-5 km, the deformations of about 10 -8 -10 -10 disturbs the thermodynamic equilibrium. Mathematical formulation of this moving boundary or Stefan problem includes heat, mass, and force balance equations. Any pressure change takes the system out of the thermodynamic equilibrium and causes solidification or melting which change the heat balance according to the Clausius-Clapeyron equation. The latent heat (sink or source) is transferred away or towards the interface by conductive-advective mechanism heating/cooling the entire rock mass leading to the energy loss and dissipation of the mechanical energy or seismic wave attenuation. An analytical solution is derived for stress distribution and heat exchange between melt inclusion and solid shell. We neglect the spatial temperature distribution in small melt inclusion. In larger inclusions, effective mixing of the low viscosity fluid probably eliminates any temperature gradients. Therefore we solve for the heat transport only in the solid phase. The analytical solution for the heat balance equation, including latent heat associated with the motion of the solid-melt interface, as well as temperature variations of the melt inclusion, gives the relation between pressure and volumetric strain oscillations. Finally, the derived compressibility of the rock and melt assemblage depends on the characteristics of the external wave-induced forcing. Wave attenuation, or quality factor (Q) is calculated from the time delay between pressure and strain oscillations, or the ratio between real and imaginary bulk moduli.

Results

Solid-melt phase transition, bulk elasticity and wave attenuation.

Numerical model simulating the effective elastic properties and results presented in Fig. 2 are obtained assuming mass conservation of melt and solid phases without any mass exchange associated with solid-melt phase transition. However, even small wave-induced pressure change at a solid-melt interface disturbs the thermodynamic equilibrium and leads to certain motion of the interface. Pressure changes move the system out of thermodynamic equilibrium, causing solidification or melting which change the heat balance according to the Clausius-Clapeyron equation. The latent heat (sink or source) is transferred away or towards the interface by conductive-advective mechanism, heating or cooling the entire rock mass, leading to energy loss and dissipation of the mechanical energy and to seismic wave attenuation. The presented model is based on the simplified cell geometry representing the partially molten rock as a set of equal elastic cells with spherical melt inclusions (Fig. 3). The simplified cell geometry allows deriving analytical relations for stress distribution and heat exchange between melt inclusion and solid shell considering three different regimes. Each cell representing the solid-melt system is at initial pressure and temperature (P 0 , T 0 ), corresponding to the thermodynamic equilibrium according to the phase diagram. The propagating pressure wave, with frequency 𝜔 and long wavelength (well above the cell size) may be considered as time-dependent pressure perturbations at the outer cell boundary, 𝑃 # 𝑡 = 𝛿𝑃 # 𝑒 )*+ (Fig. 3). The elastically deformed solid part of the cell transforms this pressure to the melt inclusion with a certain pressure factor 𝜓, i.e., melt pressure changes proportionally to the pressure wave amplitude, 𝛿𝑃 -𝑒 )*+ = 𝜓 • 𝛿𝑃 # 𝑒 )*+ , (see Appendix B for complete mathematical derivations). This pressure change violates the thermodynamic equilibrium, which is re-established by the motion of the solid-melt interface. According to the slope of the phase diagram (Clausius-Clapeyron slope, a) the temperature changes proportionally to the pressure, which relates the latent heat (heat of phase transformation) to the rate of pressure change with temperature (see table 1 for notations):

a r r = D × × = T L dT dP 2 (2)
Neglecting the time required to re-establish the thermodynamic equilibrium, the temperature variations at the solid-melt interface are:

𝑇 -𝑡 -𝑇 1 = 𝜓 𝛿𝑃 # 𝑒 )*+ /𝛼 (3)
Heat balance associated with moving solid-melt interface includes latent heat and heating/cooling of the melt inclusion and surrounding solid. Heat variation of small melt inclusion is proportional to its mass and temperature variation, while heat flux in the solid part of the cell is proportional to the temperature gradient at the interface. In Appendix B we provide exact solution for the temperature distribution in the solid shell with the oscillating boundary condition (3). The obtained solution, in the form of 𝑇 𝑟, 𝑡 = 6 7 𝑒 )87 𝑒 )*+ (4) allows calculating the diffusive heat flux towards the solid shell surrounding the melt inclusion:

𝑑𝐻 = 4𝜋𝜅𝜓 >? @ A 𝑅 𝑒 )*+ 𝐴 D + 𝑖𝐴 G 𝑑𝑡 (5)
where

𝐴 D = 1 + 𝜔𝑅 G 2𝐷 -𝐶𝑜𝑠 𝜔 2𝐷 (2𝑆 -𝑅) 1 + 2𝜔𝑆 G 𝐷 + 𝑆𝑖𝑛 𝜔 2𝐷 (2𝑆 -𝑅) 2𝜔𝑆 G 𝐷 𝑒 R * GS (G#RT) (6) 𝐴 G = 𝜔𝑅 G 2𝐷 + 𝑆𝑖𝑛 𝜔 2𝐷 (2𝑆 -𝑅) 1 + 2𝜔𝑆 G 𝐷 -𝐶𝑜𝑠 𝜔 2𝐷 (2𝑆 -𝑅) 2𝜔𝑆 G 𝐷 𝑒 R * GS (G#RT)
Solving the heat balance equation, which also includes latent heat associated with the motion of the solid-melt interface, as well as temperature variations of the melt inclusion, gives the relation between pressure and volumetric strain oscillations. Their ratio (𝐾 = 𝐾 T + 𝑖𝐾 V = -𝛿𝑃 # /𝛿𝜀 X ) is the complex bulk modulus:

𝐾 = D Y Z + [\] A^ _ @ 𝜌 -𝐶 ? + 𝐴 G -𝑖𝐴 D ab *T c RD (7)
The ratio K R /K I between the real and imaginary parts of the bulk modulus (7) defines the bulk quality factor, Q, of the media (see Appendix B)

𝑄 = A^_ @ *T c e[\] Y Z _ f g h *T c e[\]Y Z 6 c ab Y Z [\]6 i ab (8)
These derivations demonstrate that moving solid-melt phase boundary due to waveinduced pressure oscillations leads to additional reduction of the bulk modulus, K e , estimated from the elastic deformation and provides efficient mechanism for wave attenuation. Most of the material properties in (6, 7, 8) are well defined with relatively small uncertainty range (e.g., [START_REF] Ricard | 02 -Physics of Mantle Convection[END_REF]. The exception is the product of the Clausius-Clapeyron slope and latent heat values, which according to (2) is expressed as:

𝛼𝐿 = A c k ∆_ _ c (9)
Therefore, we ignore the uncertainty in most of the parameters presented in Table 1, and study the sensitivity of the model predictions to the reasonable variations of the Clausius-Clapeyron slope values 𝛼 = 3, 6, 9 MPa/ o K [START_REF] Bina | Phase transition Clapeyron slopes and transition zone seismic discontinuity topography[END_REF]. Larger slope values, from ~7.5 MPa/ o K at low pressures up to ~30 MPa/ o K at 10 GPa were reported by [START_REF] Hirschmann | Mantle solidus: Experimental constraints and the effects of peridotite composition[END_REF] for various compositions of mantle Peridotites.

Fig. 4 shows reduction of the bulk modulus as a function of the melt fraction calculated using eq. 7. We define the cell-size, as the typical distance between melt pockets. In brittle rheology, these melt pockets are typically grain-size away from each other, while for ductile rheology they would be further separated. For small cell-size (Fig. 4a for S 0 =0.5 mm) the reduction is larger than for relatively large cell-size (Fig. 4b for S 0 =5.0 mm). In both cases smaller Clausius-Clapeyron slope values lead to more significant effect, since the amplitude of the temperature oscillations at the solid-melt interface (3) is proportional to 1/a. The most pronounced change above 80% reduction occurs at low melt fraction values below 2% for 𝛼 = 3 and 6 MPa/ o K (Fig. 4a) and then remains about constant. Larger cell-size (Fig. 4b) causes gradual reduction of the bulk modulus. Values of the quality factor (Q) for 1 Hz p-wave frequency, calculated using eq. 8, also strongly depend on these factors (Fig. 5). For small grain size (S 0 =0.5 mm) the reduction is practically abrupt at very low melt fraction values and then Q remains practically constant within the values shown by the red polygon in Fig. 5. For larger S 0 =5.0 mm, the Q decrease is steeper for small a-values and approach about constant value at melt fraction values above ~5%. Since the presented model ignores wave attenuation in the solid phase, the Q-values go to infinity for zero melt fraction value (b=0).

Wave attenuation predicted by the model, is strongly frequency-dependent (Fig. 6). Less attenuation or higher Q-values are obtained at very low and high frequency values. If the oscillation period is very large, the whole cell volume is heated to the equilibrium temperature and the delay associated with heat transfer is small. Similarly, heat transfer can be neglected at very small period values. The highest attenuation occurs when the period of oscillation is close to the heat transfer time scale. In this condition, any heating that occurs due to Solid-melt phase transition during the oscillation is affecting the entire cell and therefore causing high attenuation. The time needed to transfer heat over the cell size or the characteristic time scale of conductive heat transfer over the cell size is defined as 𝑆 1 G 𝐷. Therefore, the largest attenuation occurs when 𝑓 • 𝑆 1 G /𝐷 ~ 1. For small cell-size (Fog. 6a) maximum attenuation is expected at frequency range, f, between 1 and 10 Hz.

Q-values are very low and weakly depend on the melt fraction. One order of magnitude increase in the cell-size causes two orders of magnitude (𝑆 1 G ) frequency change (0.01 -0.1 Hz) for maximum attenuation (Fig. 6b). In both cases, the increase in Q-values at low frequencies only weakly depends on the melt fraction, while it is significantly more sensitive at high frequencies.

Discussion

Understanding how magma or molten rocks are stored at depth is one of the major challenges in volcanology, igneous geology, mantle dynamics, and plate tectonics. Many of the ongoing geodynamic key processes and related natural hazards require information on the distribution of partial melt storage system and reservoirs in the Earth's crust and upper mantle.

Seismological methods are often used to show areas of significant reduction in seismic wave velocities and enhanced wave attenuation, and are interpreted as revealing pockets of partially molten rocks. Interpretation of the seismological observations requires understanding of the effects of the properties of the partially molten rocks on the melt fraction and on the topology of melt inclusions.

Previous studies show that viscous dissipation due to flow or squirt of the melt with realistic viscosity range is small and expected attenuation is well below the typical values for the partially molten rocks. The developed analytical model shows that phase transition is an efficient mechanism for strong attenuation and for the reduction of wave velocity. The model adopts cell geometry for spherical melt inclusions in partially melted rocks (Fig. 3). Previous results discussed by Zhu et al. (submitted) and shown in Fig. 2 demonstrate that in spite of a complex melt topology, this simple geometrical model provides surprisingly good approximation for the effective rock elasticity. This important result allows using analytical solutions, assuming spherical inclusions, for estimating effective elastic properties and modeling wave attenuation.

The most realistic bounds giving the narrowest possible range of bulk elasticity leaving solid-melt phase transition out of the consideration is the Hashin-Shtrikman-Walpole model.

Their model predicts relatively minor reduction in both pressure and shear wave velocities with melt fraction (Fig. 7). However, even for relatively large cell-size (S 0 =5 mm) and intermediate Clausius-Clapeyron slope value, 𝛼 = 5 MPa/ o K, the reduction in the bulk modulus associated with the phase transition is dominant (Fig. 7). The analytical solution considers only spherical symmetry and volumetric deformation, and formally cannot relate phase transition to shearing, which in perfect material is decoupled from volumetric deformation. However, original studies by [START_REF] Reynolds | On the dilatancy of media composed of rigid particles in contact. With experimental illustrations[END_REF] and [START_REF] Skempton | The pore-pressure coefficients A and B[END_REF], have shown that shear stress is strongly coupled with pressure, and therefore the formal decoupling is in doubt. Taking into account the relation between S-wave induced shearing and pressure, we speculate here that propagating S-waves will also lead to similar phase boundary motion, i.e., velocity reduction and attenuation. Dashed red line in Fig. 7 shows the inferred shear modulus reduction using the same factor as for the P-waves. More accurate relations for the S-waves will be examined in future study.

In summary, the analyzed mechanism of reduction in mechanical properties associated with the phase transition predicts that the presence of significant amount of melt in the rocks would result in a strong decrease in seismic wave velocities associated with a strong increase of their attenuation. Preliminary consideration also indicates that the velocity decrease should be more enhanced for S wave (Fig. 7), resulting in the increase of the Vp/Vs ratio. This implies that simultaneous measurement of these properties could help to better constrain the melt content in studied geological objects such as, for example, volcanic systems. A key parameter for interpretation of such results is the average cell-size S 0 , that remains poorly known. Results shown in Fig. 6 indicate that this parameter cannot be too small (< 1mm) because in this case the waves at frequencies between 1 and 10 Hz (i.e. typical range of observation for volcanic earthquakes) would be completely dissipated (Q< 5). Considering the millimeter or centimeter scale of the melt distribution spatial heterogeneity would result in more appropriate predictions with Q varying between 10 and several hundreds at 1 Hz. Moreover, with such values of S 0 , the model predicts an approximately power-law increase of Q at high frequencies, i.e., the behavior reported from field observations and laboratory experiments [START_REF] Lekić | Measurement and implications of frequency dependence of attenuation[END_REF]. Overall, the predicted values are approximately in agreement with large scale seismological observations showing that seismic waves are mostly attenuated within the shallow parts of the Earth's crust and mantle, and associated with possible presence of melt, such as in mid-oceanic ridges and island arcs or beneath large volcanic systems. The model is also in agreement with observed significant correlation between the regional-scale tectonics and the geographical distributions of uppermost-mantle velocities and attenuation for S and P waves. So far, average Q being as low as ~60 in the upper mantle low velocity zone, where significant amount of melt may be expected (e.g., [START_REF] Thybo | The heterogeneous upper mantle low velocity zone[END_REF].

The presented model considers latent heat generation at the solid-melt phase transition and heat transport near the interface. Adopting infinitely fast kinetics, the model predicts the realistic body wave attenuation with the frequency range 0.1-1 Hz. The alternative approach considering solid-solid transition with finite phase [START_REF] Durand | Constraining the kinetics of mantle phase changes with seismic data[END_REF][START_REF] Ricard | Seismic attenuation in a phase change coexistence loop[END_REF], predicts maximum attenuation at significantly lower frequency range typical for free oscillations. The real case might be somewhat in between and the complete model should consider both type of the phase transitions.

Conclusion

In spite of the complex structure of partially molten rocks, effective elastic moduli can be approximated assuming spherical melt inclusions geometry, which implies that the exact topology of the melt even at high melt fraction values is not crucial. Viscous flow mechanisms for basaltic melts can cause attenuation only at very high frequency waves. Explaining observed wave attenuation in partially molten peridotite rocks using viscous flow requires adjusting unrealistically high viscosities.

We present a new mechanism for wave velocity reduction and attenuation associated with moving solid-melt interface (Stefan problem), and latent heat generation due to wave-induced pressure oscillations around thermodynamic equilibrium. The reduction in mechanical properties (seismic velocities) associated with the phase transition is shown to be significant. The highest attenuation occurs when the period of oscillation is close to the heat transfer time scale. In this condition, any heating that occurs due to Solid-melt phase transition during the oscillation is affecting the entire cell and therefore causes high attenuation. Velocity reduction (%)

  Using the reported 3-D melt geometries,Zhu et al. (submitted) recently conducted finite element deformation experiments to obtain the effective bulk and shear moduli of these experimental charges. In the numerical simulations, the numerical elements were dictated by the pixel size of the tomography image with every pixel represented by a cubic brick separated into five tetrahedron elements. The grid avoids any smoothing or averaging of the observed microstructure. The material properties of each pixel (set of five tetrahedron elements) are defined as solid or melt according to the tomographic image. Following[START_REF] Hammond | Upper mantle seismic wave velocity: Effects of realistic partial melt geometries[END_REF], we use K S = 124 GPa, K M = 64 GPa, and µ S = 40 GPa (solid Poisson ratio is 0.28). Infinite viscosity (zero fluidity) is defined for the solid phase, while typical viscosity of 10 Pa s and zero shear modulus are adopted for the melt phase. The viscous strain component act to dampen deformation and facilitates the numerical transient phase that goes from the initial un-deformed equilibrium state to the final static deformed equilibrium one. Boundary conditions for every model run are set according to the simulated loading, i.e., 3-D compaction or uniaxial load. Static stress distribution is simulated using our 3-D code utilizing Fast Lagrangian Algorithm verified in many previous studies.The actual melt fraction, calculated for every sub-volume, depends on the amount of melt generated in the original laboratory sample and also varies between different sub-volumes of the same sample representing the heterogeneity and complex topology of the melt distribution. The lack of any clustering or any preferred location of markers in Fig.2corresponding to the size of the numerical grid (100, 200, 300, and 400 pixels) confirms that the obtained results are robust and there is no meaningful grid dependency. Colored lines in Fig.2show effective elastic moduli predicted by different analytical models discussed in Appendix A. Results of numerical simulations of the effective elastic properties demonstrate that the simple geometrical model of the elastic media with spherical melt inclusions provide the most reliable approximation for the mechanical behavior of the solid-fluid system, in spite of very simplified topology. This topology of the melt distribution in a form of spherical inclusion located in the center of equal cells

Figure 1 :Figure 2 :

 12 Figure1: Upper and lower bounds for velocity reduction according to the iso-strain (Upper) and iso-stress (Lower) models, for both, Vp (blue) and Vs (red) seismic velocities.

Figure 3 :

 3 Figure 3: The cell model for spherical melt inclusions in partially melted rocks (see text for further explanation). R and S are the radii of the inclusion and cell, respectively, 𝜎 77is the radial component of the stress in the cell, P0 is the equilibrium external pressure, under which the melt inclusion co-exists with the surrounding solid shell, and 𝛿𝑃 -and 𝛿𝑃 # are the pressure oscillations at the inner and outer solid cell boundaries.

Figure 6 :

 6 Figure 6: Q versus frequency for different melt fraction values b = 0.1%, 1.0%, and 10% and Clausius-Clapeyron slope value a = 5 MPa/ o K. Shaded polygons show Q-values for S0 = 0.5 mm, 5.0 mm, and 50.0 mm.

Figure 7 :

 7 Figure 7: Velocity change predicted by the model with spherical inclusions and corrected after the phase transition effect for S0 = 5.0 mm, a = 5 MPa/ o K, and 1 Hz frequency.

  [START_REF] Ricard | Seismic attenuation in a phase change coexistence loop[END_REF] and Durand et al. (2012) investigated seismic wave attenuation due solid-solid phase transformation. They considered finite rate of kinetics, but without latent heat generation. Their model successfully predicted seismic mode attenuation associated with phase transition at 410 km depth

Table 1 .

 1 Values of the model constants. 479

	Notation	Parameter	Value/units
	K S	Solid bulk modulus	124 GPa
	µ S	Solid shear modulus	64 GPa
	n S	Solid Poisson ratio	0.28
	K M	Melt bulk modulus	40 GPa
	K e	Effective bulk modulus	GPa
	µ e	Effective shear modulus	GPa
		b	Melt fraction	0-20%
	a	Clausius-Clapeyron slope	1-10 MPa/ o K
		k	Heat conductivity	3.5 J/s m o K
	C p	Heat capacity	1 kJ/kg o K
		L	Specific latent heat	10 2 -10 3 kJ/kg
	r S	Solid density	3,500 kg/m 3
	r M	Melt density	3,000 kg/m 3
	𝐷 =	𝜅 𝜌 # 𝐶 ?

Acknowledgements

We thank Claude Jaupart and Navot Murag for discussions. The paper benefited from useful comments by two anonymous referees and the editor, Hans Thybo. This study was supported by the Russian Ministry of Education and Science (Grant 14.W03.31.0033) and by the European Research Council under the European Union Horizon 2020 research and innovation program (Grant Agreement 787399-SEISMAZE).


Guillaume Quibeuf

François Charru

Laurent Lacaze

Particle motion within the viscous sublayer of a turbulent shear flow

The motion of isolated heavy particles on the bottom wall of a turbulent shear flow is investigated experimentally (using high frequency particle tracking), and correlated to that of the fluid (from high resolution LDV measurements). Numerous particle diameter, particle density and flow velocity are considered, where the particle remain within the viscous sublayer (particle Reynolds number lower than 10). It is shown that fluid flow fluctuations within the viscous sublayer (varying Couette flow) induce particle velocity fluctuations of same magnitude -slightly smaller -. The mean particle velocity, standard deviation and Lagrangian correlation time, are found to scale with the viscous scales. Particles roll and slide without takeoff, unlike what happens in viscous laminar flow where the lift force may overcome the immersed particle weight. Arguments are given for this difference.

Probability density functions for the streamwise and spanwise velocity are found to follow Gamma and Gaussian laws, respectively, which also scale with the viscous scales.

I. INTRODUCTION

The shearing action of a water flow above a granular bed leads, for moderate shear stress, to bedload transport where particles move close to the bed. This transport has been extensively studied at a macroscopic level with emphasis on the mean sediment discharge, because the latter is the most relevant quantity for most applications (river morphodynamics, hydraulic engineering. . . ). In this context, there exists a great variety of semi-empirical sediment transport laws derived from laboratory experiments [1]. However, their applicability range is quite narrow and discrepancy between predictions and in situ measurements can be large, especially for shear-stress close to threshold [START_REF] Nakato | Tests of selected sediment-transport formulas[END_REF][START_REF] Bathurst | Effect of coarse surface layer on bed-load transport[END_REF]. These laws put forward the dominant effect of the fluid shear stress, as measured by the Shields number which represents the ratio of the hydrodynamic force on a particle to its immersed weight. However, they generally do not account for the weaker effect of particle inertia (as measured by the particle Reynolds number), nor the effect of other phenomena more difficult to control and model such as the bed packing density or heterogeneities of the particles shape and size.

Deeper insights have been achieved from studies at the grain scale, using particle tracking for the investigation of local dynamics, in viscous flows [START_REF] Charru | Motion of a particle near a rough wall in a viscous shear flow[END_REF][START_REF] Mouilleron | Inside the moving layer of a sheared granular bed[END_REF] and turbulent streams [START_REF] Lajeunesse | Bed load transport in turbulent flow at grain scale : experiments and modeling[END_REF][START_REF] Shim | Experimental study of bed-load transport using motion tracking[END_REF]. For the latter, detailed and time-resolved measurements of the near-bed fluid flow have improved our understanding of the interactions between turbulent structures and grain motion [START_REF] Nelson | Role of the near-bed turbulence structure in bed load transport and bed form mechanics[END_REF][START_REF] Lelouvetel | Experimental investigation of ejections and sweeps involved in particle suspension[END_REF] and their feedback on the flow [START_REF] Gyr | Turbulent flows over smooth erodible sand beds in flume[END_REF][START_REF] Campbell | Bed-load effects on hydrodynamics of rough-bed open-channel flows[END_REF][START_REF] Dey | Near-bed turbulence characteristics at the entrainment threshold of sediment beds[END_REF].

The need of work on the motion of solitary grains -over a bed of sticked particles -appeared early, in order to get a deeper understanding of the hydrodynamic forces experienced by the grain, interactions with the bed (rebounds, momentum transfer, friction coefficient), and resulting trajectories. Quite different situations occur, from that of large grains (with Reynolds number of a few hundreds) saltating over the bed [START_REF] Francis | Experiments on the motion of solitary grains along the bed of a water stream[END_REF], to that of small particles reptating in viscous flows [START_REF] Charru | Erosion and deposition of particles on a bed sheared by a viscous flow[END_REF]. Other works were devoted to the problem of the flow conditions necessary to dislodge a grain out of the bed, showing in particular that the relevant quantity is the flow impulse (product of the instantaneous hydrodynamic force by its duration) rather than the force itself [START_REF] Celik | Impulse and particle dislodgement under turbulent flow conditions[END_REF], or, equivalently, some energy transfer [START_REF] Valyrakis | Entrainment of coarse particles in turbulent flows: An energy approach[END_REF]. In another direction, attention has been paid to the motion of particles over smooth walls, in relation with industrial problems such as the cleaning of surfaces, filtration, or field-flow fractionation. Available studies typically provide correlations for the particle velocity as a function of the friction velocity [START_REF] Sechet | Bursting phenomenon and incipient motion of solid particles in bed-load transport[END_REF][START_REF] Julien | Velocity of rolling bed load particles[END_REF][START_REF] Campagnol | Particle motion of bedload sediment moving over a smooth bed[END_REF].

This paper reports an experimental investigation of the motion of isolated particles (small surface density) carried by a turbulent flow over a smooth bottom wall, for small particles lying within the viscous sublayer. The objective is to clarify the link between the dynamics and statistics of the near-wall fluid flow and those of individual grains, thanks to precise and reliable velocity measurements of both phases, and find their relevant scales. The experimental setup is first presented (Sec. II). Then results are given for the mean velocity, standard deviation, probability density functions and correlation times, for both the flow and the particles (Sec. III). Results are summarized in Sec. IV.

II. EXPERIMENTAL SET-UP

The experimental set-up, as sketched in Fig. 1, consists of a 5.59 m long duct with rectangular cross section with inner height H = 59 mm and width B = 120 mm. It is basically the same as that described in Ref. [START_REF] Franklin | Subaqueous barchan dunes in turbulent shear flow. Part 1. Dune motion[END_REF] with improvements given in Ref. [START_REF] Quibeuf | Étude expérimentale du transport sédimentaire hors équilibre[END_REF].

The first 4.08 meters of the duct are made in Perspex, the end section (1.59 m long) where all experiments take place is made in a good quality glass in order to make observation easier. The duct is entirely filled with tap water (closed-conduit flow i.e. no free surface) and is pressurized by a 2 m high water column to avoid air bubble entrance. A divergenthoneycomb-convergent system homogenizes the turbulence at the entrance of the duct. The flow developpes along the Perspex part and one experimentally verified that it has reached a uniform regime (downstream invariance) at the end section. The flow rate is varied via a progressive cavity pump and monitored with an electromagnetic flow meter.

The duct Reynolds number Re = ρU deb D/η may reach 3 × 10 4 , where U deb is the mean fluid flow velocity, D = 2BH/(B + H) the hydraulic diameter, ρ = 10 3 kg m -3 the water density and η = ρν = 10 -3 Pa s its viscosity. The flow field was measured by Laser Doppler Velocimetry (Dantec Dynamics high power flowlite 2D system), in the vertical symmetry plane, at 4.2 m downstream of the entrance of the duct. The LDV probe, mounted on a 3-axis traverse system, permits to scan the flow along the vertical. Measurements were taken at every 0.05 mm within the first millimeter above the bottom wall, and at every millimeter above. The typical acquisition rate was about 100 Hz, with 5000 samples collected at each point. Let us decompose the streamwise fluid velocity as u = u + u where u(y) is the timeaveraged velocity and u the fluctuation, with similar notations for the other velocity components. The wall shear stress τ w , as determined from the near-wall slope of the mean velocity, τ w = ρν(∂u/∂y) y→0 , was found to agree with the Blasius correlation, τ w /( 1 2 ρU 2 deb ) = 0.0791 Re -1/4 [START_REF] Franklin | Subaqueous barchan dunes in turbulent shear flow. Part 1. Dune motion[END_REF]. This correlation was used to compute the wall scales: the shear velocity u τ = (τ w /ρ) 1/2 and the viscous length δ v = ν/u τ . Vertical profile of the mean velocity u + (y + ), in wall units, are displayed in Fig. 2 The particles were white spherical glass or zirconium beads from Sigmund Lindner. Their density and median diameter d, after sieving, are given in Table I together with the explored ranges of the flow Reynolds number Re, friction velocity u τ , particle Reynolds number Re p = u τ d/ν and Shields number θ = τ w /((ρ p -ρ)gd) (dimensionless wall shear stress). It can be seen that Re p , which represents the particle diameter in wall units, is always below 10, which means that the particles remain within the viscous sublayer (of conventional thickness 11.6 δ ν [START_REF] Pope | Turbulent Flows[END_REF]).

For a typical experiment, a few particles were released on the bottom wall of the duct (painted in black to enhance the contrast with the white particles), with surface density of a few particles per square centimeter. Their motion, as the flow was set out, was then followed by a CCD camera (Basler ace acA2040-180km, 2000 × 2000 pixels) placed above the last third of the duct length. The horizontal field of view was 18 × 18 mm 2 and the frame rate f was between 60 and 160 Hz. The images were then processed with a particle tracking algorithm1 providing the particle trajectories and velocities. Note that the particles were observed in the central third of the duct width, where the influence of the lateral walls is quite negligible.

III. RESULTS

A. Particle trajectory and mean velocity

A typical particle trajectory in the horizontal plane, and the corresponding streamwise and spanwise velocity components, U and W , are shown in Fig. 3 for Re p = d + = 4.8.

It can be seen that the bead has an erratic motion with periods of sudden accelerations and decelerations, up to 3u τ . This behavior, which was observed for all flow conditions, shows that although the particle lies within the viscous sublayer (d + < 11.6), its motion still reflects some unsteadiness of the turbulent flow.

U + = 0.26 Re p . (1)
In other words, the mean particle velocity is 0.26 γd, where γ = u 2 τ /ν is the viscous shear rate close to the wall, i.e. it is equal to the fluid velocity at the distance 0.26 d above the wall.

B. Assessment of a viscous flow model

Let us now compare the above particle velocity (1) with that found for isolated heavy particules in laminar viscous flow, by Leighton and coworkers [START_REF] Krishnan | Inertial lift on a moving sphere in contact with a plane wall in a shear flow[END_REF][START_REF] King | Measurement of the inertial lift on a moving sphere in contact with a plane wall in a shear flow[END_REF]. For their particles with small inertia, these authors identified three regimes. First, for low shear rates, the particle rolls on the wall with constant speed, without slipping. Then, slip occurs and the particle velocity increases. Eventually, for high shear rates, the lift force overcomes the immersed bead weight, the bead takes off and stabilizes at a finite (small) distance above the p /Re s ≈ 32. For higher shear, the model predicts takeoff and higher particle velocity, whereas our observations do not reveal anything special: the velocity still scales with the shear rate u 2 τ /ν (constant U + /Re p ). The reason for the absence of any takeoff transition may be explained as follows. The lift force, as calculated in Refs 24 and 25, is the sum of six terms, each arising as the result of a pairwise combination of the translational velocity, rotational velocity, and the shear rate. These terms involve a numerical coefficient which may be positive (for four of them) or negative (for the last two). In turbulent flow, where particles experience a wide range of sliding and rolling motions, and shear rates, it can be expected that the resulting mean lift-force be very small, with either positive or negative sign. Hence the absence of takeoff transition observed in Fig. 5 for our experiences, unlike what happens in laminar flow. This explanation also corroborates direct visual observation, which seem to indicate that particles always keep contact with the wall.

C. Standard deviations of the particle velocity and fluid velocity

It is well-known that heavy particles moving on a wall in turbulent flow may experience large velocity fluctuations, see e.g. Ref. 17. We explore here more thoroughly the statistics of these fluctuations and assess their correlation with those of the fluid. (dotted line in the figure). Thus, σ + U scales with the particle Reynolds number, as the mean velocity U + , and represents a large part of the latter, of about 60%. These results show that, although the Reynolds stresses are small within the viscous sublayer, large fluctuations may however exist there. (

) 3
Closer inspection shows that the r.m.s. velocity is slightly smaller for the heavy zirconium beads than for the lighter glass ones, which may be interpreted as a filtering effect due to particle inertia. The present fluid measurements agree with previous direct measurements of the wall shear-stress [START_REF] Keirsbulck | Statistical properties of wall shear stress fluctuations in turbulent channel flows[END_REF] which showed that fluctuations may represent 40% of the mean shear. These fluctuations likely correspond to the streamwise elongated streaks of alternating fast and slow fluid that meander in the near-wall region, below y + ≈ 100 [START_REF] Smith | The characteristics of low-speed streaks in the near wall region of a turbulent boundary layer[END_REF].

Now consider the fluctuations of the spanwise velocity, σ W . Distributions (not shown) were found to be Gaussian, consistently with previous findings for a particle on a rough bed in viscous flow [START_REF] Charru | Erosion and deposition of particles on a bed sheared by a viscous flow[END_REF]. Fig. 7 shows that, in wall units, σ + W scales with the particle Reynolds number Re p (i.e, σ W scales with the viscous velocity u 2 τ d/ν, as the streamwise component). However, a particle density effect here appears: for the same Re p , the r.m.s. velocity is much smaller for the heavy zirconium beads. In other words, heavy particles, with larger inertia, resist more strongly to spanwise fluid velocity fluctuations. A general correlation accounting for this inertia effect, also reported in Fig. 7, is

σ + W = 0.21 ρ ρ p Re p . (4)

D. Probability density functions

Now turn to the probability density functions (P DF ) of the fluid and particle streamwise velocities. Figure 8(a) displays, for five values of the Reynolds number, the normalized P DF of the dimensionless fluid velocity u + at the height y = 0.31 mm (which corresponds to the diameter of the beads G2). It can be seen that all the data points fall onto a single curve, which shows that all statistics within the viscous sublayer scale with the wall units.

P DF of the grain velocities were computed, for each flow velocity, from typically 650 trajectories corresponding to a number of about 50 × 10 3 samples. the particle motion scales with the viscous units with negligible inertia effect. However, for U + /Re p 0.5, which corresponds to velocity larger than twice the mean (U + /Re p = 0.26), the P DF separate, without, however, any clear trend with the Reynolds number variations.

The solid lines in Figs. 8(a) and 8(b) correspond to the Gamma distributions with P DF :

P DF (X) = 1 b a Γ(a) X a-1 e -X/b , (5)
where X stands either for u + or U + /Re p , Γ is the Gamma function, and a = X 2 /σ 2 X and b = σ 2 X /X are the two parameters of the distribution [START_REF] Pope | Turbulent Flows[END_REF]. For the fluid (Fig. 8(a)), fitting measurements with Eq. ( 5) provides a = 5.95 and b = 0.85. For the particles (Fig. 8(b)), the parameters are a = 3.0 and b = 0.086, as deduced from Eqs. (1) and [START_REF] Nakato | Tests of selected sediment-transport formulas[END_REF]. It is noticeable that both PDF, for the fluid and the grains, are well described by Gamma distributions. The positive skewness of these distributions, very common in near-wall turbulence, is considered as the signature of sweeps [START_REF] Zarić | Statistical analysis of wall turbulence phenomena[END_REF][START_REF] Durst | Probability density distribution in turbulent wall boundary-layer flows[END_REF][START_REF] Österlund | Experimental studies of zero pressure-gradient turbulent boundary layer flow[END_REF][START_REF] Boo | Near-wall hot-wire measurements. part II : Turbulence time scale, convective velocity and spectra in the viscous sublayer[END_REF]. Regarding the grains, Gamma laws have been previously evidenced in Ref. [START_REF] Campagnol | Particle motion of bedload sediment moving over a smooth bed[END_REF] (for shear higher than here) with long tails associated with the large velocities attained during small saltation jumps. Here, particles rather roll and slide in close contact with the wall (see the previous dicussion of Fig. 5). Such Gamma distributions of the particle velocity, with close resemblance with those for the fluid, show that particles follow the fluid fluctuations, consistently with the previous result that the mean and r.m.s. velocities scale with the viscous velocity γd = u 2 τ d/ν.

E. Correlation times

A last viscous argument may now be added to the above results, involving correlation times. Let us consider the autocorrelation function

R XX (τ ) = 1 σ 2 X X(t)X(t + τ ) (6)
where X stands either for the velocity fluctuations u or U , τ is the time lag and the overbar still denotes temporal average. The higher the duct Reynolds number is, the shorter the timescale for both the fluid and the particles.

More quantitatively, a correlation time τ c can be defined as the time lag where the autocorrelation has decreased to some conventional value, here of 0.2, corresponding, in Fig. 9, to the dotted lines (we have checked that choosing other values of the cut-off do not alter the conclusions to be drawn). Reynolds number, those of the beads being about one-half that of the fluid. No clear trend can be seen from one grain type to the other, as the grain density or diameter varies. The inset in Fig. 10 displays the same times once normalized with the viscous time γ -1 = ν/u 2 τ . The dimensionless times now appear quasi constant with the duct Reynolds number, as expected, with γτ c in the range 20-25 for the fluid and 10-15 for the grains. As a hint, the distance travelled by a grain with mean velocity 0.26 γd during the correlation time 15 γ -1 provides a correlation distance of about four diameters.

IV. SUMMARY AND CONCLUSION

The experimental study reported in this paper shows that in a turbulent boundary layer, the motion of heavy particles with Reynolds number Re p = u τ d/ν lower than 10 (i.e. particles within the viscous sublayer of characteristic thickness 11.6 ν/u τ ) closely follow fluid fluctuations (unsteady and space-varying Couette flow with mean shear rate γ = u 2 τ /ν). Consequently, all velocity statistics -mean velocity, standard deviation, PDF -scale with the viscous velocity γd = u 2 τ d/ν = u τ Re p which corresponds to the mean fluid velocity at the distance d above the wall. In particular, the mean streamwise particle velocity was found to be U ≈ 0.26 γd and its standard deviation σ U ≈ 0.15 γd. Accordingly, the correlation time was found to be 15 γ -1 , corresponding to a correlation distance of about four particle diameters.

Although the above scalings indicate small particle inertia, the spanwise motion was found to be affected by the particle density, with standard deviation proportional to the density ratio ρ/ρ p (i.e. smaller for heavier particles, as expected).

The fact that the particle velocity scales with the fluid velocity γd indicates that, within the explored range of flow conditions, particles do not take off and keep contact with the bed. This result, which is consistent with visual observations, was assessed by comparison with a theoretical model for a particle in (steady) laminar flow. Indeed, this model, validated with previous experiments, predicts that beyond some particle Reynolds number, the lift force on the particle exceeds its immersed weight. Then, the particle takes off, and its velocity departs from the previous viscous scaling, which phenomenon was not observed here. Then it can be argued from the structure of the expression of the lift force (see the discussion in the previous section) that in the unsteady and space-varying Couette flow taking place in the viscous sublayer, the positive and negative contributions to the lift force nearly cancel, so that takeoff cannot occur.

As a conclusion, it can be said that further investigations are now needed in order to understand how the above conclusions may be altered for larger particles emerging out of the viscous sublayer, on the one hand, and for particles over a rough wall, on the other hand.

[1] H. Chanson, Hydraulics of open channel flow: an introduction (Arnold, 1999).

FIG. 1 .

 1 FIG. 1. Sketch of the experimental set-up. (a) Progressive cavity pump; (b) electromagnetic flow meter; (c) 2 m-high water column; (d) divergent, honeycomb and convergent; (e) duct with rectangular cross section; (f) movie camera; (g) sedimentation tank.

FIG. 2 .

 2 FIG. 2. Streamwise fluid velocity in the vertical symmetry plane, in wall units, for Re = 18300 ( ), 21600 ( ), 25000 ( ). (a) Mean velocity; ( ), u + = y + ; ( ), u + = 1/κ log (y + ) + B with κ = 0.41 and B = 5.9. (b) Standard deviation σ u = u 2 1/2 ; ( ), Eq. (13) from Ref. 22.

  velocity increases linearly from the wall within the viscous sublayer (y + 10), and then logarithmically in the inertial range 40 y + 200. The standard deviation σ u = u 2 1/2 , displayed in Fig. 2(b), reaches a maximum σ + u ≈ 3 at y + ≈ 15 and then slowly decreases, in agreement with the semi-empirical relationship given by Ref. 22 (solid line).

  FIG. 3. (a) Chronophotography (from above, negative image) of a typical particle trajectory in the horizontal plane (the scale in the right bottom corner is 2 mm long); (b) and (c), corresponding streamwise and spanwise velocity components. Beads Z1 (see Table I), u τ = 14.6 mm s -1 , Re p = 4.8, f = 80 Hz.

Figure 4

 4 Figure 4 displays the mean particle velocity U + , in wall units, for the five bead types (open symbols), as a function of the particle Reynolds number Re p . It appears that all the data points fall onto the same curve, and that the particle velocity scales linearly with the Reynolds number, as

Figure 4

 4 Figure 4 also displays velocity measurements from Refs. 18 and 19, for Re p 20 -where particles emerge out of the viscous sublayer in the buffer layer -, and a curve fitting data. This curve corresponds to a semi-empirical law for the mean fluid velocity within the buffer layer, at the height y = d above the wall, given by Eq. (7) in Ref. 18. When extrapolated down to the viscous sublayer, this law provides u + ∼ 0.25 d + , which no longer is the fluid velocity at the height y = d (where it is u + = d + ) but fits, by chance, the present particle velocities.

Figure 6 ( 2 )FIG. 5 .

 625 FIG. 5. Mean particle velocity U + as a function of the shear parameter Re 4 p /Re s , and theoretical calculation (solid line) by Ref. 24 and 25 for relative roughness = 0.005 and friction coefficient µ f = 0.25. The vertical dashed lines locate the rolling-slipping and takeoff transitions predicted by the model. ( ), beads G1; ( ), beads G2; ( ), beads G3; ( ), beads Z1; ( ), beads Z2.

Figure 6 (

 6 Figure 6(b) compares σ U to the r.m.s. fluid velocity, σ u , at the height of the center of the particle, y = 1 2 d. It can be seen that σ U increases approximately linearly with σ u , and is always smaller, according to

FIG. 6 .

 6 FIG. 6. (a) r.m.s. fluctuations σ + U (filled symbols) and mean streamwise velocity U + (open symbols), in wall units, as a function of the particle Reynolds number Re p . Symbols: see Fig. 5; ( ), Eq. (1). ( ), Eq. (2). (b) σ U versus the r.m.s. fluid velocity σ u at the height y = 1 2 d.

Figure 8 (FIG. 8 .

 88 FIG. 7. Standard deviation of the spanwise velocity, σ +W , as a function of the particle Reynolds number Re p ; the lines correspond to Eq. (4). ( ), beads G1; ( ), beads G2; ( ), beads G3; ( ), beads Z1; ( ), beads Z2.

Figure 9 (

 9 a) displays, for five Reynolds numbers, the autocorrelation function of the fluid velocity at the height y = 0.31 mm (equal to the diameter of the beads G2). Similarly, Fig. 9(b) displays the autocorrelation function of the particle velocity, averaged over numerous trajectories. All plots in both figures display similar decrease with the time lag τ , with similar timescales a few hundreds of milliseconds.

FIG. 9 .FIG. 10 .

 910 FIG. 9. Autocorrelation functions of the streamwise r.m.s. velocity. (a), fluid at the height y = 0.31 mm; (b), glass beads with d = 0.31 mm. Insets: long time behavior. ( ) Re = 18300, (+) Re = 20000, ( ) Re = 21600, (×) Re = 23300, ( ) Re = 25000.

TABLE I .

 I Beads diameter and density, and range of the hydrodynamical conditions explored.

	Label	d (mm)	ρ p (kg m -3 )	Re	u τ (mm s -1 )	Re p	θ
	G1	0.13	2500	18300 -25000	13.5 -17.7	1.8 -2.3	0.095 -0.16
	G2	0.31	2500	18300 -25000	13.5 -17.7	4.2 -5.4	0.040 -0.069
	G3	0.54	2500	18300 -25000	13.5 -17.7	7.3 -9.6	0.023 -0.040
	Z1	0.31	3800	20000 -28300	14.6 -19.8	4.5 -6.1	0.025 -0.046
	Z2	0.54	3800	20000 -26600	14.6 -18.7	7.9 -10.1	0.014 -0.024

; ( ), Eq. (

13

) from Ref.

22.

  wall. The dimensionless parameter that naturally emerges from the model, and governs both transitions, is Re 4 p /Re s where Re s = V s d/ν is the settling Reynolds number (based on the Stokes settling velocity V s = (ρ p -ρ)gd 2 /18η). The transitions to slipping and takeoff occur, with our notations, for Re 4 p /Re s ≈ 1.6 and Re 4 p /Re s ≈ 32, respectively (corresponding, with the notations of Refs. 24 and 25, to Re 2 γ /Re s ≈ 0.2 and Re 2 γ /Re s ≈ 4). The velocities and transitions measured in Refs. 24 and 25 were found to agree with an analytical calculation of the forces and torques on the particle, involving two parameters: a typical roughness size d of the bead surface, 1, and a friction coefficient µ f . Figure 5 compares the present mean velocity measurements (data points) with the model of Refs. 24 and 25 (solid line), as a function of Re 4 p /Re s . It can be seen that the agreement is good up to Re 4

from D. Blair and E. Dufresne, freely available at http://site.physics.georgetown.edu/matlab/.


Emmanuelle Lavaine

Philippe Majerus

Nicolas Treich

Health, air pollution and animal agriculture 1

Introduction

Although animal agriculture is critical to the subsistence of smallholders in some poor countries, the global detrimental impact of animal farming is now both well documented and overwhelming.

Animal farming is a primary cause of deforestation [START_REF] De Sy | Land use patterns and related carbon losses following deforestation in South America[END_REF], biodiversity loss [START_REF] Machovina | Biodiversity conservation: The key is reducing meat consumption[END_REF], antibioresistance [START_REF] O'neill | Antimicrobials in agriculture and the environment: Reducing unnecessary use and waste. The Review on Antimicrobial Resistance[END_REF] and infectious diseases emergence and amplification [START_REF] Rohr | Emerging human infectious diseases and the links to global food production[END_REF]. Moreover, it contributes significantly to water pollution, water scarcity and climate change [START_REF] Godfray | Meat consumption, health and the environment[END_REF][START_REF] Poore | Reducing food's environmental impacts through producers and consumers[END_REF][START_REF] Springmann | Mitigation potential and global health impacts from emissions pricing of food commodities[END_REF].

Additionally, the exploitation of farmed animals, especially in its widespread intensive forms, raises various moral issues. In this paper, we discuss another impact of animal farming, that on air pollution and in turn on human health. While this impact is also potentially considerable, we stress that it has been largely overlooked by regulators as well as by researchers, and in particular by economists.

Pathways and Health Impacts

Air pollution is a major environmental health issue, and the necessity to reduce this pollution is well recognized both in academia and policy making. However, emissions of pollutants in agriculture have remained relatively stagnant over decades. In the EU, the decrease in particulate matter (PM) [START_REF] Nordgren | Pulmonary health effects of agriculture[END_REF][START_REF] Cambra-López | Airborne particulate matter from livestock production systems: A review of an air pollution problem[END_REF][START_REF] Oppliger | Exposure to bioaerosols in poultry houses at different stages of fattening; use of real-time PCR for airborne bacterial quantification[END_REF]. Chronic respiratory disorders such as asthma and ODTS (organic dust toxic syndrome) are also prevalent among poultry farm workers [START_REF] Viegas | Occupational exposure to poultry dust and effects on the respiratory system in workers[END_REF]. PM can also worsen symptoms of workers caused by preexisting lung disease [START_REF] Nordgren | Pulmonary health effects of agriculture[END_REF]. Evidence also suggests that ammonia exposure from animal farming causes irritation to the mucus membranes in the eyes and the respiratory system [START_REF] Kristensen | Ammonia and poultry welfare: a review[END_REF]. In addition, dust pollution in farming has been associated with headache, irritation of eyes, nose and throat, and drowsiness [START_REF] Hartung | Occupational and environmental risks caused by bio-aerosols in and from farm animal houses[END_REF]. Finally, livestock farm workers have significantly elevated mortality for several respiratory conditions, in particular for hypersensitivity pneumonitis [START_REF] Greskevitch | Respiratory disease in agricultural workers: mortality and morbidity statistics[END_REF][START_REF] Rautiainen | Mortality and morbidity in agriculture in the United States[END_REF]. However, much remains unknown about respiratory disease among livestock workers. In fact evaluation analyses are often limited to a single compound while exposure includes a combination of animal products, dust, pathogens, and chemicals [START_REF] Dignard | Recent research on occupational animal exposures and health risks: A narrative review[END_REF].

Individuals living in surrounding areas of intensive farms are also exposed to pollution from livestock production and investigations on such populations are discussed in recent literature reviews. Some have reported positive relationships between living close to intensive farms and asthma (mainly among children), wheezing, deficits in lung function, and higher incidence of pneumonia [START_REF] Casey | Industrial food animal production and community health[END_REF][START_REF] Douglas | A systematic review of the public health risks of bioaerosols from intensive farming[END_REF][START_REF] Smit | Impacts of intensive livestock production on human health in densely populated regions[END_REF]. However, a few studies actually indicate protective associations between intensive farm exposure and respiratory health [START_REF] Douglas | A systematic review of the public health risks of bioaerosols from intensive farming[END_REF] and one systematic review finds no evidence of a consistent association [START_REF] O'connor | Updated systematic review: associations between proximity to animal feeding operations and health of individuals in nearby communities[END_REF]. 4 5 One caveat in this literature is that the pathway between air pollution and health impacts remains unclear. For instance, the majority of studies use distance-to-farm as a measure for exposure rather than emission data, and outcomes are often based on self-reported health status, while mortality is not considered. Moreover, the fact that most of the reviewed studies are cross-sectional rather than longitudinal limits the ability to make strong causal inferences. Finally, the empirical validity in some of the publications may be questioned since important confounders are not controlled for. Air pollution is not fully exogenous and likely to be systematically related to the socioeconomic status of the household. For instance, wealthy households are less likely to live next to a polluting farm.

Hence, the non-random assignment of pollution prevents causal inferences. In particular, O'Connor et al. ( 2017) note that there is a high risk of bias among the studies they reviewed, and they cannot establish a consistent dose-response relationship between exposure and respiratory diseases.

Simulation-based air quality models

Air quality models are used to estimate ambient concentration levels of primary and secondary pollutants by sector, based on emission rates, meteorological conditions, and chemical processes in the atmosphere (EPA, 2016). In contrast to the epidemiological literature, these simulations are not restricted to local impacts, but account for the capability of particles to travel long distances and impact larger populations.

In the atmospheric chemistry literature, publications mainly rely on "full-form" models, which capture the complexities of such environmental processes by including detailed representations of each mechanism. In economics, usually reduced complexity versions of these models are applied in the form of integrated assessment models (IAM). They approximate estimates from full-scale models without explicitly representing all chemical and physical processes that impact pollutant fate and transport [START_REF] Baker | A database for evaluating the InMAP, APEEP, and EASIUR reduced complexity airquality modeling tools[END_REF]. Moreover, IAMs estimate the health damages as well as their economic cost, by integrating outcomes from epidemiological concentration-response functions and economic non-market valuation for morbidity and mortality in peer-reviewed literature (National Research Council, 2010).

Several recent publications relying on such air quality models identify agriculture to be one of the largest contributing sectors to premature mortality from PM2.5 concentration today, driven by ammonia emissions [START_REF] Brandt | Contribution from the ten major emission sectors in Europe and Denmark to the healthcost externalities of air pollution using the EVA model system-an integrated modelling approach[END_REF][START_REF] Lelieveld | The contribution of outdoor air pollution sources to premature mortality on a global scale[END_REF][START_REF] Tschofen | Fine particulate matter damages and value added in the US economy[END_REF]. Global simulations with full-scale models indicate that agriculture is the largest contributing sector to PM2.5 in eastern USA, Russia, East Asia and Europe. 6 Simulations show that a 50% decrease in agricultural ammonia emissions would reduce annual PM2.5 concentration levels by 11, 8 and 5% over Europe, North America and East Asia respectively [START_REF] Lelieveld | The contribution of outdoor air pollution sources to premature mortality on a global scale[END_REF][START_REF] Pozzer | Impact of agricultural emission reductions on fine-particulate matter and public health[END_REF]. As for economic investigations, IAM-based simulations show that, largely due to animal farming, agriculture has become the biggest contributor to PM2.5 damages in the US, with the sector's cost outweighing its value added to the economy [START_REF] Tschofen | Fine particulate matter damages and value added in the US economy[END_REF]. A prior assessment considering not only mortality but also morbidity indicates that the agricultural sector is the second largest contributor to health cost in Europe [START_REF] Brandt | Contribution from the ten major emission sectors in Europe and Denmark to the healthcost externalities of air pollution using the EVA model system-an integrated modelling approach[END_REF]. 7 Finally, [START_REF] Giannakis | Costs and benefits of agricultural ammonia emission abatement options for compliance with European air quality regulations[END_REF] compare health benefits of complying with the 2020 ammonia ceiling set by the EU with the cost of 6 Lelieveld et al. (2015) estimate that in many European countries agriculture's contribution to PM2.5 concentration is 40% or higher. 7 The study uses emissions from the year 2000 for the European level. The authors also estimate that in 2008 the Danish agricultural sector was the largest contributor to this country's health cost. abatement options to reach this goal. 8 They estimate an annual economic benefit of 14,837 M€ from avoided premature mortality, compared to abatement cost ranging between 80 and 3,738 M€ per year depending on the option, suggesting that more ambitious commitments could be set with relatively low cost. 9

Though the models applied in such investigations are highly sophisticated and able to account for complex atmospheric chemistry processes, the fact that they are simulations remains an important limitation. For instance, the produced concentration levels depend on the reliability of emissions data. While point sources are usually dependable due to mandated monitoring, ground-level sources may be problematic because they are more rarely measured and often estimated [START_REF] Tschofen | Fine particulate matter damages and value added in the US economy[END_REF]. Moreover, most of the reviewed models do not disentangle livestock-related emissions from other polluting agricultural activities. Finally, there is currently still a lack of knowledge regarding how toxicity of particles from different sources varies. For instance, agricultural emissions forming mostly inorganic particles may be much less harmful than carbonaceous particles emitted by other sectors [START_REF] Lelieveld | The contribution of outdoor air pollution sources to premature mortality on a global scale[END_REF]. 10

Life cycle assessments of livestock production systems

Life cycle assessment (LCA) is a tool based on a precise methodology, used to evaluate all environmental impacts of a functional unit of a product (such as one kg of a certain meat) across all stages of the life cycle of this product. In agriculture, the analyzed cycle is usually limited to "cradle to gate", meaning the use and disposal phases are omitted [START_REF] Dumont | Rôles, impacts et services issus des élevages en Europe[END_REF]. Although in theory particulate matter formation can be assessed among these environmental impacts (again by using atmospheric models and epidemiological studies), much more attention is given to other impact areas in practice. The array of literature reviews on livestock LCA published in recent years mainly summarises impacts on global warming, land use, energy use, eutrophication and acidification, while rarely mentioning air pollution impacts on human health (DeVries et al., 2010;DeVries et al., 2015;[START_REF] Mcauliffe | A thematic review of life cycle assessment (LCA) applied to pig production[END_REF][START_REF] Mcauliffe | Environmental trade-offs of pig production systems under varied operational efficiencies[END_REF][START_REF] Zervas | Life cycle assessment of animal origin products[END_REF]. [START_REF] Mcclelland | Type and number of environmental impact categories used in livestock life cycle assessment: A systematic review[END_REF] confirm this in a systematic literature review, in which they compare the amount of different types of impact categories investigated across 173 livestock LCAs (Figure 4). They show that PM formation is one of the least assessed topics, discussed in only 4% of reviewed studies. Impacts on human toxicity are 8 They estimate cost and benefits for 16 member states who exceeded the 2020 ceiling in 2016. Abatement options include low nitrogen feed, low-emission housing, covered manure storage and urea fertilizer application. 9 For France, the estimated annual benefit was 848.78 M€ with cost ranging between 12.5 and 569.48 M€. 10 Most of the cited studies account for this uncertainty in a sensitivity analysis. from regulatory exemptions and low oversight. [START_REF] Cambra-López | Airborne particulate matter from livestock production systems: A review of an air pollution problem[END_REF] stress for instance that no legislation is in force regarding maximum PM concentrations or emissions neither in agricultural environments nor in and from livestock houses. In the EU, the main legislative instrument to achieve the 2030 objectives of the Clean Air Programme is Directive 2016/2284/EU on the reduction of national emissions of certain atmospheric pollutants. This directive provides a range of suggestions to reduce emissions from agriculture and particularly from livestock manure, but none of them is mandatory. 12 These observations are consistent with the fact that the main pollutants emitted in animal agriculture have been less reduced than in other sectors (see Figure 1). The policy ambition appears to remain relatively low, as is the case of ammonia emissions reduction targets compared to other pollutants in France (Table 1). 13 Table 1: Objectives of emissions reduction compared to emissions in 2005 in France

Pollutant

From 2020 From 2030

Sulfur dioxide (SO2) -55% -77%

Nitrogen oxides (NOx) -50% -69%

Organic compounds (NMVOC) -43% -52% Ammonia (NH3) -4% -13%

Fine particles (PM2.5) -27% -57%

Source: MEEM (2017)

The lack of stringent regulation of agricultural air pollution may be related to different factors. The agricultural sector is well known to be powerful politically, and thus difficult to regulate in general [START_REF] Bonnet | Viewpoint: How to regulate meat consumption to improve health, the environment and animal welfare? Food Policy[END_REF]. 14 Given that world agricultural markets are very competitive, environmental regulation may strongly disadvantage domestic farmers, complicating the implementation of stringent regulation. Since farmers are often poor or very poor, regulation thus also raises strong equity issues. Further, regulation may be perceived by the public as a threat to food security and a barrier to the provision of basic needs. Finally, and more specific to air pollution, one additional 12 In the US, farms are exempted from complying with the three main air pollution laws, namely the Clean Air Act, the Comprehensive Environmental Response, Compensation, and Liability Act (CERCLA), and the Emergency Planning and Community Right-to-Know Act (EPCRA). This under-regulation problem is a longstanding one [START_REF] Ruhl | Farms, their environmental harms, and environmental law[END_REF], but it has apparently worsened recently (see for instance here: https://thefern.org/2019/12/a-breathtaking-lack-of-oversight-for-airemissions-from-animal-farms/). 13 France has received a formal notice in 2020 by the European Commission to improve its air quality legislation (see https://www.lemonde.fr/planete/article/2020/05/28/pollution-de-l-air-la-france-de-nouveau-dans-le-collimateur-de-lacommission-europeenne_6040976_3244.html). Furthermore, it appears that the constraints to register large intensive animal farm units have been made easier in France, thus indirectly increasing the risk exposure of the residents (see for instance: https://reporterre.net/Le-droit-de-l-environnement-est-detricote-au-nom-de-la-simplification). 14 To illustrate this difficulty in another environmental domain, the agricultural sector has been and still is largely exempted from climate policy such as carbon tax or emissions trading systems [START_REF] Springmann | Mitigation potential and global health impacts from emissions pricing of food commodities[END_REF].

factor may be the lack of study in academia, and thus the scientific uncertainty regarding the air pollution health damage due to agriculture. A related issue is the difficulty to attribute the precise origin of an air pollutant. As an example, ammonium nitrate, an important compound of particulate matter (Figure 2) is formed from ammonia, coming almost exclusively from the animal agriculture sector, and nitrogen oxides, coming mainly from road traffic and industry.

We emphasize that modern research in economics may help fill (partially) this knowledge gap.

Indeed, novel quasi-experimental approaches, such as difference-in-difference or regression discontinuity may help disentangle the origin of pollution and provide causal estimates [START_REF] Dominici | Particulate matter matters[END_REF][START_REF] Chabé-Ferret | Évaluation des Politiques Publiques : Expérimentation randomisée et méthodes quasi-expérimentales[END_REF]. 15 However, the only economic study focusing on causal impacts that we are aware of is [START_REF] Sneeringer | Does animal feeding operation pollution hurt public health? A national longitudinal study of health externalities identified by geographic shifts in livestock production[END_REF]. She uses geographic shifts in the animal farming industry to measure the impact of pollution on infant health, finding that a doubling of animal production leads to a 7.4% increase in infant mortality. Hence, future research and the development and accessibility of data should be able to come up with exogenous variation in air pollution generated by agriculture in order to credibly causally assess the air pollution related health impacts. These results should in turn help better design environmental and agricultural policies.

We therefore conclude this discussion by emphasizing the need to develop more studies on the air pollution health impact due to agriculture, especially economic ones. We believe however that limited academic research cannot justify the current under regulation. There are two main reasons for this. First, the results that we presented in the literature survey above emphasize the existence and often the possible severity of health risks to workers, residents as well as to the general population. Current scientific knowledge thus can hardly justify that the risks are nonexistent or small enough to be ignored. Second, if anything, scientific uncertainty due to sparse research and methodological difficulties should justify policy action, and not policy inaction, consistent with the precautionary principle. Hence, the best strategy should be to err on the side of safety while stimulating more research to improve knowledge. We further emphasize that introducing regulation can boost monitoring and measures of agricultural air pollution, and in turn the production of data, thus facilitating academic research. To break the current academic and regulatory inertia, it is thus 15 Let us give a few examples. Using the congestion pricing program in Stockholm, [START_REF] Simeonova | Congestion pricing, air pollution and children's health[END_REF] show for instance that the significant reduction of nitrogen dioxide (NO2) and particulate matter (PM10) relative to pre congestion pricing levels lead to significant reductions in the incidence of childhood asthma in Stockholm in the months and years after the program was implemented. In that case, comparing the before and after difference in pollution for the group receiving the pricing program and its counterpart makes it possible to attribute the pollution reduction to the program. In the same vein, [START_REF] Schlenker | Airports, air pollution, and contemporaneous health[END_REF] use the residual exogeneous daily taxi time variation in daily airport congestion to estimate the population dose response of respiratory and heart-related outcomes to carbon monoxide (CO) exposure. Using China's Huai River Policy, which provides free or heavily subsidized coal for indoor heating during the winter to cities located in the north of the Huai River but not to those in the south, [START_REF] Ebenstein | New evidence on the impact of sustained exposure to air pollution on life expectancy from China's Huai River Policy[END_REF] find a 10μg/m3 increase in PM10 reduces life expectancy by 0.64 years. The findings are derived from a regression discontinuity design based on distance from the Huai River making it possible once again to attribute the pollution reduction to the policy.

urgent that policy makers begin to integrate air pollution in the regulation of animal agriculture. The main regulatory instrument in the EU is the Common Agricultural Policy (CAP). We thus close the discussion with a clear policy recommendation: Regulators must urgently integrate air pollution in the CAP, for instance as a part of the conditionality requirements, as it has already been advanced by several NGOs and policy makers.

  Figure 2 pr nts. We disti s, local resid ion. In Franc

  symptoms, hypersensitivity pneumonitis, and interstitial lung disease

The authors acknowledge Henrik Andersson, Sylvain Chabé-Ferret, Augustin Colette, Anaïs Durand, Johanna Lepeule, Stephan Marette, Bénédicte Meurisse, Nicholas Muller, Arnaud Reynaud and Lea Stapper for useful comments or discussions. Nicolas Treich acknowledges financial support from the ANR under grant ANR-17-EURE-0010 (the Investissements d'Avenir program), INRAE and the FDIR chair.

[START_REF] Casey | Industrial food animal production and community health[END_REF] review 33 publications, 7 of which consider respiratory health impacts.

of these are also among the 7 residential studies on respiratory health reviewed by[START_REF] O'connor | Updated systematic review: associations between proximity to animal feeding operations and health of individuals in nearby communities[END_REF].[START_REF] Douglas | A systematic review of the public health risks of bioaerosols from intensive farming[END_REF] review 17 publications on health effects for residents, covering most of the studies considered by the aforementioned authors.5 The previously mentioned reviews apply a vote-counting approach, which has been shown to be problematic as the probability of finding a real effect decreases with the number of studies included[START_REF] Hedges | Statistical Methods for Meta-Analysis[END_REF]. A preferred method would consist of a meta-analysis corrected for publication bias[START_REF] Andrews | Identification of and correction for publication bias[END_REF].

Supplementary Appendix: Search Strategy for Publications related to Air Pollution and Agriculture

Searches were executed on Web of Science for strings included in publication title, abstract or author-specific keywords. For Fig. 1a, the search string initially generated 287 results. The abstracts of these results were then manually reviewed to attribute each publication to a sector and excluded if they did not focus on pollution.

Search String Fig. 1a TOPIC: ("Air Pollution") AND PUBLICATION NAME: ("Econometrica" OR "American Economic Review" OR "Journal of Political Economy" OR "The Quarterly Journal of Economics" OR "The Review of Economic Studies" OR "American Economic Journal: Economic Policy" OR "The Economic Journal" OR "The Review of Economics and Statistics" OR "Journal of the European Economic Association" OR "Journal of Public Economics" OR "Journal of Environmental Economics and Management" OR "Ecological Economics" OR "Journal of the Association of Environmental and Resource Economists" OR "Review of Environmental Economics and Policy") For Fig. 1b, the number of results for the search strings "Climate Change" OR "Global Warming", "Biodiversity", "Water Pollution", "Air Pollution" were compared (see below for an example). Results for the string "Air Pollution" were manually reviewed to identify publications related to livestock production.

Search String Fig. 1b for "Climate Change" OR "Global Warming" TOPIC: ("Climate Change" OR "Global Warming") AND PUBLICATION NAME: ("American Journal of Agricultural Economics" OR "European Review of Agricultural Economics" OR "Journal of Agricultural Economics") For Fig. 1c, the search string including only "Air Pollution" generated 1043 results. Within these results a refined search was executed to identify publications related to agriculture or livestock production.

Search String Fig. 1c TOPIC: ("Air Pollution") AND PUBLICATION NAME: ("Science" OR "Science Advances" OR "Nature" OR "Nature Communications" OR "Plos ONE" OR "American Journal of Epidemiology" OR "International Journal of Epidemiology" OR "Proceedings of the National Academy of Sciences of the United States of America") Refined by: TOPIC: ("Agriculture" OR "Farm" OR "Livestock")