M Meo

U Tammaro

P Capuano

Influence of topography on ground deformation at Mt. Vesuvius (Italy) by finite element modelling

Ground deformations are observed in connection with volcanic activity, and therefore geodetic monitoring can provide significant indication of changes of equilibrium conditions. The aim of this paper is to study the deformation of Mount Vesuvius (Italy) caused by overpressure sources at various depths, using a commercial (Ansys) 3D finite element code, in the framework of linear elastic isotropic material behaviour. Both homogenous and heterogeneous media with carbonate basement were analyzed to understand the influence of topography on the ground deformations. The topography of the Somma-Vesuvius was taken into account, using a digital terrain model, and the carbonate basement was schematically modelled by assuming two horizontal layers with different Young moduli. The presence of a strong deviation from axially symmetric pattern of the displacement field, and of small subsidence areas, was found. These characteristics are completely unknown from the simple Mogi model and by simplified topography model, as verified by ad hoc simulations. These preliminary results, showing areas of the volcanic edifice experiencing high deformation, can improve the determination of the sources of deformations, i.e. the most relevant problem in the volcano monitoring. Moreover, the knowledge of the deformation pattern, including the topography effects, can provide significant indications to optimize the A c c e p t e d m a n u s c r i p t location of sensors and the characteristics needed to design an efficient and reliable geodetic monitoring network able to detect shallow intrusion events.

Introduction

During almost three centuries Mt. Vesuvius has been the natural laboratory where scientists came from all over the world to test theories and to experiment new methods. Mt. Vesuvius has been quiescent since 1944. Since then an uncontrolled expansion of the municipalities located on the slope of the volcano occurred; thus the dense population surrounding it and its slopes became the most notable aspect of Vesuvius among Earth's volcanoes. However, Vesuvius is a very dangerous volcano, in fact, its volcano's history suggests that the longer is the quiescence period, the more violent is the renewal of activity. This has transformed Mt. Vesuvius surroundings into an area with the highest volcanic risk in the world.

The Somma-Vesuvius (Fig. 1) is a volcano located within a graben (Campanian Plain) formed in Plio-Pleistocene. It is made by an older stratovolcano (Mt. Somma) with a summit caldera partially filled by the composite younger cone of Mt. Vesuvius. The Campanian Plain is bordered by mostly Mesozoic carbonate rocks. During the last 20 ky seven Plinian eruptions occurred. The last one is the famous "Pompeii" eruption which occurred in 79 A.D. After this eruption three sub-plinian eruption occurred (472, 512, 1631 A.D.). From 1638 to 1944 Mt.Vesuvius was in an almost continuous strombolian activity with occasionally moderate explosive eruptions. Since 1944, Mt. Vesuvius entered a quiescence stage, which is presently characterized by moderate low magnitude seismicity (up to M D = 3.6, recorded on October 1999), low temperature fumaroles at the summit crater and the absence of ground deformation.

For this reason, Mt. Vesuvius is carefully monitored. A 3-D model of the shallow structures of the volcano and an optimized monitoring network would allow a better definition of the characteristics of the precursors of eruptions, through a substantial improvement of the accuracy of location of seismic and ground deformation sources.

A c c e p t e d m a n u s c r i p t

Recently, strong efforts have been devoted to the definition of its structure and the location of magma reservoir, with the objective to help prediction of the scenario of the next eruption and to interpret the pattern of the expected precursory signals. A number of geophysical studies have been performed in the framework of a seismic tomography project, called TomoVes [Capuano et al., 2003]. The main results of this project are detailed images of the shallow volcanic structure and the Campanian Plain [START_REF] Zollo | Seismic evidence for a low velocity zone in the upper crust beneath Mount Vesuvius[END_REF][START_REF] Stefano | Active source tomography at Mt. Vesuvius: Constraints for the magmatic system[END_REF][START_REF] Natale | Seismicity and 3-D substructure at Somma-Vesuvius volcano: evidence for magma quenching[END_REF] and evidence for a 8 to 10 km deep midcrustal discontinuity associated with an extended molten or partially molten zone [START_REF] Auger | Seismic evidence of an extended magmatic sill under Mt. Vesuvius[END_REF]. This low velocity zone (at least 400 km 2 ) was interpreted as a large scale series of sills below the volcano edifice, with magma interspersed in a solid matrix, extending toward East and in the direction of the near volcanic centre of Campi Flegrei, covering a volume larger than 500k 3 . Analysis of fluid inclusion in phenocrysts provides information about the lithostatic pressure under which the Vesuvius magmas crystallize. Application of these methods to the last eruptions of Vesuvius (1906 and1944) suggests the presence of magma reservoir at shallow depth (2-5km), refilled during the eruptions by magma coming from a deeper (10-20km) reservoir [START_REF] Fulignati | New insight on the thermometamorphic-metasomatic magma chamber shell of the 1944 eruption of Vesuvius[END_REF]. Fluid inclusions in older Vesuvius products indicate crystallization depth between 4 and 10 km [START_REF] Belkin | Fluid inclusion studies of ejected nodules from plinian eruptions at Mt.Vesuvius[END_REF].

A 3D velocity model of the Mt. Vesuvius area have been reconstructed by [START_REF] Lomax | Precise, absolute earthquake location under Somma-Vesuvius volcano using a new three-dimensional velocity model[END_REF] to estimate earthquake location and fault mechanism. The main feature of this model is the complex shape of the carbonate basement surface.

In order to detect in advance significant variations in some parameters which can indicate reactivation, Vesuvius is monitored by a comprehensive system of seismic, geodetic and geochemical networks managed by the Osservatorio Vesuviano. One of the main problems in volcano monitoring is how design an optimal monitoring network to be able to detect weak precursory signals, generated by the interrelationship between sources, structures including topography and sensors. Thus, in order to investigate how the complex structure of a volcano can affect the deformation field, more realistic simulations are of vital importance to improve the The effects of topography on the deformation field have been presented in many papers

