Mathieu Giraud
email: mathieu@algomus.fr

Using Dezrann in Musicology and MIR Research

Keywords: Music annotation, Corpus encoding, Corpus sharing, User interface, Musicology and MIR research

Easily creating, using, and sharing music annotation data are challenges for digital musicology. How can we formalize music annotation elements, encode them, and use them in systematic musicology approaches? I discuss the motivations in using Dezrann, a web platform to share music analysis [10], for musicology and MIR research.

TALKING ABOUT MUSIC

Informed music analysis should take into many elements beyond the score, with a historical, aesthetic, and comparative perspective. These elements are very challenging for the current Music Information Retrieval (MIR) techniques. Music analysis is nevertheless often built on music annotation, considering the score as Nattiez' niveau neutre [START_REF] Nattiez | Fondements d'une sémiologie de la musique[END_REF]. Music annotation can be seen as putting labels on the score at various positions, possibly with a duration, and sometimes drawing relations between these labels.

The research in the Algomus lab (http://www.algomus.fr) is on computational analysis of tonal music. We focus on the high-level structure of music, such as fugues [START_REF] Mathieu | Computational Fugue Analysis[END_REF] or sonata forms [START_REF] Bigo | Sketching Sonata Form Structure in Selected Classical String Quartets[END_REF], and study analytical elements such as patterns, rhythms, and harmonies. As many researchers in musicology and in MIR, we often engage with scores through reading, annotating, and analyzing. We also discuss other people's analyses. When designing MIR algorithms looking for analytical elements ("Find the cadences"), we frequently go back to scores to compare the results of the algorithms against reference annotations.

More generally, many people analyze, annotate, or talk about music -music teachers and students, researchers, performers. They may need to browse several scores, jumping across sections and commenting music. Music analysis is even not reserved for professionals: Everyone likes to talk about music, to describe what they hear, as demonstrated by presence of music on generic social networks or on platforms such as soundcloud. Can software help people from different backgrounds annotate or discuss music? There are already software for rendering (and sometimes annotating) scores on the web [1,2,[START_REF] Fober | Partitions musicales augmentées[END_REF][START_REF] Hoos | The GUIDO Music Notation Format[END_REF][START_REF] Pugin | Verovio: A library for Engraving MEI Music Notation into SVG[END_REF][START_REF] Stuart | Verovio Humdrum Viewer[END_REF] or for analyzing music [START_REF] Couprie | iAnalyse : un logiciel d'aide à l'analyse musicale[END_REF][START_REF] Huron | Music Information Processing Using the Humdrum Toolkit: Concepts, Examples, and Lessons[END_REF][START_REF] Lepetit-Aimon | INScore Expressions to Compose Symbolic Scores[END_REF]. Some of them have quite remarkable features, but many lack easy usage for everyday use or for people with limited programming background. Our goal here is thus to propose software components for helping people talk about music.

DEVELOPING AND USING DEZRANN

Dezrann (http://www.dezrann.net) is a web platform to share music analysis, mostly developed by Emmanuel Leguy, with inputs and tests from other members the Algomus lab (L. Bigo, L. Feisthauer, R. Groult, F. Levé, MG) as well as student projects. Music is presented either on continuous staves or on waveforms. The user may create, edit and move labels on the score -either on a staff or on spaces above or under the score (Figure 1).

Design goals were to be open-source, multi-platform (no installation, web platform, both desktop and mobile) and simple -kids can use them. Compared to the existing software, we favor the simplicity over the number of features. The usage of the Web components model (Polymer.js) enabled here a clean and modularized design.

Our article at TENOR 2018 [START_REF] Leguy | Dezrann, a web framework to share music analysis[END_REF] details the architecture and the implementation of Dezrann. I will discuss three points at DLfM, focusing on how such a tool may support research in musicology or in MIR, unlocking challenges on creating, using and sharing annotation data:

• Corpus Annotation. People annotating music on paper with a pen make very simple gestures. Dezrann was developed in this spirit, looking for the most natural way to add labels on the score when annotating music. A set of 30 movements in Mozart string quartets was annotated within Dezrann, focusing on sonata form structure and harmony, mainly by Laurent Feisthauer.

• MIR or Musicology Research. We internally use Python scripts based on music21 [START_REF] Scott | music21: A Toolkit for Computer-Aided Musicology and Symbolic Music Data[END_REF], but the underlying annotation format is simple enough to be generated by other languages. During development of MIR algorithms, we frequently load into the platform the annotation results and discuss them -even remotely. For systematic musicology approaches by noncomputer scientists, a perspective towards semi-automated analysis is to add the possibility to call music analysis software within Dezrann.

• Corpus Sharing. As analyses -either manual or automated -can be saved, the platform may be used to share either annotated corpora or the output of annotation algorithms to other people.

Even if the platform is still under development, it is already available and can be used from the public server. Some of its components can also be embedded in another web applications. Beyond usage in Musicology or MIR research, we believe that such platforms enabling to browse and tag music content may have applications for the public. We will soon test Dezrann in secondary schools with children that mostly do not read music.

  Digital Libraries for Musicology (DLfM 2018), 2018, Paris, France Unlocking Musicology Challenge -This work is licensed under the Creative Commons Attribution 4.0 International (CC-BY 4.0) License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA..

Figure 1 :

 1 Figure 1: Annotation with Dezrann on the fugue in C minor BWV 847 by J.-S. Bach, showing both labels on the staves (subjects, counter-subjects, patterns) and on the score (harmonic sequence, degrees, cadence, pedal). Permanent links such as http://dezrann.net/#/bwv847#27 allow jumping at a given place of the score. Annotations are synchronized to the musical time. This synchronization is kept across several views: Here a score generated by Lilypond (top), the autograph D-B Mus.ms. Bach P 202 (middle), and a waveform generated from a performance by Kimiko Ishizaka (bottom, http://welltemperedclavier.org).


Dr Adrien Quintard

Adrien Quintard Young

archive ouverte pluridisciplinaire HAL, est destinée au dépôt

Dr. A. Quintard Optimal organic synthesis is crucial for our future through providing society with access to optimized drugs or materials and also for finding solutions to the energy crisis, for example. In an attempt to find solutions to this challenge of eco-compatible synthesis, our research is centered on the development of new catalytic tools for rapidly transforming easily available molecules into elaborate molecular architectures. To obtain the best possible synthetic sequence while limiting both steps and waste generation, we do not restrict our research to the use of a single catalytic activation mode. Instead, we try to find the best possible catalysts able to selectively activate the desired chemical functions while achieving maximum efficiency. As a result, by combining several compatible catalysts, we can design innovative cascades and rapidly access complex scaffolds. With these im proved synthetic tools in hand, we are able to rapidly elaborate original frameworks with designed functions that should allow us to extend our research in the future to other fields, such as the study of biological interactions, for example.

SYNFORM When did you get interested in synthesis?

Dr. A. Quintard To be honest, after finishing high school, I was more interested in biology than in chemistry. However, despite my poor level in chemistry, I was lucky enough to enter an excellent technical university institute (IUT of Castres in the south of France) where the teachers shared their passion for organic chemistry. At 18 years old, discovering organic synthesis is exciting because you understand that you can build up new materials and drugs while playing, to try to create them with the best possible efficiency. This passion was confirmed during my subsequent internships where I had

Young

SYNFORM What do you think about the modern role and prospects of organic synthesis?

Dr. A. Quintard The current trend is to considerably diminish the role of organic synthesis. Science is full of politics and people who consider that synthesis is a mature field, only here to assist biologists or physicists; they think that it is out of date and that it will soon be replaced by synthetic biology. Synthetic biology is an excellent tool but there are no Swiss Army knives able to perform all required transformations. For now, even though organic synthesis has outstanding potential and can create a wide range of complex structures relatively rapidly, we are still far from what can be considered as ideal synthesis, notably when looking at waste generation and time consumption. Many chemical architectures cannot be constructed on scale and cheaply because of the limitations of our current toolbox. For this purpose, we need to im prove our understanding of the interactions between molecules and the way to activate and transform chemical functions, no tably using catalysts with enhanced efficiency. These improvements at the fundamental level will lead to optimized eco-com patible syntheses and subsequently help physicists and biologists in their own research.

SYNFORM Your research group is active in the area of stereoselective synthesis and natural products synthesis. Could you tell us more about your research and its aims?

A127 Scheme 1 Research overview

Dr. A. Quintard Our research is dedicated to the dis covery of innovative methodologies and eco-compatible synthetic routes to access complex molecular architectures rapidly from simple building blocks. Towards this goal, we focus on the use of catalysis in a wide sense, notably using organocatal ysis (aminocatalysis, chiral bases, thioureas) and metal catalysis based on abundant non-precious metals such as Fe or Cu. In order to go beyond the limitations inherent to the use of a single activation mode, a major part of our efforts focus on the development of multi-catalyzed processes allowing us to achieve the best possible reaction efficiency. Combined with cascade reactions, such an approach at the frontier between metal and organocatalysis allows the solving of key synthetic chemical problems and to considerably shorten some existing routes to crucial synthetic building blocks (Scheme 1).

SYNFORM What is your most important scientific achieve ment to date and why?

Dr. A. Quintard Within our research program, the devel opments on enantioselective borrowing hydrogen is to our eyes among the most interesting. In this work, we have shown that a combination between an iron borrowing hydrogen catalyst and a chiral organocatalyst can promote the enantioselective functionalization of allylic alcohols (Angew. Chem. Int. Ed. 2013, 52, 12883-12887;ACS Catal. 2016, 6, 5236-5244).

Apart from the considerable synthetic economies such a process brings to complex molecule preparation, it was also the first use of an iron catalyst in a borrowing hydrogen. Besides this hydrogen transfer processes, we have also shown more recently that we could rapidly construct valuable keto diols, key direct precursors of 1,3,5polyols, a class of molecules of great interest. By a selective combination between an organocatalyst and a copper complex, a single cascade allowed the enantioselective fluorination of aldehydes and their in situ derivatization by the copper catalyzed aldolization, constructing up to four acyclic stereogenic centers in one single cascade (ACS Catal. 2017, 7, 5513-5517).

  . SYNFORM regularly meets young up-and-coming researchers who are performing exceptionally well in the arena of organic chemistry and related fields of research, in order to introduce them to the readership. This Young Career Focus presents Dr. Adrien Quintard (Aix-Marseille University, France). Biographical Sketch Adrien Quintard studied chemistry at the University of Toulouse (IUT, France) and then at CPE-University of Lyon (France). He completed his PhD in 2011 under the supervision of Prof. A. Alexakis at the University of Geneva (Switzerland). He subsequently moved to the University of Stanford (USA) for a first postdoctoral stay working with Prof. B. M. Trost. After moving to Aix-Marseille University (France) in 2012 for a second post-doctoral position, he obtained in 2013 an ANR research grant allowing him to start his own research program on the development of multi-catalyzed transformations before becoming in 2014 CNRS researcher. Recently, he was recognized by the young researcher 'Emergence' award from the organic division of the French Chemical Society and he was a recipient of a Thieme Chemistry Journals Award. Dr. A. Quintard the opportunity to discover the difficult but stimulating world of research, where the term 'playing' is perfectly appropriate to describe this three-dimensional Lego game chemists are playing.

© Georg Thieme Verlag Stuttgart• New York -Synform 2018/08, A126-A128 • Published online: July 18, 2018 • DOI: 10.1055/s-0037-1609781

L'


Laiyu Lu

Mathieu Chekroun

Odile Abraham

Géraldine Valerie Maupin

Valérie Maupin

Géraldine Villain

Valérie Maupin

Mechanical properties estimation of functionally graded

Keywords: Surface waves, FGM, Inversion, Mortar PACS: 43.20.Hq, 43.35.Cg, 43.35.Zc

INTRODUCTION

The interest in non-homogeneous material systems with gradually varying properties, often referred to as Functionally Graded Materials (FGMs), has recently received considerable attention in several fields of research [START_REF] Chen | On calculating dispersion curves of waves in functionally graded elastic plate[END_REF][START_REF] Baganas | Wave propagation and profile reconstruction in inhomogeneous elastic media[END_REF][START_REF] Carpinteri | Cracks and re-entrant corners in functionally graded materials[END_REF][START_REF] Erdogan | Fracture mechanics of functionally graded materials[END_REF]. As opposed to the abrupt changes encountered in media with piecewise homogeneous layers, the gradual variation in FGM material properties is known to improve failure performance while preserving the intended thermal, tribological and/or structural benefits from combining dissimilar materials [START_REF] Carpinteri | Cracks and re-entrant corners in functionally graded materials[END_REF]. Since certain performance requirements cannot be practically met with spatially uniform or multi-layered material compositions, FGMs have been enjoying widespread use in the fields of aeronautics, astronautics, etc [START_REF] Carpinteri | Cracks and re-entrant corners in functionally graded materials[END_REF][START_REF] Erdogan | Fracture mechanics of functionally graded materials[END_REF].

Another material that can be considered as an FGM is concrete. Concrete is characterized by its extremely heterogeneous nature since, as a common construction material, its composition includes cement or asphalt along with other materials, such as aggregates and water. The potential for a non-destructive evaluation of concrete is of major importance in monitoring the durability of civil engineering structures, especially as regards cover concrete, which is directly exposed to aggressive attack from external sources. Ultrasonic waves are often introduced to characterize such material properties; experiments have shown that wave dispersion and damping are influenced by both grain size variation and the water-to-cement ratio [START_REF] Philippidis | Experimental study of waves dispersion and attenuation in concrete[END_REF][START_REF] Jacobs | Affect of Aggregate Size on Attenuation of Rayleigh Surface Waves in Cement-Based Materials[END_REF][START_REF] Owino | Attenuation Measurements In Cement-Based Materials Using Laser Ultrasonics[END_REF]. It is difficult however to generate a quantitative estimation for these properties based on measured ultrasonic signals, especially for attenuation. The main reason behind this constraint is the multiple-scattering of waves within such a heterogeneous medium, as a result of the random distribution of pores, air bubbles and aggregates when the wavelength has the same order of magnitude as the dimension of heterogeneities. The general conventional method, i.e. homogenization, has been widely employed to estimate equivalent or effective material properties [START_REF] Aggelis | An iterative effective medium approximation (IEMA) for wave dispersion and attenuation predictions in particulate composites, suspensions and emulsions[END_REF]. Even though such estimations provide a reasonable overall prediction of mechanical behavior, they are still insufficient to accurately predict local behavior, as experiments on damaged concrete have also revealed [START_REF] Aggelis | Experimental study of surface wave propagation in strongly heterogeneous media[END_REF][START_REF] Dave | Graded viscoelastic approach for modeling asphalt concrete pavement[END_REF]. The modulus of elasticity for pavement, composed mainly of bituminous concrete, is not constant and in fact varies with depth due to a number of factors including aging, moisture content and temperature [START_REF] Dave | Graded viscoelastic approach for modeling asphalt concrete pavement[END_REF][START_REF] Ullidtz | Pavement analysis, developments in civil engineering[END_REF].

Given the growing interest in FGMs, the need for non-destructive techniques to measure their mechanical properties has become more acute. Shear-wave (S-wave) velocity and attenuation are usually considered key parameters for characterizing the mechanical properties of materials. As a non-invasive method, surface wave analysis has proven efficient in evaluating these parameters; it has been extensively used in many fields at various scales, ranging from earth imaging in geophysics [START_REF] Woodhouse | Mapping the upper-mantle 3-dimensional modeling of Earth structure by inversion of seismic waveforms[END_REF][START_REF] Curtis | Probing the earth's interior with seismic tomography[END_REF][START_REF] Leveque | Waveform inversion of surface wave data:test of a new tool for systematic investigation of upper mantle structures[END_REF][START_REF] Maupin | Love-Rayleigh wave incompatibility and possible deep upper mantle anisotropy in the Iberian Peninsula[END_REF][START_REF] Dziewondki | Parametrically simple Earth models consistent with geophysical data[END_REF] to the exploration of pavement structures in geotechnical engineering [START_REF] Ryden | Guided wave propagation in three-layer pavement structures[END_REF], as well as the estimation of material elastic properties in ultrasonic NDT (Non-Destructive Testing) [START_REF] Wu | Inverse determinations of thickness and elastic properties of a bonding layer using laser-generated surface waves[END_REF][START_REF] Cho | Non-destructive testing of high strength concrete using spectral analysis of surface waves[END_REF][START_REF] Cho | Spectral analysis of surface wave response of multi-layer thin cement mortar slab structure with finite thickness[END_REF]. At most scales, and particularly at the geotechnical scale [START_REF] Park | Multichannel analysis of surface waves[END_REF][START_REF] Alam | Spectrum analysis of seismic surface waves and its applications in seismic landmine detection[END_REF], an analysis of velocity or attenuation constitutes the main tool for extracting useful information from data. Material parameters are generally obtained by means of an inversion technique, which yields an optimal model minimizing the difference between the predicted and measured dispersion (and/or attenuation) curves of surface waves.

In the geotechnical field, this method has been successfully applied in order to infer the properties of homogeneous or multi-layered media [START_REF] Park | Multichannel analysis of surface waves[END_REF], although applications with continuous variations over the depth are still lacking. On the other hand, the inversion of surface wave dispersion for a continuous profile of the Earth's mantle parameters has been frequent in seismology for many years [START_REF] Leveque | Waveform inversion of surface wave data:test of a new tool for systematic investigation of upper mantle structures[END_REF][START_REF] Maupin | Love-Rayleigh wave incompatibility and possible deep upper mantle anisotropy in the Iberian Peninsula[END_REF][START_REF] Dziewondki | Parametrically simple Earth models consistent with geophysical data[END_REF]. Due to the typically small variations of parameters in the mantle (i.e. just a few percent), dispersion curves are inverted using linear methods. In this paper, we will demonstrate that such a method may be extended to shear-wave velocity and attenuation inversions in FGMs featuring much greater variations and contrasts. This demonstration will be conducted first with a synthetic experiment, then by application to non-destructive concrete testing.

THEORETICAL FORMULATION AND METHODOLOGY

Forward Problem

To calculate the predicted dispersion and attenuation for a forward model, one important step involves solving the eigenvalue problem of Rayleigh waves for the visco-elastic model. In linear viscoelasticity, this solution can be obtained by employing the same governing equations as those found in the corresponding elastic problem with identical boundary conditions, by simply replacing real variables with corresponding complex variables. In this paper, the dual conditions of linear viscoelasticity and weak dissipation have been assumed; under such an assumption, the eigenvalue for a visco-elastic model can be derived from the eigenvalue of the pure elastic case through perturbation theory [START_REF] Aki | Quantitative Seismology[END_REF]. As a consequence, the software developed for elastic media can easily be modified to the visco-elastic case.

The eigenvalue problem for Rayleigh waves propagating in purely elastic, layered structures can be described by differential equations in matrix form [START_REF] Aki | Quantitative Seismology[END_REF], i.e.:

df(z) dz = Af(z) ( 1 )
where z is the depth coordinate, and A a 4 × 4 matrix associated with the propagator matrix, which is a function of the depth-varying compressional wave velocity V p , shear-wave velocity V s and density ρ . f(z) is a column vector that comprises the displacement vector r 1 , r 2 and two elements r 3 , r 4 of the stress tensor. Combined with boundary conditions, the precondition for existence of a solution to Eq.( 1) yields the eigenequation of Rayleigh waves, which offers the following general form:

f (c 0 , ω; V p , V s , ρ) = 0 (2)
where c 0 is the phase velocity of Rayleigh waves, and ω the angular frequency. The subscript 0 in this context denotes the elastic case. Both analytical and numerical methods have been developed to treat seismic wave propagation in FGMs [START_REF] Chen | On calculating dispersion curves of waves in functionally graded elastic plate[END_REF][START_REF] Baganas | Wave propagation and profile reconstruction in inhomogeneous elastic media[END_REF]. A general treatment is based on the transfer matrix approach, according to which the medium with continuouslyvarying inhomogeneity is regarded as a stack of many thin, piecewise homogeneous layers [START_REF] Aki | Quantitative Seismology[END_REF].

An explicit formulation of the transfer matrix for layered media can then be obtained.

Another approach to calculating the transfer matrix of FGMs is based on an exact solution in the form of the Peano series of multiple integrals [START_REF] Baron | Impact of localized inhomogeneity on the surface-wave velocity and bulk-wave reflection in solids[END_REF][START_REF] Shuvalov | General formalism for plane guided waves in transversely inhomogeneous anisotropic plates[END_REF]. In this manner, continuous inhomogeneity serves to replace the exponential solution to the wave equation by the Peano integral expression. We adopted a Runge-Kutta scheme herein to numerically integrate the system of differential equations, derive the propagator matrices, and calculate dispersion and attenuation [START_REF] Aki | Quantitative Seismology[END_REF][START_REF] Takeuchi | Seismic surface waves. In seismology: surface waves and earth oscillations[END_REF].

After solving for eigenvalues, specifically the phase velocity for each frequency, the partial derivatives of phase velocity with respect to model parameters can then be obtained from the eigenfunctions by employing the variational principle [START_REF] Aki | Quantitative Seismology[END_REF]. The derivatives with respect to the three parameters ρ, V p and V s are given by the following expressions:

∂c 0 ∂ρ = 1 2ρ ∂c 0 ∂V s V s + ∂c 0 ∂V p V p - 1 2k 2 UI ω 2 (r 2 1 + r 2 2 ) ( 3
)
∂c 0 ∂V p = ρV p 2k 2 UI kr 1 + dr 2 dz 2 ( 4
)
∂c 0 ∂V s = ρV s 2k 2 UI kr 2 - dr 1 dz 2 -4kr 1 dr 2 dz ( 5
)
where U and k are the group velocity and wave number, respectively. I is the energy integration of the Rayleigh wave and equals:

I = ∞ 0 ρ r 2 1 + r 2 2 dz (6)
In this paper, the software package developed by Saito has been used to calculate the Rayleigh wave eigenvalues and partial derivatives of phase velocity for the isotropic elastic model [START_REF] Saito | Disper80: a subroutine package for the calculation of seismic normal mode solutions[END_REF].

For a visco-elastic material, the modulus M * and, by extension, both the body wave and surface wave velocities can be represented by complex quantities:

M * (ω) = M 1 (ω) + iM 2 (ω) ( 7 )
The degree of dissipation is often characterized by the quality factor Q, which can be expressed as:

Q = M 1 M 2 (8)
The complex modulus and velocities depend on frequency ω since the relationship between stress and strain depends on time, as a result of visco-elasticity. In addition, the real and imaginary parts of the modulus are not independent. The relationship known as Kramers-Kronig dispersions [START_REF] Carcione | Wave fields in real media: wave propagation in anisotropic, anelastic and porous media[END_REF], which states that visco-elastic materials are inherently dispersive, must be satisfied. In mathematical terms, this implies that M1 and M2 are Hilbert transform pairs, and this relationship constitutes the necessary and sufficient condition for M to satisfy the fundamental principle of causality [START_REF] Carcione | Wave fields in real media: wave propagation in anisotropic, anelastic and porous media[END_REF].

Laboratory experiments show that over a broad bandwidth(10 -2 -10 7 Hz), Q can be considered as independent of frequency at very low strain levels [START_REF] Toksoz | Attenuation of seismic waves in dry and saturated rocks: I. Laboratory measurements[END_REF]. One commonly used form of the dispersion relation that is able to satisfy the Kramers-Kronig relationship with Q remaining nearly constant is the one developed by Liu et al. [START_REF] Liu | Velocity dispersion due to anelasticity: Implications for seismology and mantle composition[END_REF], which can be written as [START_REF] Aki | Quantitative Seismology[END_REF][START_REF] Kjartansson | Constant Q-Wave propagation and attenuation[END_REF]:

V (ω) V (ω ref ) 1 + 1 πQ ln ω ω ref ( 9
)
where ω ref denotes a reference circular frequency and V the real part of the P-wave or Swave velocity. The dispersion relation is only applicable for weakly dissipative media, e.g. Q > 10, in which the dispersion caused by intrinsic dissipation remains small. Frequencydependent modifications to the velocities introduced by attenuation are then mapped onto Rayleigh wave phase velocity variations, yielding a first-order expression using the variational principle:

c = c 0 + 1 π ln ω ω ref ∂c 0 ∂V s V s Q -1 s + ∂c 0 ∂V p V p Q -1 p dz ( 10
)
where c 0 is the phase velocity in the corresponding purely elastic model. It would appear from this expression that phase velocity dependence on the quality factor of P-or S-waves

is: ∂c ∂Q -1 t = 1 π ln ω ω ref ∂c 0 ∂V t V t ( 11
)
where the subscript t = P, S denotes the compressional or shear wave, respectively.

The spatial damping of Rayleigh waves can be characterized by the dissipation factor

Q -1 R .
For a plane Rayleigh wave propagating with an attenuation coefficient of α :

u = u 0 e -αr e -ikr (12) Q -1 R is related to α by Q -1 R = αc πf (13) Q -1
R can be obtained by measuring and processing the Rayleigh wave amplitude at various distances. Q -1 R depends, at the first order in attenuation, on the variations in quality factors with respect to depth via the following relation:

Q -1 R = 1 c 0 ∂c 0 ∂V s V s Q -1 s + ∂c 0 ∂V p V p Q -1 p dz ( 14
)
This expression yields the partial derivatives of Rayleigh wave dissipation with respect to the quality factors:

∂Q -1 R ∂Q -1 t = 1 c 0 ∂c 0 ∂V t V t (15)

Inverse Problem

For our problem, the relation between data and depth-dependent parameters can be summarized as:

Q -1 R (f ), c(f ) = g(V p , V s , ρ, Q -1 s , Q -1 p ) (16)
where In this paper, the generalized nonlinear inversion technique for continuous problems, as developed by Tarantola and Valette [START_REF] Tarantola | Generalized nonlinear inverse problem solved using the least squares criterion[END_REF], has been used to invert both the velocity and attenuation profiles. This method was designed to minimize the square of the differences between predicted and observed data on dispersion and/or attenuation.

Q -1 R (f )
In expressing Eq. ( 16) in the general form:

d = g(p) (17)
where d and p are the data and parameter sets respectively, then the inverted model at iteration k + 1 will be given according to the least square solution for discrete nonlinear inverse problems (as proposed by Tarantola and Valette [START_REF] Tarantola | Generalized nonlinear inverse problem solved using the least squares criterion[END_REF]) by:

p k+1 = p 0 + C p0 • G T k • (C d0 + G k • C p0 • G T k ) -1 • {d 0 -g(p k ) + G k • (p k -p 0 )} ( 18
)
where G is the matrix of partial derivatives with respect to the model parameters, i.e.:

G iα = ∂g i /∂p α (19)
For the problem treated herein, which can be calculated from the formulation described in Section 2.1, G T is the transpose of matrix G, p 0 the a priori model, d 0 the data vector, g(p k ) the data predicted from the model p k , C p0 and C d0 are the a priori covariance matrices of the parameter and data, respectively.

The equivalent to Eq. ( 18) for problems with continuous variations in model parameters is expressed by:

p k+1 (z) = p 0 (z) + dz i i j C p 0 (z, z ) • G i k (z ) • (S -1 ) ij • d j 0 -g j (p k ) + dz • G j k (z ) • [p k (z) -p 0 (z)] (20)
where the matrix S k is given by

S ij k = (C d0 ) ij + dz dz G i k (z ) • C p0 (z, z ) • G j k (z ) (21)
Theoretically speaking, for a property that varies continuously with depth, we need to perform an inversion for the property at an "infinite" number of depth points in the model. In practice however, we are obviously required to sample the functions at a finite number of depth points, yet we are also intent on maintaining inversion as independent of sampling, which means introducing certain a priori information to constrain the inversion process. As undertaken in [START_REF] Leveque | Waveform inversion of surface wave data:test of a new tool for systematic investigation of upper mantle structures[END_REF] [START_REF] Leveque | Waveform inversion of surface wave data:test of a new tool for systematic investigation of upper mantle structures[END_REF], Maupin and Cara (1992) [START_REF] Maupin | Love-Rayleigh wave incompatibility and possible deep upper mantle anisotropy in the Iberian Peninsula[END_REF] and Debayle and Sambridge (2004) [START_REF] Debayle | Inversion of massive surface wave data sets:Model construction and resolution assessment[END_REF], we introduced the Gaussian-shaped function as the a priori covariance function of the a priori model p 0 :

C p0 (z, z ) = σ(z)σ(z ) exp( -(z -z ) 2 2L 2 ) ( 22
)
where z and z are two depth points, L the correlation length and σ the variance at depth z. This set-up acts as a spatial filter to smooth the model, thereby imposing a correlation between points separated by a distance on the order of L, with σ(z) controlling the amplitude of allowable model perturbation at z. Since this approach insures that the inversion result remains independent of discretization with depth, it is no longer necessary to sample the model with equal spacing, but rather the spacing should be greater than the correlation length.

If more than one parameter requires inversion (e.g. in this case, the shear-wave velocity V s and P-wave velocity V p ), albeit with a certain relationship between the two being expected (perhaps through an expected Poisson's ratio), then the covariance matrix can be used to impose a correlation between the variations of the two parameters:

C p0 [V s(z), V p(z )] = σ Vs (z)σ Vp (z ) exp( -(z -z ) 2 2L 2 )C sp ( 23
)
where C sp is the coupling coefficient between parameters V s and V p , which is capable of varying between 1 and -1. The respective variation for each parameter is controlled by its standard deviation.

In the applications that follow, we will only invert for V s , Q -1 s and/or V p . Our software can also simultaneously invert for ρ or Q -1 p , with V s , V p and Q -1 s , although tradeoffs must be introduced whenever too many parameters are involved in the inversion process. Since

Rayleigh waves are less sensitive to ρ and Q -1 p , their variations are assumed to be negligible during an inversion.

SYNTHETIC EXPERIMENTS

The synthetic model shown in Figure 1 will be discussed in this section. We will begin by calculating the velocity and attenuation of the Rayleigh wave propagating in this model.

Velocity and attenuation will then be adopted as the 'measured' data used for model property inversions. Since the model is known, it proves helpful to validate the algorithm and investigate the effects of input parameters on inversion results, such as initial model and correlation length. As is customary in the seismological literature , Q -1 , i.e. the inverse of quality factor Q, will be used to characterize the material attenuation.

The V s and Q -1 s profiles of this model are defined by Eq.( 24) below: where z is the depth. The Poisson's ratio equals 0.22 and V p is obtained by:

⎧ ⎪ ⎨ ⎪ ⎩ φ(z) = φ(d) 1 + 1 2 φ(0) -φ(d) φ(d) tanh[a(1 -2z/d)] tanh a + 1 0 ≤ z ≤ d φ(z) = φ(d) z > d (24)
V p = 2(1 -σ) 1 -2σ V s (25)
Next, let a = 2, d = 5 in Eq.( 24), with the model being assumed homogeneous below 5cm. [START_REF] Baron | Impact of localized inhomogeneity on the surface-wave velocity and bulk-wave reflection in solids[END_REF][24] adopted the function in Eq.( 24) in order to model a transition layer, in which material properties vary continuously without any abrupt change at the edge points. These authors also discussed the forward problem for the elastic case by introducing the Peano series and went on to offer an analytical expression of the dispersion relation. This model has been extended here to the visco-elastic case in order to investigate the inverse problem. We set

V s (0) = 1.8km/s, V s (d) = 2.3km/s. It is assumed that Q -1 s exhibits the same variation as V s and that Q -1 s (0) = 0.02, Q -1 s (d) = 0.06. Q -1
p is also assumed to satisfy the equation in [START_REF] Anderson | Theory of the Earth[END_REF], i.e.:

Q -1 p (z) = 4 3 V 2 s (z) V 2 P (z) Q -1 s (z) (26)
Figures 1a andb show the velocity and Q -1 s profiles of this model, while Figures 1c andd 151 indicate the velocity and Q -1 of the Rayleigh wave. The velocity for the pure elastic case is 152 also given in Figure 1c.

In civil engineering, the Poisson's ratio, instead of the V p profile, of materials can sometimes be approximately estimated as a priori information. Figure 2 displays the inversion results with a known Poisson's ratio (0.22) as the a priori information. The initial model, inverted profiles and true model are all presented in this figure. The initial value of V s consists of the profile 1.1c -0.5λ . In early geotechnical engineering practice, this profile was often used to approximate the V s profile. The two coefficients (1.1 and 0.5) may at times be modified depending on the expected Poisson's ratio of the medium. We performed a large number of calculations for various initial models, and this profile proves to be a better initial model than a constant profile. The initial values of Q -1 s can be chosen depending on the Q -1 of Rayleigh waves. A morderate value for Rayleigh wave Q -1 turns out to be a good initial model for Q -1 s . The initial V p is generated from V s by means of a known Poisson's ratio. Figure 2 reveals a set of good inversion results obtained for Q -1 s , V s and hence V p . In contrast, without any a priori information on the Poisson's ratio and V p , inverting for V s and V p as the independent parameters may be preferred. Figure 3 presents the inversion results for this particular case. Compared to Figure 2, Figure 3 shows that larger differences between the true and inverted velocity profiles are observed within the 4cm -5cm range for V s and Q -1 s . This case features less a priori information and a higher number of inverted parameters, which in turn increases the level of inversion uncertainty. For P-wave velocity V p , the difference between the inverted and true models is considerable. As a consequence of the low sensitivity of Rayleigh waves to V p , the inverted V p profile has in fact changed very little compared to the initial model. Even though the inverted V p is far from the true model, the inverted phase velocity and Q -1 of the Rayleigh wave display good agreement with the data, which implies that the inverted V p is not reliable and moreover that the same dispersion curves and Q -1 can allow for multiple solutions of V p . Numerical results for the other synthetic test also support this finding. In a practical application therefore, we should seek a priori information on the Poisson's ratio or P-wave velocity of the materials using the other method, such as reflection and refraction method. It should be pointed out however that the incorrect V p profile does not significantly bias the inverted V s and Q -1 s profiles and that these profiles can be properly recovered even in the absence of accurate knowledge on V p or Poisson's ratio.

APPLICATION: DETERMINATION OF MORTAR PROPERTIES FROM LASER MEASUREMENTS

In this section, the method described above will be applied to perform an inversion for the shear-wave velocity (correspondingly the Poisson's ratio) and for the attenuation of a mortar sample, based on surface wave data collected using a laser interferometer.

Experimental set-up and measurements

Experimental measurements have been carried out on mortar samples with a maximum grain size of D max. = 4 mm. Two series of mortar slabs M1 and M2, differing in just their water/cement ratio, were considered: the M1 series has a low water/cement ratio (w/c = 0.35), while the M2 ratio is higher (w/c = 0.65), thus inducing higher porosity. CEM 1 52.5N CE CP2 NF cement has been used and the granulate are silico calcareous. The slabs were held underwater in between experiments to ensure remaining fully saturated at all times.

Each series was composed of 5 identical slabs with dimensions 600 mm × 600 mm × 120 mm.

The 120 mm specimen thickness was considered sufficient to ensure that signals received at the surface corresponded to Rayleigh waves, thus avoiding the generation of Lamb wave modes.

A piezoelectric transducer equipped with a wedge was used as a source to generate Rayleigh waves in the mortar slabs. The source function was a Ricker wavelet with a central frequency equal to 120 kHz. Reception was performed with a laser interferometer (Tempo from Bossa Nova Tech), which acquires the normal displacement at the slab surface according to a non-contact protocol. The laser beam position was controlled by a robot to allow for an acquisition every 1 cm at a distance from the source ranging from 10 cm to 45 cm, to yield the equivalent of a common-shot gathers in seismology. The precision of the laser beam position is better than 0.01mm and the data acquisition card has a sampling rate equal to 10MHz and a 16 bits resolution [START_REF] Bretaudeau | Small-scale modeling of onshore seismic experiment: a tool to validate numerical modeling and seismic imaging methods[END_REF].

To take into account the heterogeneous nature of the mortar (i.e. size and position of sand, presence of bubbles and other surface inhomogeneities [START_REF] Chekroun | Caractérisation mécanique des premiers centimètres du béton avec des ondes de surface[END_REF]), a total of 36 similar common-shot gathers were collected at different positions on the 5 slabs for each series; a spatial average could then be established to obtain the coherent field.

Geometrical spreading was corrected from the measurements by multiplying all signals by √ r, where r is the distance from the source. We used a pω transform to extract the phase velocity dispersion curves [START_REF] Mokhtar | Seismic velocity and Q model for the shallow structure of the Arabian shield from short-period Rayleigh waves[END_REF], where p represents the slowness of the waves (p = 1/c) and ω the angular frequency. This method transform the multi-channel data wave field into the slowness-frequency domain. In p-ω domain, the maximum will be reached at the eigenvalues of the Rayleigh wave. The algorithm proposed by Herrmann is used to extract the velocities for each frequency and also provides error bars [START_REF] Herrmann | Computer Program in Seismology[END_REF]. The attenuation factor is estimated from the decrease of the amplitude spectrum of the coherent field during propagation. Damping factor vs. frequency was evaluated by performing a linear fit of the natural logarithm of the spectral amplitude [START_REF] Chekroun | Analysis of coherent wave dispersion and attenuation for non-destructive testing of concrete[END_REF].

Figure 4 provides the experimental phase velocity and error bars. Due to limitation in the transducer frequency band as well as the signal-to-noise ratio, only data in the 60-180 kHz bandwidth could be introduced. This frequency range corresponds to wavelengths ranging from approximately 10mm to 40mm [START_REF] Chekroun | Caractérisation mécanique des premiers centimètres du béton avec des ondes de surface[END_REF]. The spectrum modulus of the Rayleigh waves at this bandwidth is larger than -20dB. The interval between two adjacent frequencies is 2440Hz. Details about the experiment and extraction of the dispersion curves can be found in reference [START_REF] Chekroun | Analysis of coherent wave dispersion and attenuation for non-destructive testing of concrete[END_REF].

Dispersion of Rayleigh waves may arise from two phenomena: the variation with depth of the properties of the media and the dispersion of P and S waves related to multiple scattering in heterogeneous media. The scattering produces in addition attenuation. It is therefore in theory possible to distinguish between the two effects by a combined analysis of dispersion and attenuation. However, the frequency range used here corresponds to wavelength varying from 10 to 40 mm, while in the mortar series M1 and M2 the maximum grain size is D max = 4 mm. Then the wavelength is much bigger than the heterogeneities of mortar, inducing negligible scattering effects. The scattering-related dispersion of P and S waves in the mortars used here and in associated concrete samples has been studied in reference [START_REF] Chekroun | Caractérisation mécanique des premiers centimètres du béton avec des ondes de surface[END_REF].

They find that for the concrete samples, which in essence contain bigger heterogeneities (coarse aggregates of D max = 20 mm) than the mortar samples, scattering related to heterogeneities generates a noticeable additional dispersion compared to the mortars, in the actual frequency range. In addition, the scattering by aggregates produces an attenuation about three times the one observed for mortars. We conclude that for the mortar used here, the major part of the dispersion is likely to originate from depth-dependence of the structure.

Inversion results

0.6 0.8 Figure 5 shows the inversion results for the M1 series. The correlation length has been set at 0.5cm, and the a priori variance is 0.002km/s for V s and V p and 2 × 10 -4 for Q -1 . In comparison with synthetic models, a greater number of iterations was required here before obtaining convergence. We conducted a total of 100 iterations and typically reached convergence after between the 50th and 60th iteration. In order to compare results for differenta priori information, the inversion procedure was carried out in two ways: first with an unknown Poisson's ratio ν, where V s , V p and Q s were inverted simultaneously; and second with a known Poisson's ratio given as a priori information, where only V s and Q s were inverted.

No real significant differences are observed in Figure 5, except for the larger a posteriori error obtained on the inversion with an unknown Poisson's ratio. This finding is not difficult to explain since the constraint imposed on the inversion procedure without an a priori Poisson's ratio was less stringent. Similar trends and orders of magnitude were observed for these series of results. Figure 4 indicates that these two inverted dispersion and Q -1 curves fit the experimental data quite well. Based on this example and in recalling the discussion on numerical models in the previous section, it is considered an appropriate choice to perform the inversion with an a priori Poisson's ratio estimated from experience. Poisson's ratio set equal to 0.2. As expected, it was observed that the shear-wave velocity of mortar M1 exceeded that of mortar M2, since M2 is more porous. Similarly, the quality factor Q s of M1 was greater than that of M2 at all depths. For both mortars, a slight dispersion in the phase velocity dispersion curve is visible, hence indicating that material properties vary with depth. At low frequencies, phase velocity increases such that a higher body-wave shear velocity with depth can be expected. This observation is in accordance with common knowledge held on concrete properties. The first few millimeters contain fewer aggregates than lower depths due to a wall effect: the proportion of large aggregates becomes constant after a depth equal to the radius of the largest aggregate [START_REF] Zheng | Aggregate distribution in concrete with wall effect[END_REF]. To be able to explain this increase of wave velocity, we performed density profile measurements by gammadensitometry [START_REF] Villain | Gammadensimetry: A method to determine drying and carbonation profiles in concrete[END_REF]. In the slabs the density is increasing towards the surface(See Figure 8 for the density variation). We suppose that this phenomenon is due to the use of a wood form that is coated with bakelite so that bleeding water is not absorbed by this formwork on the contrary to wooden formwork classically used for concrete structures. As a consequence, this higher quantity of water available near the surface, during setting, is increasing density by chemical reaction in this area compared to the depth. Furthermore, the grain size is surely modifying the density near the surface, and this phenomenon is competing with pore size and pore distribution near the surface together with water gradient that can extend to different depth [START_REF] Villain | Gammadensimetry: A method to determine drying and carbonation profiles in concrete[END_REF]. For the inversion step, a constant Poisson's ratio of 0.2, i.e. a classical value for mortar, was chosen. It is apparent on the shear-wave velocity profiles that below 1.5 cm, the mortar can be considered homogeneous with a shear-wave velocity equal to 2,660 m.s -1 for M1 and 2,330 m.s -1 for M2. Figure 7 shows the experimental phase velocity and attenuation with error bars, together with the corresponding inverted one.

The depth limits for the two inverted profiles are slightly different because the higher phase velocity in M1 yields slightly larger wavelengths, and therefore larger penetration depths, than in M2, for the same frequency range. The inverted V s profiles for M1 and M2 tend to display a typical characteristic, namely increasing smoothly with depth over the first few millimeters near the surface, with an inflection point observed at 1.4 -1.6cm, then tending to a constant value at greater depths. This finding is a result of aggregate size, the presence of air bubbles and a saturation rate capable of differing nearer the surface and deeper due to the wall effect and exposure to air. Consequently, the density may feature similar characteristics. To verify these inverted profiles, we measured the mortar sample density at various depths using gamma-densitometry. It should be pointed out this does not mean we attribute the velocity variation only to the density variation. The shear wave velocity variation is related to several parameters such as density, Young's modulus and water content, often in competition.

Figure 8 shows the density profiles with errors measured by means of gamma-densitometry for mortars M1 and M2; average density equals 2,256±5 kg/m 3 for M1 and 2,151±17kg/m 3 for M2. A difference in density is noticed near the edges, as a result of the skin effect. For the sake of comparison, we used the function with exponential attenuation to fit measurements (counting from the surface):

f = f 0 + Ce -x/a ( 27
)
where f and x denote density and depth, respectively. C and a are the constants to be determined. The function and its graph are shown in Figure 8. For M1, the profile is nearly symmetrical on both sides. At a depth of 1.5-2.0cm from the surface (as denoted by two circles drawn in a dashed line), a similar inflection point can be observed. The density tends to a constant between the two inflection points. The profile shape and inflection point location closely match the V s profile. For M2, the density profile is not symmetrical from both sides of the surface. This outcome is likely due to the presence of air bubbles near one side (where depth equals 120mm on Fig. 8). Air bubbles are usually removed by hammer blows for different samples on both sides of all slabs. For mortar M2 however, the protocol perhaps differed on the side of the slab used for gamma-densitometry coring samples, which could explain the difference in density near the 120mm deep side. For mortar M2 near the other side (i.e. 0mm depth), the same density variations as for mortar M1 can be observed, which is why density measurements have, to some extent, confirmed our inverted Vs profile and thereby offer a reliable verification for our method.

DISCUSSION AND CONCLUSION

We have presented herein a method for inverting surface wave dispersion and attenuation in terms of the depth-dependent visco-elastic parameters of Functionally Graded Materials (FGMs). The particularity of this method lies in the fact that the visco-elastic parameters are allowed to vary continuously with depth and that both the forward problem (involving calculation of dispersion and attenuation) and the inverse problem (evaluation of viscoelastic parameters) effectively take into account the continuous nature of these variations.

The forward problem solves the equation of elastic motion using a Runge-Kutta integration scheme. The viscous part is treated as a first-order perturbation to the elastic structure, which limits method applicability to materials with weak attenuation.

The inverse problem is treated using the nonlinear solution to continuous inverse problems developed by [START_REF] Tarantola | Generalized nonlinear inverse problem solved using the least squares criterion[END_REF] [START_REF] Tarantola | Generalized nonlinear inverse problem solved using the least squares criterion[END_REF]. This set-up allows for continuous variations, without having to define a set of functions over which the profile is to be decomposed. The inversion step is dependent on a prescribed correlation length for the particular profile and not on its sampling with depth. This protocol is very flexible and enables representing a large set of models very easily. Our software introduces interfaces at prescribed depths, yet this option has not been adopted in the present application. The starting model is defined by a simple formula directly related to the dispersion curve. Problem non-linearity is taken into consideration by iterating linear inversions; we have shown that convergence towards the correct profile is obtained after just a few iterations in the purely elastic case, but a larger number of iterations is required in the case of attenuation.

Dispersion and attenuation have been inverted simultaneously for the S-wave velocity, attenuation and, in some cases, for P-wave velocity. The dispersion due to attenuation has also been incorporated, thus avoiding erroneous mapping as a result of depth-dependent, S-wave velocities. P-wave velocity variations can also be accounted for either by including them as a parameter to be inverted or by coupling their variations to those of the S-wave velocity through a fixed Poisson's ratio value. Partial coupling via a correlation coefficient between 0 and 1 is also possible.

The Rayleigh wave phase velocity and attenuation measured by means of laser ultrasonic experiments are used to infer the depth-dependence of S-wave velocity and attenuation on mortar samples. It has been found that mortar inhomogeneity is confined to the first 1.5cm of depth. This depth-dependence compares well with that of the sample density inferred from gamma-densitometry.

The procedure implemented has been designed to analyze the variations in elastic param-eters with depth; moreover, we assumed herein that the material is homogeneous in the two horizontal directions. This method however may still be combined with other techniques in order to recover the depth-dependence and lateral variations of the elastic parameters. In the case of randomly distributed lateral heterogeneities, such as those in the mortar samples, stacking several recordings has yielded information on the average structure. In the case of lateral variations that are consistent with respect to the length of analyzed surface waves, it can be shown that the dispersion measured between two points depends on the average structure between points [START_REF] Woodhouse | Surface waves in a laterally varying layered structure[END_REF]. Average structures along many paths can then be interpreted within a 3-D structure using various tomographic techniques [12-16, 33, 43]. As an alternative, local measurements of dispersion and their inversion with depth can serve to map lateral variations more directly [START_REF] Lin | Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array[END_REF].

Finally the possibility of using surface wave to investigate variation of properties of the first centimeter of concrete (cover concrete) are underway. A major step will consist in dealing with the strong scattering on our measurements as the wavelength will be of the same order of grandeur as the aggregates (few centimeters).

  and c(f ) are respectively the dissipation factor and phase velocity of the Rayleigh wave at each frequency f . First-order variations in the data compared to variations in model parameters are given in Equations 3, 5, 11 and 15. The inversion step consists of identifying the model or class of models that predicts measured data as closely as possible.

Figure 1 :

 1 Figure 1: a) velocity profiles; b) profile of the inverse of quality factor Q s ; c) phase and group velocity of the Rayleigh wave propagating in this model; and d) inverse of the Rayleigh wave Q-value

Figure 2 :

 2 Figure 2: Inversion results with a known Poisson's ratio. (a),(b),(c) and (d) are the profiles of Poisson's ratio, S-wave velocity, P-wave velocity and Q -1 , respectively.

Figure 3 :

 3 Figure 3: Identical to Fig.2, yet without knowing the value of Poisson's ratio.

Figure 4 :

 4 Figure 4: Rayleigh wave phase velocity (a) and Q -1 (b) of both the experimental and inverted models for mortar M1

Figure 5 :

 5 Figure 5: Inversion results for mortar M1 -two results are shown, one with an unknown Poisson's ratio and the other with a known estimated Poisson's ratio equal to 0.2

Figure 6 :

 6 Figure 6: Inversion results for mortars M1 and M2.

Figure 6

 6 Figure 6 presents the inversion results for M2 and M1, with a known and constant

Figure 7 :

 7 Figure 7: Rayleigh wave phase velocity (a) and Q -1 (b) of both the experimental and inverted models for mortar M2

Figure 8 :

 8 Figure 8: The density profiles for M1 and M2 measured by gamma-densitometry.

ACKNOWLEDGMENTS

The authors would like to thank the Electricit de France (EDF) utility company for financing a portion of the research. This work effort was also partially supported by the European INTEREG IIB program MEDACHS and the French ANR project SENSO, as well as by the National Natural Science Foundation of China (40974032).


Julien Gigault
email: julien.gigault@univ-rennes1.fr

Marianne Balaresque

Hervé Tabuteau
email: herve.tabuteau@univ-rennes1.fr

Estuary-on-a-chip: unexpected results for nanoparticles fate and transport

Keywords: nanoparticles, fate, estuary, aggregation, microfluidic, lab-on-a-chip

The aim of this work is to evaluate the influence of a spatial salinity gradient, as encountered in an estuarine when fresh water meets salty one, on the aggregation features of fullerene nanoparticles aggregates (nC 60 ). To model these environmental conditions, we design a specific microfluidic device. In literature, physical phenomena occurring in estuarine area, such as the salinity gradient and the flow conditions, were never considered. Our results suggest that even a short time exposure (couple of seconds) to a salinity gradient, i.e. dynamical spatial conditions, profoundly affects the aggregation properties of nanoparticles and therefore their environmental behavior and fate. It appears that a salinity gradient could enhanced the stability of nanoparticles, especially at high ionic strength (closed to the seawater level one). Contrary to the common approach of particle aggregation in the lab in environmental conditions, our study shows the need to reconsider the spatial and temporal variation of abiotic parameters (ionic strength, pH, organic matter concentration) in the understanding and the evaluation of nanoparticles environmental fate and transportation mechanisms.

Introduction

Because the main route of contaminant transport is river, estuarine systems are one of the most targeted environmental compartment. Those ones are known to play a decisive role in the conditioning of the colloidal materials life cycle [START_REF] Mayer | Aggregation of Colloids in Estuaries[END_REF][START_REF] Lasareva | Transport of the Colloid Matter of Riverine Runoff through Estuaries[END_REF] . Estuarine systems are characterized by specific physical and chemical conditions: salinity gradient, wastewater income, organic matter, and human pressure (industrial and agricultural activities). From all these parameters the ionic strength, characterized by the salinity, and the organic matter are the main parameters that affect the stability of micro-nanoscale materials [START_REF] Mayer | Aggregation of Colloids in Estuaries[END_REF][START_REF] Hotze | Nanoparticle Aggregation: Challenges to Understanding Transport and Reactivity in the Environment[END_REF] .

Over the last ten years, the influence of the ionic strength (IS) on the physical and chemical properties of nanoparticles is one of the most studied parameter, for both carbon-based (fullerene, carbon nanotubes) [START_REF] Yang | Ionic Strength and Dissolved Organic Matter Alter Aggregation of Fullerene C60 Nanoparticles Suspensions in Wastewater[END_REF][START_REF] Chen | Aggregation and Deposition Kinetics of Fullerene (C60) Nanoparticles[END_REF][START_REF] Chen | Influence of Humic Acid on the Aggregation Kinetics of Fullerene (C60) Nanoparticles in Monovalent and Divalent Electrolyte Solutions[END_REF][START_REF] Brant | Aggregation and Deposition Characteristics of Fullerene Nanoparticles in Aqueous Systems[END_REF][START_REF] Chen | Aggregation Behavior of Single-Walled Carbon Nanotubes in Dilute Aqueous Suspension[END_REF][START_REF] Gigault | Size Characterization of the Associations between Carbon Nanotubes and Humic Acids in Aqueous Media by Asymmetrical Flow Field-Flow Fractionation Combined with Multi-Angle Light Scattering[END_REF] and metal-based (Ag, Au, TiO 2 ) nanomaterials [START_REF] Badawy | Impact of Environmental Conditions (PH, Ionic Strength, and Electrolyte Type) on the Surface Charge and Aggregation of Silver Nanoparticles Suspensions[END_REF][START_REF] Domingos | Aggregation of Titanium Dioxide Nanoparticles: Role of a Fulvic Acid[END_REF][START_REF] Delay | Interactions and Stability of Silver Nanoparticles in the Aqueous Phase: Influence of Natural Organic Matter (NOM) and Ionic Strength[END_REF][START_REF] Rahman | Transport of Aluminum Oxide Nanoparticles in Saturated Sand: Effects of Ionic Strength, Flow Rate, and Nanoparticle Concentration[END_REF][START_REF] Merschel | Contrasting Impact of Organic and Inorganic Nanoparticles and Colloids on the Behavior of Particle-Reactive Elements in Tropical Estuaries: An Experimental Study[END_REF] . To do so, the kinetic of nanoparticles aggregation is generally determined by monitoring their size evolution according to the time after an addition of a known concentration of electrolyte in the bulk. However, since this measurement is made while the bulk salt concentration is supposed to be homogeneous, it cannot be representative to what happens in the natural environment. Indeed, these experiments are performed in static or large scale dynamic reactors (such as columns) that are not representative of the natural conditions in which there are spatial and temporal variations of the salt concentration. Indeed, the major issue in the environmental media studies (in bulk liquid and porous media) is the limited applicability to the real dynamic process affecting the transport of particles at environmental scale. More generally NPs behaviors in the environment results from various physical and chemical heterogeneities that the NPs come across. For estuarine system, the salinity gradient is one of these spatial and temporal heterogeneities, and up to now nobody has considered its influence on the fate of NPs.

Recently, new experimental methods and capacities, based on microfluidic devices, are opening the door to environmental science for better understand the fate of species in dynamic systems, closer to real natural conditions [START_REF] Sun | Mapping the Salinity Gradient in a Microfluidic Device with Schlieren Imaging[END_REF][START_REF] Drescher | Biofilm Streamers Cause Catastrophic Disruption of Flow with Consequences for Environmental and Medical Systems[END_REF][START_REF] Dunphy Guzman | Influence of Surface Potential on Aggregation and Transport of Titania Nanoparticles[END_REF][START_REF] Stocker | Rapid Chemotactic Response Enables Marine Bacteria to Exploit Ephemeral Microscale Nutrient Patches[END_REF] . The aim of this work is to evaluate the fate of nanoparticles, using fullerene aggregates (nC 60 ), in a salinity gradient. To do that we developed a specific microfluidic device to simulate an environmentally representative salinity gradient. Our results, compared to classical batch experiments, demonstrates a different and specific aggregation kinetics, which raise new questions and concerns of the existing knowledge on the environmental fate of nanoparticles in estuarine system. Experimental

Sample preparation

nC 60 were prepared according to classical procedure previously described [START_REF] Brant | Characterizing the Impact of Preparation Method on Fullerene Cluster Structure and Chemistry[END_REF][START_REF] Labille | Hydration and Dispersion of C60 in Aqueous Systems: The Nature of Water-Fullerene Interactions[END_REF][START_REF] Gigault | Selection of an Appropriate Aqueous Nano-Fullerene (NC60) Preparation Protocol for Studying Its Environmental Fate and Behavior[END_REF] . Typically, 80 mg of C 60 powder (VWR, Fontenay-sous-Bois) were placed in an Erlenmeyer flask with 200 mL of deionized (DI) water (18.2 MΩ.cm at 25°C) and continuously stirred for 2 weeks.

Finally, the final fullerene aggregates solution (nC 60 ) was filtered at 0.45 µm to remove large aggregates. The final concentration of nC 60 is 7.2±0.7 µg g -1 and was characterized by both Total Organic Carbon analyzer (Shimadzu, TOC-Vcph) and UV-Vis (Perkin Elmer spectrophotometer, Lambda 650). The corresponding size distribution of the nC 60 initial stock solution is d zH =180.7 nm (with a polydispersity index, PDI=0.07) and is determined by in-situ dynamic light scattering (Vasco-Flex, Cordouan Technologies, France). Preparation of nC 60 solutions were replicated 8 times with similar size distribution (see Electronic Supplemental Information, Fig. S1).

Microfluidic reactor design and operation

We design a microfluidic device to simulate the estuarine salinity gradient (figure 1). We used standard soft lithography to make this microfluidic channel [START_REF] Duffy | Rapid Prototyping of Microfluidic Systems in Poly(Dimethylsiloxane)[END_REF] . First we used a negative photoresist resin (Su8 2050 Microchem) that was spin coated on a silicon wafer, heated and cured with UV light through a mask on which the outline of the channel is printed. The resin is again heated and developed in a bath of PGMEA to remove the uncured part of the resin.

Thereafter we cast polydimethylsiloxane (PDMS, Sylgard 184 silicone elastomer kit, Dow Corning) on the SU8 master. Once cured at 70°C for an hour the PDMS slab becomes solid, it is removed from the mold and punched to get the inlets and outlets holes of our channel. Then the channel is sealed by gluing the slab on a glass plate using a plasma cleaner (Harrick Plasma). This process renders the surface of the PDMS channel hydrophilic. To recover a hydrophobic surface we put the microfluidic device in an oven at 120°C for two hours. We inject fluid in the channel thanks to syringe pumps (KDS) at the inlets. In a typical experiment, we first saturate the channel with pure water and then connect pure water or salty on one inlet and a dilute fullerenes suspension on the other. The flow rates were controlled by syringe pump from 0.1 to 1.0 ml/h. At those rates the channel is filled in less than a couple of seconds.

The microfluidic device is composed of 5 mm long and 200 µm wide inlets and outlets ("A out" and "B out") arms on either side of a 3cm long channel with a width of 800 µm and a height of 75 µm (figure 1). In the following we focus on a given flow rate equals to 1.0 ml/h but we get similar results over a wide range of flow rates (0.1 to 1.0 ml/h). The Figure 1 illustrate the microfluidic device used to simulate the estuarine salinity gradient.

Each microfluidic experiment are performed in three steps:

1-Salt calibration. We first inject pure water and a salty blue dye solution (Patent Blue from Sigma Aldrich, #198218). This configuration allows us to visualize the salinity gradient on the image that can be analyzed (Fig. S2 top). We checked that the dye does not affect the salinity gradient. The average salt concentration at the two outlets were determined thanks to the UV-spectra of the NaCl stock solution for various salt concentration (Fig. S2 middle and bottom). To this end we collect samples at both outlets during one hour and then we determine absorbance spectrum of the salty solutions.

2-Diffusion nC 60 in the MD. We inject, at a flow rate of 1 ml/h, nC 60 dispersed in water in one arm and pure water in the other. We collect samples at the two outlets during one hour, we stop the flow and then we determine the size distribution of the NPs. We found the same distribution for both outlets (Fig. S3). This check step confirms that nC 60 diffuse around as it flows through the microchip, i.e, the particle concentration is symmetric with respect to the middle vertical plane of the straight channel.

3-nC 60 aggregates in the salinity gradient. We inject at a flow rate of 1 ml/h in each inlet arm respectively dilute suspension of nC 60 in DI water and a salty dye solution (patent blue dissolved in 500 mM of NaCl), a salinity that is encountered in the vicinity of Estuary mouths (∼30 g L -1 ). In this experiment the average NaCl concentration is 42 mmol L -1 and 427 mmol L -1 for A and B outlets, respectively. The residence time of the nC 60 in the salinity gradient is lower than 3s in average. Again, after one hour of sample collection at the two outlets the flow rate was stopped and we start measuring the temporal evolution of the nC 60 size distribution at the two outlets for two hours. Thereafter we compare the temporal evolution of the size distribution of the nC 60 aggregates with that of a batch experiment as described in literature, i.e. dropwise concentrated NaCl dispersed into a nC 60 solution and then manually stirred. All the experiments were replicated five times using for each replicate a new nC 60 stock solution and NaCl preparation.

Size characterization by Dynamic Light Scattering

We performed dynamic light scattering measurements within the collection vials connected to both outlets in order to both monitor and determine in real time the presence and the size distribution of fullerene nanoparticles aggregates transported through the microchip.

Thanks to an optical fiber remote head, the system used in this work is a dedicated DLS system for contact less in-situ measurement (see ESI, Fig S4). For the nanoparticles setup, DLS specifications were: laser of 658 nm with a 60mW power; the scattering angle is 170°, the working distance is fixed at 80 mm.

For nC 60 size kinetics determination with the NaCl, the intensity-weighted mean diameter (i.e. called the z-average hydrodynamic diameter, d zH ) variation was used and was derived using the Cumulants algorithm from the auto-correlation function (ACF) [START_REF] Koppel | Analysis of Macromolecular Polydispersity in Intensity Correlation Spectroscopy: The Method of Cumulants[END_REF][START_REF] Hallett | Particle Size Analysis by Dynamic Light Scattering[END_REF] . Cumulants is the simplest analysis of the ACF and is ideal to determine the average particle size of a population. For obtaining size distribution from the autocorrelation curve (ACF function), two algorithms were used: Sparse-Batesian-Learning (SBL) and the Cumulants. Cumulants algorithm is one of the most famous and recognized method for determining size distribution of polydisperse sample (Kopel). In the case of environmental and/or highly polydisperse sample, Cumulants method presents several disadvantages. Indeed, this algorithm doesn't permit independent determination of the long-time baseline of the ACF function and can lead to inconsistent results when different numbers of data points are included in the fit. Using the baseline as a floating parameter makes it possible to detect problems in the data and to fit data when a little bit of noise is present, which our case. The d H size distribution was determined using the Sparse Bayesian learning (SBL) algorithm for the Laplace transform inversion [START_REF] Nyeo | Sparse Bayesian Learning for the Laplace Transform Inversion in Dynamic Light Scattering[END_REF] . This algorithm is adapted to identify, discriminate and determine the presence of several polydisperse population from the ACF coming from the scattered light [START_REF] Nyeo | Sparse Bayesian Learning for the Laplace Transform Inversion in Dynamic Light Scattering[END_REF] . More specifically, it generates a set of sparse solutions of several domains but of the same dimension. Then, it superimposes these solutions to give a global solution with its dimension treated as a regularization parameter. An optimal solution providing a reliable description of a colloidal suspension is determined by an L-curve for selecting the suitable value of the regularization parameter. So, the final result corresponds to the most probable physical solution of the algorithm which position is easily defined as the inflexion of a L-curve representation. This method was proved to be more robust than classical methods used to determine size distribution of polydisperse sample.

Results and Discussion

Here our objective is to show that the aggregate growth of nanoparticles depends strongly on the homogeneity of the salt concentration. To do this we compare on Figure 2 the temporal variation of the average hydrodynamic diameter, d zH, of nC 60 aggregates either obtained by dropwise NaCl into the bulk, at two concentrations 85 and 500 mmol L -1 , (named bulk experiments, BE), or collected at the outlets of our microfluidic device (named MD), A out and B out , with an average salt concentration of 42 and 427 mmol L -1 , respectively. Both salt concentration in the BE have been chosen to be very close to those obtained at the two outlets of the MD. We choose the two salt concentrations in order to be on both sides of the critical coagulation concentration (CCC) of nC 60 , which is between 120 and 260 mM of sodium chloride [START_REF] Yang | Ionic Strength and Dissolved Organic Matter Alter Aggregation of Fullerene C60 Nanoparticles Suspensions in Wastewater[END_REF][START_REF] Chen | Aggregation and Deposition Kinetics of Fullerene (C60) Nanoparticles[END_REF][START_REF] Chen | Influence of Humic Acid on the Aggregation Kinetics of Fullerene (C60) Nanoparticles in Monovalent and Divalent Electrolyte Solutions[END_REF][START_REF] Meng | Aggregation Rate and Fractal Dimension of Fullerene Nanoparticles via Simultaneous Multiangle Static and Dynamic Light Scattering Measurement[END_REF] according to the different formation pathways [START_REF] Gigault | Selection of an Appropriate Aqueous Nano-Fullerene (NC60) Preparation Protocol for Studying Its Environmental Fate and Behavior[END_REF] .

BE is generally performed to mimic the evolution of the aggregation features of nanoparticles in salty natural environments, i.e., in the ocean or in Estuaries. In this case, the aggregation of nC 60 takes place in a spatially homogeneous salt concentration. Indeed the salt dispersed at the beginning of the experiment diffuse so quickly throughout the vial that no significant aggregates growth of nC 60 is observed during this homogenization process. In the MD experiment the fullerene aggregates (nC 60 ) quickly passed through a NaCl gradient and then they can grow into a vial at rest, in the same conditions than those of the BE.

-Below the CCC of (nC 60 ) there is a slight increase of the aggregates sizes at the beginning of the experiment (t< 15 min) for the BE and then there is no temporal evolution of d zH anymore (Fig. 2a, blue circle). The screening of the surface charge is limited at electrolyte concentrations below the CCC and thus nC 60 aggregation is limited. In this electrolyte conditions, the nC 60 suspension remains quite stable.

We found similar results for the microfluidic experiment in the outlet A out (Fig. 2a, black circle). This is further confirmed by looking at the global size distribution summed over the total kinetics time for both types of experiment (Fig. 2b), which are very similar. The mean aggregate size is a bit larger in the BE than in the MD experiment since the salt concentration is slightly higher in the former case. However, it is worth noting that there is a small micron size population (d zH >1µm) in the MD experiment, which is not the case in the BE. We will get back to this point later on.

-Above CCC of (nC 60 ) the situation is completely different whether the nC 60 aggregates encounter a spatial salt gradient or not (Fig. 2c). In the BE d zH increases over time, up to 400 nm. A similar trend was previously observed in literature, which validates our bulk experiments strategy [START_REF] Meng | Aggregation Rate and Fractal Dimension of Fullerene Nanoparticles via Simultaneous Multiangle Static and Dynamic Light Scattering Measurement[END_REF] . As expected, above the CCC, the nanoparticle surface charge is completely screened, thus eliminating the energy barrier to aggregation. Surprisingly, for roughly the same amount of salt at the outlet B of the microchip there is only a limited aggregate growth while we are well above the CCC, d zH increasing slightly from 170 nm and levels of around 260 nm. We do not observe any further growth after two hours. Two distinct aggregate size distributions (Fig. 2d) were obtained for Bout and 500mM for the BE. For the BE, the first population at 170 nm corresponds to the initial population of nC 60 , while the second population at 450 nm corresponds to the nC 60 aggregates induced by NaCl. This means that one part of the nC 60 aggregates that belongs to this initial population has not yet partake to the aggregation process. If the experiment would have last longer this first population would have disappeared progressively. In the MD experiment we also have two aggregate populations. However in the MD experiment, contrary to BE, the initial population is predominant in the size distribution. There is only a slight shift of the average aggregate size, from 170 to 250 nm. Therefore in this case the aggregation mechanism is quite different than before, for the BE. Moreover, as for A out , a micro-sized population is still observed but with much larger contribution to the overall population of B out . We suppose that there is a release of micron-sized aggregates of nC 60 formed in the MD. Nevertheless, their quantity is relatively low compare to the nano-sized population (170 nm) due to the dependence of intensity of light scattered (I) with size (r 6 ) and concentration c. Several hypotheses can be made to explain the difference between the particle aggregation in the bulk and the microfluidic experiment.

o First, the difference between microfluidics and the bulk experiments could be explained by the Diffusion-Limited Colloidal Aggregation (DLCA) and Reaction-Limited Colloidal Aggregation (RLCA) mechanisms [START_REF] Tang | Fractal Morphology and Breakage of DLCA and RLCA Aggregates[END_REF] . Indeed, it is well known that in batch mode, DLCA is the predominant behavior, where no repulsive barrier exists between the particles and the aggregate growth rate is entirely driven by Brownian motion. DLCA induces the formation of open structure aggregates, with fractal dimension (D f ) generally ranging between 1.5 and 2.0 with a linear increase of the aggregate size over time [START_REF] Chen | Aggregation and Deposition Kinetics of Fullerene (C60) Nanoparticles[END_REF][START_REF] Meng | Aggregation Rate and Fractal Dimension of Fullerene Nanoparticles via Simultaneous Multiangle Static and Dynamic Light Scattering Measurement[END_REF] , as we also obtained here (Fig. 2.c). In the other case, where there is a repulsive energy barrier between fullerene particles, RLCA is dominant with a limited aggregate growth, contrary to DLCA, as it was observed in case of the microfluidic experiment presented in this work. This result suggests that the salt gradient causes the nC 60 aggregation according to the RLCA mode and thereafter remains stable in time over our experimental window.

o Secondly, during the transportation of nC 60 through the microchip, after several hours of experiments and replicates larger aggregates can be formed in the salt front localized at the center of the channel corresponding to the black line in fig. 3. Indeed, this black-like deposition is characteristic of fullerene aggregates and is localized at the middle of the salt gradient. If too large aggregates could be formed, the micron-sized fraction of those sediment into the MD (Fig. 3). This hypothesis would confirm that micron-sized particles released in A out and B out are aggregates formed mainly within the salinity gradient in the MD. Nevertheless, this hypothesis cannot explain why the other part of the size distribution, small particles still present in B out at high [NaCl] (i.e. >CCC), do not form greater aggregates over time, as big as those observed in bulk experiments. Due to the NaCl concentrations used in this study, no representative electronic microscopy image of the nC 60 aggregates can be obtained. But, despite the different hypotheses described above, it appears clearly that nC 60 are transported and diffused through the salinity gradient.

During this short time transportation within the salinity gradient, a clear difference on the aggregation behavior is identified, compared to the classical results in bulk experiments.

Environmental implication

Based on these results and comparing to the literature concerning the environmental behavior of manufactured nanomaterials, several questions raised up. Indeed, both incidental and manufactured nanoparticles are expected to be released massively closed to urban system that will be eventually released in the river and finally transport through the estuarine system characterized by an increasing ionic strength. The expectations are that nanoparticles (NPs) will settle at different rate when the ionic strength increases. In literature, physical phenomena occurring in estuarine area, such as the salinity gradient and various flow conditions, were never considered. But, our results suggest that a short exposition of NP such as nC 60 within the salinity gradient affects definitely their aggregation mechanisms, and therefore their environmental behavior. Secondly, the scientific community on colloids and nanoparticles agrees to the fact that an ionic strength closed to the seawater ones leads to an irreversible aggregation and a change in size (up to micrometer). Our results suggest that the salinity gradient could enhanced the stability of nanoparticles, especially at high ionic strength (closed to the seawater one). This first study shows the need to reconsider the spatial and temporal variation of abiotic parameters (ionic strength, pH, organic matter concentration) in the understanding and the evaluation of nanoparticles environmental fate and transportation mechanisms.

Figure 1 :

 1 Figure 1: Scheme of the microfluidic dynamic and operation. At the entrance of the MD, two inlets are dedicated for the injection of nC 60 (nanoparticles) in DI water and for the NaCl with salinity closed to Estuary mouth (500 mM). Two outlets are identified, Aout and Bout, and are symmetric to the two entrances in term of geometry.

Figure 2 :

 2 Figure 2: (a,c) illustration of the z-average hydrodynamic diameter (nm), determined by the Cumulants algorithm, according to the ionic strength using bulk experiment and the use of MD (A out and B out ) (b,d) the average hydrodynamic diameter distribution determined by SBL algorithm and integrated along the whole duration of the kinetics (8500 and 4500 s for bulk an MD experiments, respectively). On (c) for 500 mM (BE), there are two linear fits of the data (d zH = a.T, with a=0.064 and a=0.023 for larger values).

Figure 3 :

 3 Figure 3: Illustration of the MD with deposition of nC 60 micron-sized aggregates in the center of the channel (dark lines) after several cycles of operation and replicates

Acknowledgements

This work was financially supported by the Interdisciplinary Mission (MI) of the French National Center for Scientific Research (CNRS; LaCUNE and LaCUNE-2 projects) and the Brittany Council of France (SAD program, NANOSALT project).