showing that the mechanical model of volcano [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF] were unsuccessful in the estimation of ground deformations because of the assumption of flat topography and half space solutions. In particular, [START_REF] Mctigue | Displacements and tilts from dip-slip faults and magma chambers beneath irregular surface topography[END_REF], investigated the effects of topography on surface displacements by using a plane-symmetric approximation. [START_REF] Cayol | Effects of topography on the interpretation of the deformation field of prominent volcanoes, Application to Etna[END_REF] considered an axial symmetrical volcano with average slopes of flanks ranging from 0° to 30° and they showed that the flattening of the deflation in the summit area is probably a topographic effect: the steeper the volcano, the lower the central uplift. [START_REF] Williams | The effects of topography on magma deformation models: application to Mt. Etna and radar interferometry[END_REF] showed that the predicted deformation field varies significantly when topography is included. [START_REF] Trasatti | effects of topography and rheological layering on ground deformation in volcanic regions[END_REF] investigated the effects of topography and rheological layering on ground deformation at Mount Etna, by assuming the approximate topography as axial symmetrical.

They showed, in agreement with results of [START_REF] Cayol | Effects of topography on the interpretation of the deformation field of prominent volcanoes, Application to Etna[END_REF], the existence, on the topography surface, of an annular region of maximum uplift with a central depression. [2004] developed 3D asymmetric and 2D axisymmetric finite element numerical models with heterogeneities of country rocks at Mt. Vesuvius. Their model took into account the presence of volcanic edifice: the volcanic structure was axially symmetric with a summit at about 1200 m above sea level and flanks drop to sea over 7 km.

Russo and Giberti

Finally, [START_REF] Lungarini | Finite element modelling of topographic effects on elastic ground deformation at Mt. Etna[END_REF] analyzed topographic effects on elastic ground deformation at Mount Etna with a model that use the real 3D topography of volcano. Their work showed, through the comparison of FEM and half-space models results, the importance of using a full 3D topography when studying ground deformation at volcanoes.

All the cited papers showed that real 3D topography can produce significant effects on the ground deformation pattern in the case of central volcanos with strong slopes of their flanks. In this paper, we calculated the ground deformation patterns of Vesuvius surface due to an overpressure

A c c e p t e d m a n u s c r i p t

source at various depths, by developing for the first time a full 3D finite element model based on the real topography. This model can be used for the detection of the sources of deformations and for the design of an optimal geodetic monitoring network.

Seismic tomography investigations

The idea of realizing a high resolution seismic experiment on Mt. Vesuvius was driven by the need of a detailed reconstruction of the volcano structure, including the location and size of the magmatic reservoir(s), in order to predict quantitatively the possible scenarios of the next eruption.

Although a hundred local microseismic (M D < 3.6) events occur in average every year at Mt.

Vesuvius, most of them are located in a clustered area just below the crater and they are shallower than 4 km making the passive tomography approach unfruitful. Then, the proposed project relied mostly on active energy sources.

During TomoVes experiment, five multi-two-dimensional profiles (Fig. 1) were analysed with a total of 17 shot points, recorded by about 140 3-component seismographs. Seismic signals of each shot were recorded simultaneously along the longitudinal and quasi-orthogonal directions (for more details see [START_REF] Gasparini | Looking inside Mount Vesuvius[END_REF].

Firstly, a 1D non linear inversion τ-p method approximating the Earth's topography with spherical shells and including body waves was performed [START_REF] De Matteis | 1-D P-velocity models of Mt Vesuvius Volcano from the Inversion of TomoVes96 First Arrival Time Data[END_REF] to provide a reliable 2D initial model for perturbative tomography inversion. The choice of a non-linear approach, instead of a linearized tomography inversion was due to the strong heterogeneity of the investigated medium that make the definition of a reliable initial velocity model challenging, despite of the availability of subsurface information. In this case, a linearized approach may fail if the used a priori reference velocity model is too far from the actual velocity structure.

In particular, a Bayesian approach [START_REF] Zollo | Bayesian estimation of 2-D P-velocity models from active seismic arrival time data: Imaging of the shallow structure of Mt. Vesuvius (Southern Italy)[END_REF] for mapping errors on parameters and for selecting the optimal medium parameterization. The latter was evaluated by applying the Akaike Information Criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF]) which takes into account the concepts of best misfit and simplicity of the model (small number of parameters). The search for the maximum likelihood model, maximizing the pdf, in the multidimensional parameter space was performed by a non-linear optimization method which uses, at consecutive steps, the Genetic [START_REF] Goldberg | Genetic algorithm, in Search, optimization and machine learning[END_REF]) and Simplex [START_REF] Press | Numerical recipes[END_REF]) algorithms. A 'multi-scale' approach was used. This procedure, which was introduced for velocity estimation by [START_REF] Lutter | Inversion of crustal structure using reflections from the ASSCAL Ouachita experiment[END_REF], allowed us to determine the large-scale components of the velocity model and to estimate progressively the smaller-scale components. This is roughly equivalent to moving from a low to high wavenumber description of the velocity field. The cubic spline was used to compute the velocity values at every point of the 2D gridded medium. The cost function used in the inversion was a least-squares L 2 norm, defined as the sum of the weighted squared differences between observed and computed travel times.

Fig. 2 shows 2D tomography images relative to the 5 analyzed profiles. The inferred images show a high variability in the shallow structure, with values ranging from 1.5 km/s to 6 km/s, with a very sharp increase observed at depth ranging from few hundreds meters to 3-4 km. It was reasonably assumed that the top of the limestone formation is marked by the P-velocity isoline at 5000 m/s. The checkerboard analysis of the retrieved images shows that the best resolved areas are located beneath the volcanic complex.

Using a circumferential interpolation procedure a 3D P-velocity model and the topography of the top of the limestone basement were computed [START_REF] Lomax | Precise, absolute earthquake location under Somma-Vesuvius volcano using a new three-dimensional velocity model[END_REF]. This model was used to refine earthquake location occurred in the Vesuvius area, using a non-linear global search location, using the probabilistic inversion presented by [START_REF] Tarantola | Inverse problems=quest for information[END_REF].

A c c e p t e d m a n u s c r i p t

Finally, tomography studies (Capuano et al., 2003) revealed the absence of low velocity zone compatible with magma reservoirs up to 5-6 km depth, with a lateral resolution of approximately 0.5km. Based on this results and with the findings of a low velocity zone at 8-10km depth, as the possible primary source of the shallower magma chamber, we simulated overpressure source located at 2, 3, 5 km bsl with a lateral extension smaller than 0.5km, using reservoirs with a parallelepiped shape.

FEM model

A 3D linear-elastic finite element model to study ground deformation at Mt. Somma-Vesuvius was developed, using ANSYS commercial code. The validity of this approach is related to the time scale of the event we want to simulate. In general, the time scale of overpressure in magma chamber is much smaller than the scale of material relaxation in a volcano. Our model was built taking into account real topography of the volcano as inferred from a detailed digital elevation model. When modelling the behaviour of large scale structures such a volcano, various modelling issue must be taken into account. To describe a volume, either hexahedral or tetrahedral shaped can often be used; however, a combination of the two elements is not generally recommended. In our studies the choice of using hexahedral instead tetrahedral elements was dictated by higher accuracy and performance of these elements. The hexahedral element used was defined by eight nodes having three degrees of freedom at each node: translations in the all nodal directions.

One of the most important modelling issues is the assessment of the appropriate mesh density. Using a coarse mesh can generate inaccurate results. On the other hand, using a fine mesh leads to an increase of run times. To avoid such problems, in our work a sensitivity analysis of the appropriate mesh density was performed.

Another issue for a correct analysis was the level of details necessary to describe adequately the structural behaviour of the volcano. The mesh was constructed by using elements with smaller edge length in regions of high stress gradient while a coarser mesh was used away from the main The size of the domain is clearly a function of the loading conditions. In our case, various dimensions were considered for the various loadings considered to avoid edge effect. The domain size was chosen in such a way that no significant changes of the vertical displacement in volcano area were recorded for the loads imposed. Domain size was increased until nearly identical results for succeeding meshes were obtained. A tolerance of 5% was chosen between two different domain sizes. Finally, the model size is 37000x26000 m with a depth size of 25000 m (Fig. 3).

Regarding the boundary conditions, different constraints were imposed on the faces of the domain (Tab. 1).

A Pre-conditioned Conjugate Gradient iterative equation solver was used to perform the linear static analysis due to the reduced disk file space needed compared to other solvers and for its speed in solving large scale models. Sensitivity analysis of solver tolerance value was conducted and a value of 1.0•10 -6 was used.

Our objective is the computation, using a FEM approach, of the ground deformation at Vesuvius volcano, evidencing the strong influence on the deformation pattern due to the real complex topographic structure. Three different overpressure sources were considered located at various depths, 2, 3 and 5 km in a homogenous medium and in heterogeneous medium with carbonate basement. In particular, two half-space models to compute deformation at surface were used: a) an homogeneous model where the reservoir was represented by a prolate overpressure source parallelepiped shaped with a horizontal side of 480 m and thickness of 1.8 and 2.0 km, for the shallow and deeper source respectively, locate along the crater axis; b) an heterogeneous model schematically obtained assuming two horizontal different layers, the first underneath the crater and down to 2 km below sea level with Young's modulus equal to 10 GPa and the second layer,

A c c e p t e d m a n u s c r i p t

representing the carbonate basement, with a Young's modulus of 40 GPa [START_REF] Russo | Numerical modeling of surface deformations on Mt. Vesuvius volcano (Italy) in presence of asymmetric elastic heterogeneities[END_REF].

The Poisson's ratio was assumed to be equal to 0.3. The morphology of at the interface, located at 2 km depth [START_REF] Zollo | Seismic evidence for a low velocity zone in the upper crust beneath Mount Vesuvius[END_REF] was considered flat in this work. The source geometry was chosen taking into account the tomography images and its resolution [Capuano et al., 2003], as explained in the previous section. To obtain an adequate source volume we used a prolate parallelepiped shape, differently from [START_REF] Russo | Numerical modeling of surface deformations on Mt. Vesuvius volcano (Italy) in presence of asymmetric elastic heterogeneities[END_REF] that used prolate ellipsoidal or cylindrical shape.

The shape size was also constrained by the number of elements and by the type of element (hexahedral).

For each model the three components of the ground deformation were computed at the nodes of FEM mesh, corresponding to the real topographic surface. A uniform overpressure of 10

MPa was applied to the source walls. Different overpressure values will affect the ground deformation magnitude but the deformation pattern will keep unchanged, due to the linear elastic material model used in the simulation.

Results and conclusion

Results relative to the source located at 2, 3 and 5 km are reported below for homogenous and heterogeneous medium.

Homogenous medium

The surface deformation produced by a source located at 2 km depth is shown in Figure 4a.

A strong deviation from an axial symmetrical pattern of the displacement field is evident due to the presence of the real topography. As pointed out by [START_REF] Lungarini | Finite element modelling of topographic effects on elastic ground deformation at Mt. Etna[END_REF] at Etna volcano, the peculiar effect of real topography is to increase displacements at the point of lower elevation.

Maxima are located on the slope of the volcano (Colle Umberto, Colle dei Canteroni and Cognoli di Trocchia area), while the maximum uplift is located in SE area. Deformation is concentrated in a radial distance of about 3 km from the crater axis. Another relevant characteristic is the presence of small area in subsidence at middle-high elevation, with a strong variability in the summit area. The

A c c e p t e d m a n u s c r i p t

complex deformation pattern is still present when the source is located deeper, even if with a smaller amplitude. This is shown in Fig. 5a,b, where the deformation for each component and each source depth, along N-S and E-W profiles are plotted. To verify the effects of real topography we firstly computed the ground deformation produced by a Mogi model without topography and then by assuming an axial symmetrical shape of the volcano. The axial-symmetrical volcano was realized by rotating a simulated profile obtained considering a height of 1230m (maximum elevation of the volcano) and a constant slope to reach sea level at a distance of 5km (distance of the crater axis to the coast); the sources were located 3 km depth below sea level and their volume was equivalent (0.4 km 3 about).

Figure 6 shows the vertical and radial ground displacements for both simulations showing the complexity added by the real topography, which obviously decrease with the source depth. The Mogi model and also the simplified topography were not able to define the influence of the real topography and therefore they can provide incorrect localisation of the sources of the observed deformations and incorrect data to signal a possible volcano reawakening. The real topography produced also non-zero values of the NS component along EW direction and vice versa (Fig. 5a,b).

The horizontal components of the deformation field show a complex pattern in the summit area, evidencing also some shortening areas.

Heterogeneous medium

The surface deformation produced by a source located at 2 km depth is shown in Figure 4b.

Also in this case with the presence of a carbonate layer, the ground deformation is strongly affected by the topography with similar characteristics. Moreover, maximum vertical displacements are greater with respect to the homogeneous case in the crater region. In this case the source was located above and below the interface between the two layers (carbonate and volcanic coverage) and the differences in Young's modulus creates additional vertical forces in the area above the Sources located at greater depth showed similar behaviour with smaller amplitude. The maximum horizontal displacements show similar increased amplitude and coupling between the two horizontal components (Fig. 5c,d).

Conclusion

In a volcanic area, ground deformations are generally studied by using a simple models, taking into account both the topography and the heterogeneity of the medium. This approach can provide incorrect measured displacements and the deformation scenario when used in hazard analysis. Thus, simulation of the expected deformation pattern with a more complex model was needed. We started to investigate the influence of real topography at the most hazardous central volcano in the world: the Somma-Vesuvius structure. We performed simulations of overpressure sources whose geometry (shape, depth and lateral extension) were constrained mainly by results of seismic tomography studies carried out in the last ten years.

The main results of our simulation, performed using a real topography, was the presence of a strong deviation from the axial symmetrical deformation pattern, around the crater axis, with the peculiar effect of enhancing displacements at the point of lower elevation. This behaviour was also characterized by small subsidence areas and by the presence of non-zero horizontal component (x or y) along the orthogonal direction (y or x). Both these characteristics were completely unknown in the simple Mogi model and also by model created using simplified topography. These results can be very helpful in the understanding the position of the sources of ground deformation during volcano reawakening and for a better definition of the geometry and characteristic of a geodetic monitoring network, particularly in the summit area, that experienced strong deformation differences with respect to the simplified models.

A c c e p t e d m a n u s c r i p t

For a more accurate modelling that could help in unbiased analysis of deformation produced by possible uplift events on the summit of Vesuvius, and in the definition of a more accurate monitoring network, further refinements of the model are needed. This will include a finer mesh, the introduction of the real morphology of a carbonate basement as shown by the tomography [START_REF] Zollo | Seismic evidence for a low velocity zone in the upper crust beneath Mount Vesuvius[END_REF], a detailed description of the elastic characteristics of the volcanic coverage for an improved modelling of a complex elastic structure, the introduction of the bathymetry of the gulf of Naples, and finally the simulation of different pressure source shapes, compatible with the tomography constraints.

  sources of surface deformations and to find optimal measurement sites location.

  to image the shallow structure of Mt Vesuvius volcano and the surrounding plain was developed. The problem of velocity field estimation was solved through the computation, during the model space exploration, of the a A c c e p t e d m a n u s c r i p t posteriori probability density function (pdf) which relates the observed to predicted arrival time data. In this way, a maximum likelihood velocity model was obtained, which is particularly useful

  the main area of interest (with a side of 10 km) enclosing the volcano, brick elements with an edge length of 80 m were used; then a coarser mesh was used for the remaining domain with elements having a decreasing (maximum) edge length of 240 and 480 m. The final model was composed of 643525 nodes and 612880 elements.

  regions away from the crater axis, higher displacements were experienced by the model with homogenous material.

Fig. 1

 1 Figure captions and Table

Fig. 2

 2 Fig. 2 2D P-velocity model relative to 5 profiles. The dashed line delimitates at depth the best

Fig. 5

 5 Fig. 5 Normalized vertical and horizontal displacements profile with source located at 2 (black), 3

Fig. 6

 6 Fig. 6 Vertical (red) and radial (black) displacements for the Mogi model (dashed line) and an

Figure 1 A

 1 Figure 1

  Figure 5a

  Figure 5c

Acknowledgments

We are grateful to the Laboratorio di Geomatica e Cartografia (INGV-OV) for providing the DTM of the area. This work was partially funded by Istituto Nazionale di Geofisica e Vulcanologia and Dipartimento della Protezione Civile. The reviewer's comments improved strongly this paper.


Guillaume Renaud
email: guillaume.renaud28@etu.univ-tours.fr

Samuel Callé

Jean-Pierre Remenieras

Marielle Defontaine

Nonlinear acoustic

Keywords: Acoustic waves interactions, Harmonic generation, Nonlinear elasticity, Microdamage, Bone PACS: 43.25.Ba, 43.25.Zx, 43.80.Qf

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

Introduction

It is possible in the skeleton to distinguish between cortical bone with low porosity (less than 30%) and trabecular bone with high porosity (more than 70%). As most bone fractures in the elderly occur in trabecular bone, this study focused on this type of bone, particularly the heel bone. Trabecular bone is a spongy bi-phasic material with connected porosity. Its solid phase consists of a connected network of bone tissue rods known as trabeculae. Porous spaces are filled with marrow (mostly fat). The clinical diagnosis of bone fragility is currently performed using X-ray measurements which provide imaging of bone mineral density. This quantitative assessment, similar

A c c e p t e d m a n u s c r i p t

to porosity measurement, accounts for 60-70% of variation in bone strength [START_REF] Ammann | Bone strength and its determinants[END_REF]. However measuring the amount of material in a porous structure is not enough to evaluate its strength. It also depends on the spatial arrangement of the material and on the mechanical characteristics of the bone tissue. Furthermore, it is well known that cracks typically of tens to hundreds of µm in size are generated in bone tissue by daily loading. Figure 1 shows different types of common bone microdamage. In healthy bone, damage is repaired by continuous bone remodeling. When an imbalance in this process occurs, cracks can accumulate. Crack surface density, measured by histological analysis using staining techniques [START_REF] Burr | Validity of the bulk-staining technique to separate artifactual from in vivo bone microdamage[END_REF], and the associated standard deviation increase with age [START_REF] Fazzalari | Cancellous Bone Microdamage in the Proximal Femur: Influence of Age and Osteoarthritis on Damage Morphology and Regional Distribution[END_REF] [13] [START_REF] Zioupos | Ageing human bone: factors affecting its biomechanical properties and the role of collagen[END_REF].

Whereas crack surface density is below 1 crack per mm 2 in young healthy people, it can reach values above 10 cracks per mm 2 in elderly people. Whether the increase in crack density is a cause or a consequence of the decrease in the bone toughness properties with age is still not resolved [START_REF] Burr | Bone microdamage and skeletal fragility in osteoporotic and stress fractures[END_REF] [23] [START_REF] Zioupos | Accumulation of in vivo fatigue microdamage and its relation to biomechanical properties in ageing human cortical bone[END_REF].

Linear ultrasound imaging techniques have been developed to image trabecular bone velocity and attenuation [START_REF] Gomez | In vivo performance of a matrix-based quantitative ultrasound imaging device dedicated to calcaneus investigation[END_REF]. These techniques are nonetheless almost insensitive to progressive induced damage [START_REF] Nicholson | Quantitative ultrasound does not reflect mechanically induced damage in human cancellous bone[END_REF]. On the other hand, it is now well established that inhomogeneous materials, with soft mesoscopic flaws embedded in a rigid matrix, exhibit strong nonlinear quasi-static and dynamic behavior known as non-classical or hysteretic nonlinearity [START_REF] Guyer | Nonlinear mesoscopic elasticity : evidence for a new class of materials[END_REF] [10] [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF]. Non-destructive evaluation techniques are being developed to detect localized damage in materials [START_REF] Van Den Abeele | Microdamage diagnostics using nonlinear elastic wave spectroscopy (NEWS)[END_REF] [27] [START_REF] Zaitsev | Novel nonlinearmodulation acoustic technique for crack detection[END_REF]. Nevertheless very few studies have been performed on the relationship between crack density and level of nonlinearity , where ρ 0 is the density of the medium and p the acoustic pressure) was found to range between 80-120 in bovine trabecular bone [START_REF] Donskoy | Nonlinear acoustic parameter of trabecular bone[END_REF], whereas β does not exceed 7.5 (PMMA) in homogeneous material [START_REF] Johnson | Nonlinear tube waves[END_REF]. Hysteretic nonlinearity has been clearly observed in bovine and human cortical bone using nonlinear resonance ultrasound spectroscopy [START_REF] Muller | Mise au point d'une technique ultrasonore d'acoustique non linaire pour l'etude de l'endommagement osseux (Development of a nonlinear ultrasound technique to assess bone damage[END_REF] [START_REF] Muller | Nonlinear resonant ultrasound spectroscopy (NRUS) applied to damage assessment in bone[END_REF]. One interesting point is that the nonlinear hysteretic parameter α, related to the resonance frequency downshift, increases overall with age as well as its dispersion. This trend is similar to the increase in crack density with aging in trabecular [START_REF] Fazzalari | Cancellous Bone Microdamage in the Proximal Femur: Influence of Age and Osteoarthritis on Damage Morphology and Regional Distribution[END_REF] and cortical bone [START_REF] Schaffler | Aging and Matrix Microdamage Accumulation in Human Compact Bone[END_REF]. Crack density is probably a mirror of the "biomechanical health" of bone tissue and the increased dispersion with age probably reveals disparities in nutrition, physical activity and diseases. Several features have been identified to improve understanding of the relationship between the level of nonlinearity and crack density. These include the effects of preferential orientation of cracks and the relative nonlinear contributions of different types of damage in bone.

We propose three different modalities to measure trabecular bone nonlinearity in this context in order to image crack density in future studies. Such noninvasive methods would be valuable for in vitro biomechanical studies, to understand how microcracks impact on the biomechanical properties of bone and for in vivo monitoring of crack accumulation (to date only possible by biopsy). The following sections comprise a brief technical presentation and experimental results obtained from harmonic generation, parametric emission and parametric reception using phase modulation.

Method 1 : Harmonic generation

One way to investigate the nonlinearity exhibited by a material consists of increasing monofrequency excitation and observing the dependence of harmonic amplitudes on the fundamental amplitude. According to theory, in the case of weak nonlinear effects classical nonlinear materials have n th harmonic amplitude increasing with power n of the fundamental amplitude. Typical manifestation of microcracks in damaged solids is the dominance of odd harmonics and their almost quadratic dependence on pump amplitude [START_REF] Van Den Abeele | Nonlinear Elastic Wave Spectroscopy (NEWS) Techniques to Discern Material Damage, Part I: Nonlinear Wave Modulation Spectroscopy (NWMS)[END_REF].

Experimental method

Figure 2 describes the experimental set-up. We used a focused transducer to generate a monofrequency wave at 400 kHz in a water tank. The width of the focused beam at -6 dB was less than 1 cm. The amplitude of the transducer

A c c e p t e d m a n u s c r i p t

excitation was increased from 30 to 160 V peak-peak. The fundamental and harmonic pressures were measured by a needle hydrophone placed at 45 mm along the axis of the transducer, where the fundamental, second and third harmonics reach their maximum levels in water (figure 3). The experiment was first performed in water without any sample. Marrow was then removed from the bone sample and the experiment was conducted with a slice of defatted heel bone (24 mm thickness) and a slice of porous ceramic (20 mm thickness) (see Table 1 for more information on the samples). The samples were placed about 20 mm from the transducer so that the generation of harmonics occurred only in the sample (figure 3). Because of the strong attenuation in both trabecular bone and porous ceramic, the resonance effects in the sample were expected to be very weak. Moreover, the reflection coefficient for acoustic pressure between trabecular bone (considered as an effective medium) and water is 0.1.

A c c e p t e d m a n u s c r i p t 2.2 Results and interpretation

The pressures of the fundamental, second, third and fifth harmonics measured in water, through trabecular bone and through porous ceramic are plotted as functions of the transducer voltage excitation in figures 4, 5 and 6, respectively. sample) was 112 kPa, thus the hypothesis of weak nonlinear effects is valid. The results obtained in water clearly showed classical nonlinear behaviors : second, third and fifth harmonics behaved with slopes 2.08, 3.08 and 4.34, close to the theoretical values of 2, 3 and 5, respectively. In contrast, trabecular bone and porous ceramic exhibited hysteretic nonlinearity characteristics : amplitude of the third harmonic higher than that of the second harmonic (despite very high attenuation due to scattering above 800 kHz, see Table 1) and a power law close to 2, especially for the third and fifth harmonics. In trabecular bone the second harmonic amplitude increased with power 3 of the fundamental amplitude. This astonishing behavior was observed in damaged materials. Some studies have reported that progressively induced damage increases the power law relating the second harmonic with the fundamental amplitudes, whereas the overall level of the third harmonic was increased, without impacting on the power law [START_REF] Van Den Abeele | Microdamage diagnostics using nonlinear elastic wave spectroscopy (NEWS)[END_REF]. The common velocity dispersion measured in trabecular bone is -26 m/s/MHz in the frequency range 200-800 kHz [START_REF] Wear | Measurements of phase velocity and group velocity in human calcaneus[END_REF]. In view of the small thickness of the samples, dispersion effects may have little effect on harmonic generation. Moreover harmonic distortion can be induced by the pressure-dependence of the flow resistivity above a certain acoustic velocity threshold. After a certain propagation distance, odd harmonics have higher amplitudes than even harmonics. For example, Kevlar with porosity of 96% exhibits that behavior when acoustic velocity exceeds 1 m/s [START_REF] Nelson | Harmonic generation, propagation and attenuation for finite-amplitude tones in an air-filled porous material[END_REF]. In the present study, the acoustic velocity did not exceed 8 cm/s and the mean pore diameter was larger in trabecular bone (0.5-1 mm) than in Kevlar. We therefore assumed that this effect was negligible in our experiments.

A c c e p t e d m a n u s c r i p t 3 Method 2 : Parametric emission

An alternative way to investigate nonlinear characteristics of bone is parametric generation. We propose a method based on the combination of two focused primary beams emitted at two slightly different frequencies (f 1 and f 2 ) using an annular array [START_REF] Fatemi | Ultrasound-stimulated vibro-acoustic spectrography[END_REF]. These primary waves, focused in the bone, create an amplitude modulation giving rise to a low-frequency (LF) wave at f LF (equivalent to the frequency difference of the two primary frequencies). An internal stress field results from this LF excitation, and the acoustic response of the object to this mechanical excitation is used to form the image. The LF acoustic field resulting from the vibration of the medium is detected by a sensitive hydrophone and then used to form the image. Two main physical phenomena are responsible for this LF wave generation: the acoustic radiation force (due to reflection and attenuation) [START_REF] Fatemi | Ultrasound-stimulated vibro-acoustic spectrography[END_REF] and the nonlinear interaction [START_REF] Callé | Presence of nonlinear interference effects as a source of low frequency excitation force in vibroacoustography[END_REF]. The overall low frequency field q LF can be expressed as q LF = q rad + q NL , where q rad represents the LF pressure originating from the radiation force and q NL represents the contribution of the LF nonlinear interaction (parametric generation).

Experimental configuration

In a previous study [START_REF] Remenieras | A simple calculation approach for the parametric sound field generated by a focused annular array[END_REF] we reported that the choice of the focused annular array depends on several parameters (e.g. the curvature of the annular array, the frequency downshift ratio f 1 /f LF and the number of rings). Finally, a 1 MHz array with 4 rings and a focusing gain (G) of 16.14 per element was an acceptable compromise between sensitivity and resolution for parametric generation [START_REF] Callé | Presence of nonlinear interference effects as a source of low frequency excitation force in vibroacoustography[END_REF]. The two odd rings emit at f 1 and the two even rings emit at f 2 . The focal distance and the probe diameter were 6 cm and 8.32 cm, respectively. In this configuration, the LF field resolution at the focus is in the same range as the primary field resolutions [START_REF] Callé | Application of nonlinear phenomena induced by focused ultrasound to bone imaging[END_REF]. The LF resolution (1 mm) corresponds to the lateral resolution of our system, which is promising for image construction. Figure 7 presents the experimental set-up: the bone slice faces the transducer, and the LF hydrophone linked to the transducer is beside the object.

In vitro heel bone imaging

Due to both the overall resonance frequencies and the local bone properties, the LF pressure amplitude does not vary simply according to the choice of f LF [START_REF] Callé | Application of nonlinear phenomena induced by focused ultrasound to bone imaging[END_REF]. The optimal frequency of 76 kHz was therefore chosen for the later heel bone images. Figure 8 represents the image of the LF pressure amplitude in the medium plane of a 19 mm thick human heel bone slice (figure 9). The spatial resolutions related to displacement steps were 300 µm in both directions. The central frequency was 800 kHz to minimize attenuation effects in the bone. The image clearly shows the trabeculae framework. The cortical outlines of the heel are also clearly demonstrated. Moreover, we can see that the white areas of the image (large parametric amplitude) at the bottom of the heel bone correspond to the region in contact with the ground, which undergoes daily loading. Such areas could contain many microcracks.

Interpretation and discussion

The radiation force resulting from the combined primary fields at frequencies f 1 and f 2 is one source of LF stress in the object [START_REF] Fatemi | Ultrasound-stimulated vibro-acoustic spectrography[END_REF]. However, we have already reported that a LF field is created by the nonlinear interference of the two primary beams inside the propagating medium [START_REF] Callé | Presence of nonlinear interference effects as a source of low frequency excitation force in vibroacoustography[END_REF]. This field, directly related to the coefficient of nonlinearity β, may then propagate, scatter and be reflected on the LF contrast sources of the object and thus provide a measurable radiating field. In the above study the LF acoustic response of tissue varied with tissue type and with the choice of the difference frequency, indicating that this is potentially a valuable tool for bone investigation. The

A c c e p t e d m a n u s c r i p t

main advantages of this method compared to other ultrasound bone imaging methods are: (i) the non-transmission mode, which for instance allows in vivo imaging of the neck of femur, and (ii) its high lateral resolution. Moreover, we reported that the creation of this internal force in our conditions depended on the attenuation and impedance (radiation force term) and the nonlinear coefficient (nonlinear interference term). In practice, the contribution of the preponderant phenomenon (radiation force or NL interaction) has not been clearly established, since it strongly depends on the type of tissue being studied and is particularly difficult to estimate in bone. Thus, although the LF signal is partly linked to the coefficient of nonlinearity β and hence is larger for bone regions affected by microcracks, this technique does not exclusively provide information about microdamage.

Method 3 : Parametric reception using phase modulation

The third modality proposed to measure nonlinearity in trabecular bone is parametric reception using a phase modulation method.

Experimental method

In this method a low-frequency wave (f LF = 150 Hz) and a high-frequency wave (f HF = 800 kHz) are transmitted in the medium in perpendicular directions and cross in the sample immersed in water (figure 10). The LF wave is generated by an 8 cm diameter circular piston fixed to a shaker while the HF wave is emitted by a plane transducer. Nonlinear interactions in the water and in the sample generate two parametric waves at the sum and difference of the two primary frequencies f LF and f HF . If those parametric frequencies belong to the bandwidth of the HF receiver, this can be interpreted as a modulation of the HF wave phase. The propagation velocity of the HF wave in the medium is modulated by the LF wave due to nonlinearity.

Because of the high primary frequency ratio f HF /f LF , and if the interaction distance is much less than the LF wavelength, the LF acoustic disturbance appears to be quasi-static for the HF wave. Considering two plane waves interacting along the distance L in a homogeneous medium, the phase modulation is given by integration of the wave number along the interaction distance :

φ(t) = 2πf HF Lβ c 2 v LF (t),
where c, β and v LF are the propagation velocity, the classical quadratic nonlinear parameter and the LF acoustic velocity, respectively. β is defined as , where ρ 0 is the density of the medium and p the acoustic pressure. Thus, according to classical nonlinear theory, the phase modulation has the same temporal form as the LF acoustic velocity. If we insert a sample (as presented on figure 10) of different density, velocity and nonlinear parameter in a water tank, the phase modulation is now expressed by :

φ(t) = 2πf HF (L -L sample )β water c 2 water v LF water (t) + 2πf HF L sample β sample c 2 sample T water-sample v LF water (t),
where T water-sample is the transmission coefficient for acoustic velocity given by :

T water-sample = 2Z water Z water + Z sample ,
where Z water and Z sample are the acoustic impedance of water and the sample, respectively. The experiment consisted of measuring the increase in phase modulation as a function of the LF wave amplitude. Repeating the protocol in water (without sample) and with a sample provided the ratio β sample /β water . Because the LF pressure amplitude is proportional to the inverse of the distance from the piston, the sample was placed as close as possible to the piston and the HF beam axis was thus 2.5 cm below the piston. The maximum LF pressure amplitude measured in water, corresponding to a 300 mV peak-peak excitation of the transducer, was 2000 Pa in the cross-section of the LF and HF beams. The pressure under the LF piston was nearly uniform at 2.5 cm under the piston, except near the edges. The HF wave was received by a plane transducer after propagation through the sample. The phase modulation amplitude was obtained by numerical analysis. The signal was band-pass filtered around 800 kHz to suppress noise, then the phase of the received signal was

A c c e p t e d m a n u s c r i p t

calculated by Hilbert transform and phase unwrapping. Finally the slope of ωt was subtracted and the modulation contribution remained. The mean amplitude of the phase modulation was calculated using Fast Fourier Transform analysis. The phase noise level was below 0.5 mrad.

Results

By performing the experiment in water, i.e. without a sample, and with a sample placed in the cross-section of the two beams, it is possible to evaluate the nonlinear parameter of the sample. The phase modulation amplitude is reported as a function of the LF wave amplitude on figure 11. We found β = 14 for Plexiglas (PMMA), whereas the usual value is 7.5. This over-estimation was probably due to reflections inside the Plexiglas sample. Attenuation was weak (Table 1) and the reflection coefficient for acoustic velocity was 0.38. Cumulative effects can occur and enhance the phase modulation. These effects may be negligible in trabecular bone and porous ceramic because of strong attenuation. Considering the bone and ceramic samples as effective media, we found β values of 36 and 3400, respectively. As the nonlinear effects happened mainly in the solid phase, the corresponding β had to be estimated using the solid phase properties instead of the effective media properties (T water-sample , c sample and L sample ) (see Table 1). We then obtained β values of 440 and 24000 for trabecular bone and porous ceramic, respectively. Porous ceramic exhibited a much higher nonlinear param-

A c c e p t e d m a n u s c r i p t

eter than trabecular bone. Bone tissue comprises 70% calcium hydroxyapatite crystals and 30% organic collagen fibers (in volume) whereas porous ceramic comprises only industrially made calcium hydroxyapatite crystals. This material used in surgery may not contain microcracks. It is composed of juxtaposed grains (like rocks) of typical size in the order of 1 µm. For higher LF wave amplitudes (above 2000 Pa) phase modulation increased with different nonlinear behaviors that require further investigation.

Discussion and conclusion

The aim of this preliminary study was to define the best method(s) to measure and scan trabecular bone acoustic nonlinearity. The measurement of second and third harmonic amplitudes as functions of the fundamental amplitude provided valuable information and the second harmonic slope and the level of third harmonic may provide qualitative and/or quantitative information about microdamage. The fundamental frequency must not exceed 400 kHz to limit attenuation and velocity dispersion. Nevertheless, this technique needs fairly large fundamental amplitudes, especially for in vivo application. Soft tissue and the dense cortical shell of the heel bone (which was removed in this study) add an attenuation of 10 to 15 dB before the wave reaches the trabecular bone area, meaning that three to six times more than the 112 kPa

A c c e p t e d m a n u s c r i p t

used in the present study is necessary. To avoid biological thermal effects, the acoustic pressure amplitude must not exceed 0.5 to 1 MPa for a long 400 kHz burst (tens of milliseconds), although the threshold was in fact not clearly established. Moreover, harmonics are also attenuated by 10 to 15 dB when propagating through the cortical shell and soft tissue, i.e. from trabecular bone to water surrounding the heel. This means that the fundamental pressure amplitude required (in water before propagating in the heel) would be around 1 MPa. Although this configuration is not appropriate for in vivo application, it is a simple experimental set-up for in vitro studies.

The parametric generation configuration presented the same drawback. Resonance effects can be used to enhance the level of the parametric wave (as in the present parametric emission experiment). Nonetheless, evaluation of the gain due to resonance for each point is necessary to interpret two different values measured in two regions of the bone. Perhaps an ad hoc parametric frequency should be chosen for each point of the scan. The maximum variation in ultrasound velocity measured in a heel bone is approximately 300 m/s (typical velocities are in the range 1700-2000 m/s). However, we previously reported that in whole heel bone the surrounding cortical surfaces directly influence propagation of the primary wave and provide an image that is more related to relief than to the internal structures [START_REF] Callé | Application of nonlinear phenomena induced by focused ultrasound to bone imaging[END_REF]. Several effects related to reflection and attenuation on complicated bone shape make the images smoother and more difficult to interpret. However, it is currently difficult to prevent these effects.

In contrast to the other two methods, parametric reception using the phase modulation method does not need a high level of energy. Measurement of the phase modulation at a given LF pressure showed a small standard deviation, but the accuracy and reproducibility of the parametric reception method require evaluation. Moreover, the nonlinear behaviors of the increase in the phase modulation need to be investigated for higher LF pressures. Measurement of the second-order phase modulation induced by parametric frequencies f HF ±2f LF must also be studied. The amplitude of the second-order phase modulation is proportional to the hysteretic nonlinear parameter α [START_REF] Van Den Abeele | Nonlinear Elastic Wave Spectroscopy (NEWS) Techniques to Discern Material Damage, Part I: Nonlinear Wave Modulation Spectroscopy (NWMS)[END_REF] and may provide different information on bone microdamage. Impact can be used instead of a LF wave and this was tested for nondestructive evaluation applications [START_REF] Van Den Abeele | Microdamage diagnostics using nonlinear elastic wave spectroscopy (NEWS)[END_REF]. Further study is also needed to establish which of the two methods is the most accurate. Parametric reception using phase modulation is thus the best configuration to measure in vivo bone acoustic nonlinearity in order to assess crack density. Phase modulation will be measured using focused rather than plane HF transducers to study spatial variations. The relationship between levels of acoustic nonlinearity and crack density needs to be validated by histological analysis.

Fig. 1 .

 1 Fig. 1. Four types of damage in trabecular bone (grey areas represent marrow and white areas represent bone tissue): (1) isolated crack, (2) cross-hatched cracks, (3) diffuse damage and (4) broken trabeculae. The field of view is about 1.5 mm wide.

[ 19 ]A c c e p t e d m a n u s c r i p t

 19 [START_REF] Van Den Abeele | Microdamage diagnostics using nonlinear elastic wave spectroscopy (NEWS)[END_REF] [START_REF] Hamiel | Stable and unstable damage evolution in rocks with implications to fracturing of granite[END_REF] [START_REF] Nazarov | Nonlinear elastic constants of solids with cracks[END_REF].Moreover, very few studies have been undertaken on bone nonlinearity. The classical nonlinear quadratic parameter β (defined as ρ 0 c 0 ∂c ∂p

Fig. 2 .

 2 Fig. 2. Experimental configuration for harmonic generation.

Fig. 3 .

 3 Fig.3. Measurements of the pressure amplitude on the axis of the focused transducer in water. Amplitudes of 2nd and 3rd harmonics are multiplied by 10 for better readability.

Fig. 4 .

 4 Fig. 4. Pressure levels measured in water at 45 mm for different levels of excitation.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Pressure levels measured in water at 45 mm behind the trabecular bone sample for different levels of excitation.

A c c e p t e d m a n u s c r i p tFig. 7 .

 7 Fig. 7. Experimental set-up of the parametric emission method.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Amplitude of the parametric LF wave received by the hydrophone.Fig. 9. Photograph of the human heel bone slice.

A c c e p t e d m a n u s c r i p t ρ 0 2 ∂c 2

 2 ∂p 0

Fig. 10 .

 10 Fig. 10. Experimental configuration for parametric reception.

Fig. 11 .

 11 Fig. 11. Results for phase modulation obtained for water (without sample), Plexiglas, trabecular bone and porous ceramic. The noise level was measured around 1 mrad.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Table 1

 1 Characteristics of samples and experimental results.

						water	Plexiglas	trabecular bone	porous ceramic
	porosity (%)							90	75
	solid phase density (kg/m 3 )		1000	1180	1900	3100
	solid phase velocity (m/s)			1480	2700	3000	5500
	effective density (kg/m 3 )					1090	1525
	effective velocity (m/s)					1650	4000
	β (effective medium)			2.5	14	36	3400
	thickness (mm)					24	24	20
	β (solid phase only)					440	24000
	thickness	of	solid	phase	=			11	13
	thickness×porosity 1/3 (mm)
	sample attenuation, 400 kHz (dB/cm)	0.0003	2	7	18
	sample attenuation, 800 kHz (dB/cm)	0.0013	3	20	36
	sample attenuation, 1.2 MHz (dB/cm)	0.003	5	30	50

A c c e p t e d m a n u s c r i p t