Julianne De

Castro Oliveira

Jean-Baptiste Feret

Jorge Flavio

Yann Ponzoni

Jean-Philippe Nouvellon

Otavio Camargo Gastellu-Etchegorry

José Luiz Campoe

Luiz Stape

Estraviz Carlos

Gueric Rodriguez

Le Maire

Julianne De Castro Oliveira

Jean-Baptiste Féret

Jorge Flávio

Guerric Rodriguez

Simulating the canopy reflectance of different eucalypt genotypes with the DART 3-D model

1



Abstract-Finding suitable models of canopy reflectance in forward simulation mode is a prerequisite for their use in inverse mode to characterize canopy variables of interest, such as Leaf Area Index (LAI) or chlorophyll content. In this study, the accuracy of the 3D reflectance model DART was assessed for canopies of different genotypes of Eucalyptus, having distinct biophysical and biochemical characteristics, to improve the knowledge on how these characteristics are influencing the reflectance signal as measured by passive orbital sensors. The first step was to test the model suitability to simulate reflectance images in the visible and near infrared. We parameterized DART model using extensive measurements from Eucalyptus plantations including 16 contrasted genotypes. Forest inventories were conducted and leaf, bark and forest floor optical properties were measured. Simulation accuracy was evaluated by comparing the mean top of canopy (TOC) bidirectional reflectance of DART with TOC reflectance extracted from a Pleiades very high resolution satellite image. Results showed a good performance of DART with mean reflectance absolute error lower than 2 %. Inter-genotype reflectance variability was correctly simulated, but the model didn't succeed at catching the slight spatial variation for a given genotype, excepted when large gaps appeared due to tree mortality. The second step consisted in a sensitivity analysis to explore which biochemical or biophysical characteristics influenced more the canopy reflectance between genotypes. These results present perspectives for using DART model in inversion mode.

Index Terms-DART, 3D modeling, eucalypt, radiative transfer model, remote sensing

I. INTRODUCTION

MONG the different methods to estimate biophysical or biochemical characteristics of forest plantations, the analysis of the images measured by Julianne de Castro Oliveira and Luiz Carlos Estraviz Rodriguez are with University of São Paulo, ESALQ/USP, Brazil (e-mail: juliannecoliveira@usp.br; lcer@usp.br), Jean-Baptiste Féret is with IRSTEA, UMR TETIS, BP5092 Montpellier, France (e-mail: jeanbaptiste.feret@teledetection.fr), Flávio Jorge Ponzoni is with INPE, Brazil (email: flavio@dsr.inpe.br), Yann Nouvellon is with CIRAD, UMR ECO&SOLS, F-34398 Montpellier, France and with University of São Paulo, ESALQ/USP, Brazil (e-mail: yann.nouvellon@cirad.fr), Jean-Philippe Gastellu-Etchegorry is with CESBIO, France (e-mail : jean-philippe.gastelluetchegorry@cesbio.cnes.fr), Otavio Camargo Campoe is with Federal University of Santa Catarina, UFSC, Brazil (e-mail : otavio.campoe@ufsc.br), José Luiz Stape is with Suzano Pulp and Paper, Brazil (e-mail: stape@suzano.com.br) and Guerric le Maire is with CIRAD, UMR ECO&SOLS, F-34398 Montpellier, France and with NIPE, UNICAMP, Campinas, Brazil (e-mail: guerric.le_maire@cirad.fr).

sensors on orbital platforms is appropriate for large spatial scales studies. Images are converted into reflectance values for each spectral band of the image, and later used to retrieve biophysical parameters of the forest through empirical relationships, or through radiative transfer models (RTM) inversion [START_REF] Le Maire | Calibration of a species-specific spectral vegetation index for leaf area index (LAI) monitoring: example with MODIS reflectance time-series on Eucalyptus plantations[END_REF] - [START_REF] Le Maire | Calibration and validation of hyperspectral indices for the estimation of broadleavedforest leaf chlorophyll content, leaf mass per area, leaf area index and leafcanopy in[END_REF].

RTM explicitly take into account stand structural characteristics (tree dimensions and positions, leaf area index, leaf angle distribution, crown cover, among others) and can simulate the quantitative value of the reflectance spectra of the canopy as observed on top of the canopy or by a sensor onboard a plane or a satellite. They are based on the knowledge of the physical laws that control the transfer and interaction of solar radiation in a vegetative canopy, in interaction with the soil [START_REF] Gastellu-Etchegorry | A modeling approach to assess the robustness of spectrometric predictive equations for canopy chemistry[END_REF]. The DART -Discrete Anisotropic Radiative Transfer -model [START_REF] Gastellu-Etchegorry | A simple anisotropic reflectance model for homogeneous multilayer canopies[END_REF], [START_REF] Gastellu-Etchegorry | Discrete Anisotropic Radiative Transfer (DART 5) for modeling airborne and satellite spectroradiometer and LIDAR acquisitions of natural and urban landscapes[END_REF] is a comprehensive threedimensional model that simulates bidirectional reflectance and enables new possibilities of data analysis to evaluate, for example, canopy structure [START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF], radiative budget [START_REF] Gastellu-Etchegorry | DART: a 3D model for simulating satellite images and studying surface radiation budget[END_REF], [START_REF] Demarez | Modeling of the radiation regime and photosynthesis of a finite canopy using the DART model. Influence of canopy architecture of canopy architecture assumptions and border effects[END_REF], photosynthesis [START_REF] Demarez | Modeling of the radiation regime and photosynthesis of a finite canopy using the DART model. Influence of canopy architecture of canopy architecture assumptions and border effects[END_REF], chlorophyll content [START_REF] Malenovský | Retrieval of spruce leaf chlorophyll content from airborne image data using continuum removal and radiative transfer[END_REF], [START_REF] Demarez | A modeling approach for studying forest chlorophyll content[END_REF], Leaf Area Index (LAI) [START_REF] Banskota | An LUT-Based inversion of DART model to estimate forest LAI from hyperspectral data[END_REF], [START_REF] Banskota | Investigating the utility of wavelet transforms for inverting a 3-D radiative transfer model using hyperspectral data to retrieve forest LAI[END_REF], among others.

Eucalypt plantations in Brazil cover 5.6 million ha, which accounts for 71.9 % of planted forests in Brazil [START_REF]Indústria Brasileira de Árvores[END_REF]. Currently, most areas are planted with several genotypes, mainly on clonal plantations, which have been tested and selected for distinct widespread soils and climatic Brazilian conditions [START_REF] Gonçalves | Integrating genetic and silvicultural strategies to minimize abiotic and biotic constraints in Brazilian eucalypt plantations[END_REF]. These genotypes provide different phenotypes, with distinct canopy structure, leaf morphology and biochemical compounds and biomass production. Due to their high economic importance in Brazil, the understanding of how biophysical parameters of planted forests could explain the spatial-temporal growth dynamics and the estimation of such parameters through remotely-sensed images is of paramount importance [START_REF] Le Maire | Calibration of a species-specific spectral vegetation index for leaf area index (LAI) monitoring: example with MODIS reflectance time-series on Eucalyptus plantations[END_REF], [START_REF] Le Maire | Leaf area index estimation with MODIS reflectance time series and model inversion during full rotations of Eucalyptus plantations[END_REF].

Eucalyptus plantations in Brazil present particular structures: they are planted at high densities (e.g. 1700 trees/ha), they generally have a low leaf area index compared to other dense forests, and they are planted in rows of different spacing (anisotropy). One supplementary difficulty comes from the variability of eucalypts species and genotypes that are planted in Brazil. The different genotypes can have different structural and biophysical properties, even at the same age, and these parameters may change the canopy reflectance in different magnitude. It is therefore necessary to understand better the drivers of the reflectance differences between genotypes to further assess if their estimation through inversion procedures is possible.

2

Despite the successful use of physical approach of DART to retrieve canopies characteristics from inversion procedures, e.g. in [START_REF] Banskota | An LUT-Based inversion of DART model to estimate forest LAI from hyperspectral data[END_REF], [START_REF] Yáñez-Rausell | Estimation of spruce needle-leaf chlorophyll content cased on DART and PARAS canopy reflectance models[END_REF] - [START_REF] Gastellu-Etchegorry | An interpolation procedure for generalizing a look-up table inversion method[END_REF], few detailed studies have tested the efficiency of this 3D reflectance model in forward mode in forest canopy ecosystem [START_REF] Couturier | A modelbased performance test for forest classifiers on remote-sensing imagery[END_REF], [START_REF] Schneider | Simulating imaging spectrometer data: 3D forest modeling based on LiDAR and in situ data[END_REF]. The first assumption of inversion procedure is the suitability of the RTM to simulate accurately the reflectance for a range of canopy characteristics corresponding at least to the range of application conditions. In this study, we parameterized DART model using an extensive in situ measurement dataset. Eucalyptus plantations of 16 different genotypes were used to test the accuracy of the simulations generated by DART when compared with experimental images acquired from a very high spatial resolution satellite, Pleiades. In a second step, we performed a sensitivity analysis using the parameters variability as they were measured in situ to quantify the effect of the main stand parameters (inter-genotype variability) on the canopy reflectance. We finally discussed the use of DART for inversion studies for these particular ecosystems.

II. DATASET DESCRIPTION

A. Study site

The study site is located in Itatinga Municipality, in the state of São Paulo, southeastern Brazil, 22°58'04''S and 48°43'40''W (Fig. 1), as part of the IPEF-Eucflux project. A genotype trial experiment of eucalypt was installed in November 2009 with 16 genotypes comprising several genetic origins from different eucalypt growing companies and regions in Brazil (G1, G2, G10: E. grandis; G3-G9, G11-G13, G15: E. grandis x urophylla; G14: E. saligna; G16: E. camaldulensis x grandis). Fourteen of these 16 genotypes were clones and two (G1 and G2) had seminal origin. Planting rows were mainly east-west oriented, with plant arrangement of 3 m × 2 m (1666 trees per hectare). The experiment comprised 9 blocks, each having 16 treatments (genotypes) randomly distributed within a 4 × 4 subplot grid of 192 trees each (each subplot comprised 12 lines of 16 trees). Only the 10 lines and 10 rows central part of the subplot was analyzed (100 trees, 20 m × 30 m area).

B. In-situ measurements

Forest inventories were carried out at 6, 12, 19, 26, 38, 52, 62 and 74 months of age. During these inventories, trunk diameter at breast height (DBH) and tree height were measured. Close to most of these dates, 10-12 trees were cut for each genotype to compute the biomass per compartment (leaves, branches, trunk and bark) to generate allometric relationships between trunk DBH and tree height, height to the base of the live crown, crown diameter and leaf area, as classically done in other studies in the same area [START_REF] Laclau | Mixed-species plantations of Acacia mangium and Eucalyptus grandis in Brazil -1. Growth dynamics and aboveground net primary production[END_REF] - [START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF]. All these allometric relationships presented good adjustments (e.g. R 2 ~ 0.72, 0.70 and 0.88, respectively, for crown diameter, crown height and leaf area) and included the age as an explanatory variable, allowing their application for each tree at each inventory date. LAI was calculated as the sum of the leaf area of each tree inside the plot divided by the plot area. Leaf angle distribution (LAD) was estimated from the leaf angles measured in the field for each genotype (as described in [START_REF] Le Maire | Leaf area index estimation with MODIS reflectance time series and model inversion during full rotations of Eucalyptus plantations[END_REF]) and adjusted with an ellipsoidal leaf angle density function. In each tree, a clinometer was used to measure the inclination of 72 leaves selected according to their position within the crown to be representative of the tree-scale distribution. The eucalypt stands were analyzed at the date of May, 2014 (54 months), corresponding to the date of satellite image acquisition, using interpolation of the field measurements between inventories at 52 and 62 months. For the leaf area, auxiliary leaf area index values retrieved from more frequent measurements on one of the genotypes allowed to improve the interpolation by considering a common seasonal variation.

Leaves, trunks and forest floor optical properties were measured on October 2015 with an ASD Field SpecPro (Analytical Spectral Devices, Boulder, Colorado, USA) spectrometer in the spectral range from 400 to 2500 nm with 1 nm intervals at 71 months after planting (in October 2015). In these dates, three trees per genotype were selected and for each tree, leaves were collected randomly at three crown layers (bottom, middle and top, divided by exact height proportions) and two horizontal positions in each layer (near and far from trunk), totaling two leaves per crown layer, six leaves per tree and 18 leaves per genotype. These leaves were kept cold and in the dark for less than one hour. Adaxial leaves reflectance and transmittance were measured in the laboratory using an integrating sphere (LI-COR 1800, LI-COR, Inc., Lincoln, Nebraska, USA). Forest floor and bark reflectance were measured using a Contact Probe (ASD, Inc., Boulder, Colorado) on five different points for each genotype, in the same week without rain.

The spectral measurement date occurred more than one year after the satellite image acquisition. However, these component spectra have probably not evolved a lot during this interval: for leaves, there were no significant difference between months 52 and 72 for specific leaf area, water content, and SPAD values (measured with the SPAD-Minolta device) (data not shown). For trunk and forest floor, we assumed no changes, which seem a reasonable hypothesis for these components.

C. Pleiades satellite images

Very high spatial resolution multispectral scenes including four bands (blue: 430-550 nm, green: 490-610 nm, red: 600-720 nm and near infrared: 750-950 nm) from Pleiades satellite were used to validate DART simulations. The image (four bands) was acquired on May 2014, at 13:36 GMT, with the following angles: view azimuth φ 𝑣 = 180.03°, view zenith θ 𝑣 = 13.40°, sun azimuth φ 𝑠 = 33.43° and sun zenith θ 𝑠 = 44.48°. The image was orthorectified and projected. Polygons of each internal plot extension (20 m × 30 m) were used to extract the radiance of the plots in each band of the Pleiades image. Transformation to TOA reflectance was performed, followed by an atmospheric correction to compute the reflectance of the top of canopy (TOC) of the scenes using the 6S model and default atmospheric parameterization for this location [START_REF] Vermote | Second simulation of the satellite signal in the solar spectrum, 6S: An overview[END_REF].

III. ANALYSES AND DART PARAMETERIZATION

A. DART parameterization

DART was used in the ray tracing method and reflectance mode [START_REF] Gastellu-Etchegorry | A simple anisotropic reflectance model for homogeneous multilayer canopies[END_REF], [START_REF] Gastellu-Etchegorry | DART: a 3D model for simulating satellite images and studying surface radiation budget[END_REF] to simulate TOC bidirectional reflectance images. Simulations with DART were conducted on 4 wavebands corresponding to Pleiades sensor relative spectral response.

The input solar angles (θ 𝑠 andφ 𝑠 ) were computed knowing the local latitude, date and hour of satellite overpass. Image acquisition geometry (θ 𝑣 , φ 𝑣 ) was obtained from metadata of Pleiades images. All DART simulated scenes were created using individualized positions and dimensions of the 192 trees of each subplot, but the output stand reflectance computation was restricted to an internal plot of 20 m × 30 m (100 trees), to avoid any border effect. One scene was simulated for each of the 16 genotypes and 9 blocks at 54 months (corresponding to date May, 2014), with computing cubic cells of 0.50 m edge.

Input parameters related to the trees positions (coordinates x and y in the plot), dimensions (e.g. crown diameter and height, DBH and total height), LAI and LAD for each tree were all insitu measurements (described on Section II.B.). For simulating tree crowns, we used a half ellipsoid shape, which typically fit well with the shape of eucalypts crown. Optical properties of the leaves were prescribed in function of the crown layer for each tree (upper, middle and lower) and in function of the genotype, such as the bark and forest floor reflectance. In these canopies, the branches are very thin and represent a very small absorbing surface in comparison to leaves and barks, and therefore they were not simulated.

B. Comparison between simulated and satellite images

The accuracy of the simulated reflectance TOC scenes from DART was checked against the TOC reflectance obtained from Pleiades scenes, for all 4 broadbands (blue, green, red, and NIR), 9 blocks and 16 genotypes. The overall accuracy level for simulating eucalypt plantations was expressed by the mean absolute error (MAE) of each spectral band [START_REF] Willmott | Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance[END_REF]:

𝑀𝐴𝐸 𝜆 = 1 𝑛 ∑ |𝑅 𝑃𝑙é𝑖𝑎𝑑𝑒𝑠,𝜆 (𝑖)-𝑅 𝐷𝐴𝑅𝑇,𝜆 (𝑖)| 𝑛 𝑖=1 , (1)
where𝑅 𝑃𝑙é𝑖𝑎𝑑𝑒𝑠,𝜆 is the reflectance measured by Pleiades satellite for spectral band λ, 𝑅 𝐷𝐴𝑅𝑇,𝜆 is the reflectance simulated by DART for the same spectral band, and 𝑛 is the number of samples (𝑛=144 plots, product of 9 blocks by 16 genotypes). The systematic error (BIAS), root mean square error (RMSE) and the determination coefficient (R 2 ) were also computed, both at genotype scale (averaged by blocks, so n=16 for each band), or for each genotypes for inter-block variability (so n=9 for each band and each genotype).

C. Sensitivity analysis of DART for eucalypt plantations

A simple sensitivity analysis was performed to better understand the effect of inter-genotype differences in structure, biophysical and biochemical parameters on the simulations output. We selected one of the genotype (the G3, that represents the main genotype planted around the experimental area), grown in one of the block (B2, where the plots shows good growth and health) as an example. For each of the parameter listed below, we exchanged one by one the G3 value by the value of another genotype of the same block B2. The range and variation of these values reflected therefore the real inter-genotype variability as it appears on in situ measurements, which enabled more realistic description and analysis of parameters influences on the reflectance. For instance, the LAI of G3 was replaced by the one of G1, the DART reflectance in the four bands were simulated, then a new simulation was performed with the LAI of G2, etc. At the end, we computed the average, variance and produced a boxplot figure for each parameter at each reflectance band. The tested parameters were LAI, LAD, leaf, bark and forest floor optical properties (reflectance), trees dimensions (tree and crown height, crown diameter and DBH), and row azimuth. Note that for the particular case of row azimuth, we changed the orientation by using the orientation of the other blocks one by one, and this is not linked to the genotype. However, including this variability will give more precise information on the importance of this factor. This procedure allows us to better understand which parameters drive the inter-genotypes variability in reflectance.

IV. RESULTS

A. Differences between genotypes structural and biochemical properties

The main characteristics of the genotypes (DBH, height, leaf area, LAI, crown length, crown diameter, leaf angle and mortality) based on field measurements and used for DART parameterization are shown in Fig. 2, together with their interblock repetitions. Overall, we can see that the tree dimensions and structural properties are similar between genotypes having the same age, and high local variability. However, when looking closer, there are some differences between genotypes. The DBH and height values were very similar between genotypes, with higher variability for the seminal material G1 and G2, and higher growth homogeneity of the clonal materials. G16 was the most homogeneous clone. G7, G12 and G16 presented the lowest leaf area values and lowest variability between trees. LAI for all genotypes was around 3-6 m 2 /m 2 and with small spatial variability, mainly for G12 and G16. G10, G11 and G13 presented the highest and G16 the lowest LAI values. In contrast with the tree height, the crown length varied more between genotypes. Similar with tree DBH and height, the crown diameter exhibited little variability across genotypes, with a median around 3 m indicating that at this age (54 months) the trees inside the plots are exploring more or less the space they individually have (3 m × 2 m). Note that there was a small measured difference between within rows and between-rows crown diameters that was included in the simulations. The leaf inclination angle showed high between-trees variability, mostly driven by differences in tree size since there were strong canopy vertical gradients of leaf inclination angles [START_REF] Le Maire | Leaf area index estimation with MODIS reflectance time series and model inversion during full rotations of Eucalyptus plantations[END_REF]. G16 had the highest leaf inclination angles with low variability between trees. Mortality exhibited large variability across genotypes, with the highest values (reaching around 20 %) for genotypes G1, G3, G6 and G13. For other genotypes, the mortality was lower in : IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, n° 11, 2017. p.4844-4852 than 10 %, which are common values for eucalypt plantations [START_REF] Zhou | Mapping local density of young Eucalyptus plantations by individual tree detection in high spatial resolution satellite images[END_REF].

Leaves, trunks and forest floor optical properties are shown on Fig. 3 for each genotype. Leaf reflectance (shown on Fig. 3 for expanded mature leaves of the middle crown layer) exhibited high absorption peaks in the blue and red regions and high NIR reflectance for all genotypes. Note that the reflectance ranking between genotypes was conserved for all wavelengths in the visible but changed further in the NIR and MID regions. There were larges differences of bark reflectance between genotypes. Interestingly the reflectance was very high in the visible and NIR regions compared to leaf reflectance. Some spectra clearly show an absorption feature in the red region. Forest floor reflectance showed similar pattern for all genotypes, but with a high inter-genotype variability, with low reflectance in the visible region and increasing values along the spectrum, and a mild absorption peak in the water absorption band (1400 nm).

Fig. 4 shows the leaves reflectance in the green, red and near infrared bands for each crown level (bottom, middle and top) and genotypes. There was no significant difference between crowns layers, significant differences between genotypes and no significant differences for the interaction genotype × crown layers for each band (N-Way ANOVA under Matlab 2013a, α=0.05). Note that statistical analysis was done using all measured reflectance data instead of the average values shown on Fig. 4.

B. Comparison of DART simulations with Pleiades satellite image

The TOC reflectances simulated by DART and acquired by the Pleiades sensor at the four multispectral bands for each genotype are shown on Fig. 5, averaged by genotype and with standard deviation. In general, the mean TOC reflectances from DART simulations were in good agreement with the mean TOC reflectance of the Pleiades scenes for all four bands and genotypes. Discrepancies were found mainly for the blue band (430-550 nm) for all genotypes, and some discrepancies appeared in the near infrared band (750-950 nm) for some genotypes (e.g., genotypes 5, 8 and 12). A numerical comparison between the reflectance simulated by DART and acquired by Pleiades scenes was performed using the MAE, RMSE and R 2 for all blocks and genotypes in each band (Table 1). The minimum and maximum range of R 2 values computed for each genotype for the four bands are also presented.

The MAE values were low for all bands (< 0.0195), with the lowest values for the green band. Higher values were found for blue and NIR bands, which corroborates the results of Fig. 5. BIAS, that represents the average difference between Pleiades and DART reflectance, was negative and indicated that TOC reflectances simulated by DART model were, in general, slightly higher than TOC reflectances derived from Pleiades images. RMSE were also low (<0.023), mainly for the bands in the visible domain (<0.0023). NIR band had the higher value. The R 2 best performance was for red and NIR bands and worst for blue band. The R 2 for each genotype computed with the different blocks (spatial variability) in all bands showed a wide range of values. The spatial variability of some genotypes was correctly simulated whereas others were not significant. An example of the level of detail of trees parameterization on DART simulated scenes compared with Pleiades scenes is shown on Fig. 6 for the near infrared band of G14, block 4. (scenes with 0.50 m and 2 m of spatial resolution). The near infrared (2.0 m of spatial resolution) and panchromatic band of Pleiades (0.50 m of spatial resolution) scenes for the same G14 and block 4 are also presented. The Pleiades panchromatic was chosen to present this example, due to the higher spatial resolution of this band. This visual comparison illustrates how the DART model represents the canopy. We can see that DART simulations are in accordance with the image in terms of shadow proportion, gaps, row orientations, textures and object dimensions. However, the model cannot reach the level of detail for a use on a tree-by-tree analysis in this type of canopy structure.

C. Sensitivity Analysis

The results of the sensitivity analysis of the simulated reflectance for the blue, green, red and NIR bands according to stand parameters (LAI, LAD, leaf, bark and forest floor optical properties -reflectance, trees dimensions and row azimuth) are presented in Fig. 7. The behavior of real range variation of each parameter individually (without interaction between them) on the average canopy reflectance was presented together to compare their magnitude. LAI, leaf reflectance, trees dimensions and row azimuth had the highest sensitivity and explain most of the difference between genotypes in the visible bands. These variability were of the same order of magnitude as the variability due to row orientation. Bark and forest floor reflectance and LAD showed the weakest sensitivity in these bands despite their intergenotype variability being relatively high. The NIR band showed similar reflectance results among the replacing tests, but with higher inter-genotype standard deviations compared to the others bands. The LAD, bark and forest floor reflectance showed higher influence in the NIR band compared to visible bands.

V. DISCUSSION

A. Parameterization of DART

Overall, the differences between eucalypt trees of different genotypes and locations were not very large for many of the in : IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, n° 11, 2017. p.4844-4852 parameters. However, the final importance of a parameter to explain the difference in TOC reflectance between genotypes (and/or locations) is a conjunction of the inter-genotype variability (or spatial variability) of this parameter, and the sensitivity of that parameter in the model. It is therefore important, before setting some of the DART parameters to constants (and therefore not explaining the genotype or spatial variability), to model the system with the maximum precision, and simplify afterwards if possible. The model parameterization is therefore a critical step of this work.

The leaf reflectance was shown to be different between genotypes, reflecting differences in pigment contents, and internal structure of leaves. A more detailed analysis could be done to assess which leaf structural or biochemical characteristics could explain this reflectance variability, but such analysis is out of the scope of the study: here we focused our analysis on the macro-scale differences between genotypes, and leaf reflectance was therefore an input parameter of DART.

The high inter-genotypes difference of bark reflectance (Fig. 3) was expected, since their color and roughness was extremely different in the field. The absorption feature in the red is associated to the presence of chlorophyll pigments in the bark surface for some of the genotypes, as observed in many other studies (e.g. in [START_REF] Wittmann | The optical, absorptive and chlorophyll fluorescence properties of young stems of five woody species[END_REF], [START_REF] Girma | Photosynthetic bark: use of chlorophyll absorption continuum index to estimate Boswellia papyrifera bark chlorophyll content[END_REF]). There was also a high intergenotype variability on forest floor reflectance (Fig. 3), mainly in the NIR and MID regions. This behavior is due to the different composition of the forest floor materials (e.g. green or yellowing leaves just fallen and dead dry leaves, bark and branches proportion, leaf sizes), their structural variance, moisture content and decomposition stage [START_REF] Asner | Variability in leaf and litter optical properties: implications for BRDF model inversions using AVHRR, MODIS, and MISR[END_REF]; which directly influences the reflectance.

The ANOVA analysis of the leaf reflectance for bottom, middle and top crown layers in the green, red and NIR bands (Fig. 4) showed that there was no statistical significant difference between bottom, middle and top crown layers but there were differences between genotypes considering all crown layers. Therefore, the use of different spectra for upper, middle and lower part of the canopy could be unnecessary for simulating reflectance in these wavebands. However, since some genotypes showed different spectra for upper layer, which could be locally important for TOC simulation, we preferred to keep this detailed description in the simulations. Also, the leaves inside each crown layers have different combinations of development stages (juveniles and mature). Generally there is a gradient of these development stages inside crown, with more juvenile leaves in the top layer and more mature leaves in the bottom. Mature leaves have more pigments, higher mass per area than juvenile [START_REF] Stone | Spectral reflectance characteristics of eucalypt foliage damaged by insects[END_REF], and different internal structure, which directly influence the reflectance in the visible and NIR wavelengths. However, our results did not clearly show any vertical trend of reflectance between crown layers. The explanation is that the proportion of juvenile leaves in the top layer is variable between genotypes, and between trees of different heights.

B. Suitability of DART for TOC simulations

Assessing if a RTM is suitable to simulate a given ecosystem depends on the objective of the study. In this study, we can distinguish the results in function of the level of variability of the observed canopy, i.e., evaluating the degree of precision of DART for simulating i) a "typical" Eucalyptus plantation reflectance, ii) inter-genotype reflectance variability and iii) the inter-block reflectance variability for the same genotype.

Our results showed that the DART model was suitable to simulate Eucalyptus plantation in general, with their very high tree density, their tall trunks, bright forest floor, and ellipsoidal form of their crown (Fig. 5): this is especially underlined by the low MAE obtained for this ecosystem (lower than 2%). The inter-genotype variability of reflectance comes from the variability of many structural and biochemical parameters of the ecosystem, as represented in Fig. 2 and Fig. 3 (e.g. optical properties of the different components, leaf angles, dimensions of trees, etc.). This inter-genotype variability was adequately simulated as could be seen on Fig. 5 and Table 1, with coefficients of determination > 0.41 for all spectral bands, and of 0.55 in the NIR bands. Such a bias for blue band could come from residual atmospheric effects not properly taken into account in the atmospheric correction of the Pleiades images, which was based on standard atmospheric parameterization of 6S in absence of local measurements of atmospheric water, ozone and aerosols contents.

Finally, the simulations of the spatial variability between blocks, for each genotype were not adequately simulated for most genotypes. There were very low average coefficients of determination in all bands considering each genotype for all blocks. The spatial variability for a given genotype is more difficult to assess by simulation, mainly because of the low variability existing between these blocks. Therefore, the precision of the simulation is not sufficient to catch up this spatial variability. However, some genotypes had higher mortality rates (e.g. G1, G3, G6 and G13), which created large gaps in the canopy and increased the variability to a range possible to simulate (high R 2 scores). As a consequence, the use of DART model in inversions mode for these ecosystems would gain precision if the genotype is already known, and in areas where the proportion of gaps remains low. Moreover, the row orientation could also act as a confounding factor and should also be prescribed prior to inversion, i.e. a pre-analysis of row orientations needs to be carried on.

In terms of bi-directional TOC reflectance, the comparison between simulated and real satellite scenes from forest stands is a difficult task, since the reflectance image is dominated by the macroscopic properties of the illuminated and shadowed crowns as well as ground surface [START_REF] Houborg | Combining vegetation index and model inversion methods for the extraction of key vegetation biophysical parameters using Terra and Aqua MODIS reflectance data[END_REF], as illustrated in Fig. 6, at very high resolution.

Our results confirm the ability of DART to simulate remote sensing data under several eucalypt forest conditions. Some comparisons between DART simulations and forest ecosystems reflectance was also done in [START_REF] Couturier | A modelbased performance test for forest classifiers on remote-sensing imagery[END_REF], [START_REF] Schneider | Simulating imaging spectrometer data: 3D forest modeling based on LiDAR and in situ data[END_REF] and the main conclusions were that DART showed very low pixels spectral dissimilarity compared with IKONOS images and R 2 of 0.48 for a pixel-wise comparison with APEX imaging spectrometer, respectively. DART has been successfully compared with other 3D models throughout the RAdiative transfer Models Intercomparison -RAMI exercise [START_REF] Pinty | RAdiation transfer Model Intercomparison (RAMI) exercise: results from the second phase[END_REF], [START_REF] Widlowski | The fourth phase of the radiative transfer model intercomparison (RAMI) exercise: actual canopy scenarios and conformity testing[END_REF] under several conditions. Our results extend DART model validation on real measured dataset of individualized trees and stands of Eucalyptus plantations, which have particular characteristics (e.g. a high tree density but rather low LAI, lots of trunk surface but few branches).

C. Source of the inter-genotype reflectance variability

After having tested the model suitability to simulate intergenotype TOC reflectance variability, we seek to address which of the stand structural or biochemical parameters (LAI, LAD, leaf, bark and forest floor optical properties, trees dimensions and row orientation) influences more the reflectance between genotype (Fig. 7). These parameters were chosen since they are the main input parameters of DART.

The LAI was one of the most influencing parameter for explaining the difference of reflectance between genotypes. Numerous studies have proved that vegetation reflectance is strongly affected by LAI in the entire spectra, but more in the NIR [START_REF] Shi | Consistent estimation of multiple parameters from MODIS top ofatmosphere reflectance data using a coupled soil-canopy-atmosphereradiative transfer model[END_REF] - [START_REF] Le Maire | Calibration and validation of hyperspectral indices for the estimation of biochemical and biophysical parameters of broadleaves forest canopies[END_REF]. The leaf reflectance, which reflects in the visible the different leaves pigments contents, was another very important factor driving the canopy reflectance, mainly in the visible region. These results agree with [START_REF] Xiao | Sensitivity analysis of vegetation reflectance to biochemical and biophysical variables at leaf,canopy, and regional scales[END_REF], which performed a sensitivity analysis of vegetation reflectance and found more influence of leaves pigments content in the visible and LAI in the NIR regions at canopy scale. They also showed a weak effect of leaf angle at this scale.

The crown dimensions also explained the difference of TOC reflectance between genotypes (Fig. 7), as shown in other studies [START_REF] Rautiainen | The effect of crown shape on the reflectance of coniferous stands[END_REF]. This variable, jointly with the row azimuth, mainly drives the proportion of visible soil between rows and the proportion of shaded/illuminated crowns on the image. The presence of empty spaces (dead trees) in some of the plots increased even more this heterogeneity, which also increased the contribution of this parameter to the inter-genotype and spatial variability of TOC reflectance. Some of the parameters tested here showed moderate sensitivity on simulated TOC reflectance, which is the case for bark and forest floor reflectance. Therefore, average values could have been chosen for these parameters, and could simplify further inversions. In contrast, TOC reflectance showed high sensitivity to LAI, leaf reflectance, trees dimensions and row azimuth. It seems therefore important to perform genotype-specific inversion in the future, or grouping genotypes for their crown dimensions. Also, knowledge of the row orientation will be critical for inversion purposes. Further step will be to simulate a comprehensive database along eucalypts growth stages for different genotypes, and use this database to estimate some variables such as the LAI or chlorophyll content through inversion procedures. Our first sensitivity analysis can further help distinguish the inversion errors coming from the model itself or coming from the inversion methodology (algorithm, constraints, etc.).

These sensitivity analysis results confirm the relevance of using 3D models such as DART, as they are particularly suitable to explicit the influence of tree shape, leaf pigments and plot heterogeneity on the canopy reflectance of different genotypes and row orientations.

VI. CONCLUSION

In this study we tested DART model to simulate Eucalyptus plantation reflectance, their difference between genotypes and between plots for a given genotype. DART was reliable for eucalypts plantation simulation in general, and adequately simulated the difference of reflectance between 16 genotypes including the mostly planted ones in the region, and some particular genotypes (e.g. G16: E camaldulensis x grandis). However, the local difference of reflectance was correctly simulated only when the range of TOC reflectance was high for a given genotype, which occurred mainly through local mortality.

The difference of TOC reflectance in the visible bands between genotypes is mainly explained by differences in LAI, leaf optical properties and row orientation. In the NIR, the same parameters influence the TOC canopy, together with the tree dimensions. Leaf angles, bark and forest floor reflectance have a smaller effect in comparison to the other parameters, although their inter-genotype variability was large.

Successful test of DART in forward mode for simulating the TOC reflectance of these different genotypes open possibilities for parameter estimation through model inversion procedures for eucalypt plantations.

FIGURESFig. 1 .

 1 FIGURES

Fig. 2 .

 2 Fig. 2. Main stand structural characteristics (diameter at breast height -DBH, tree height, tree leaf area, leaf area index -LAI, crown length, crown diameter, leaf inclination angle and mortality) of the 16 genotypes on May, 2014. Mortality represents the percent of dead trees in each block per genotype. Lines inside boxes are the median values, inferior and superior boxes limits are the first and third quartiles, respectively; and error bars outside boxes extend from minimum and maximum values within 3 standard deviations. Variability considered here is the tree-scale variability considering all blocks. Mortality and LAI variability is inter-block variability.

Fig. 4 .

 4 Fig. 4. Leaves reflectance in the green, red and near infrared regions at bottom, middle and top crown layer for the 16 genotypes (labeled as G1 to G16).

Fig. 5 .

 5 Fig.[START_REF] Gastellu-Etchegorry | A modeling approach to assess the robustness of spectrometric predictive equations for canopy chemistry[END_REF]. DART (light gray) and Pleiades (dark gray) mean top of canopy (TOC) reflectance of four bands (B=blue, G=green, R=red, NIR=near infrared) for each genotype averaged for all blocks and subplots. Lines in each bar represent the standard deviation for blocks.

Fig. 6 .

 6 Fig. 6. Example of near infrared DART simulated scene with 0.50 m (a) and 2 m (b) of spatial resolution, panchromatic Pleiades image(c) with 0.50 m and near infrared Pleiades image with 2 m of spatial resolution for the genotype 14 in the block 4.

Fig. 7 .

 7 Fig. 7. Sensitivity analysis of the reflectance in blue, green, red and near infrared bands relative to stand parameters (respectively, LAI, LAD, leaf, bark and forest floor reflectance, trees dimensions and row azimuth). Boxplot definition is given in Fig. 2. Dashed green line represents the TOC reflectance of the genotype 3 (reference). Numbers above each boxplot are the standard deviation. Red crosses are the outliers values.Parameters

TABLE 1 MEAN

 1 ABSOLUTE ERROR (MAE), SYSTEMATIC ERROR (BIAS), ROOT MEAN SQUARE ERROR (RMSE) AND DETERMINATION COEFFICIENT (R 2 ) FOR SIMULATED BANDS (BLUE, GREEN, RED AND NIR) IN RELATION TO PLEIADES BANDS, AVERAGED BY GENOTYPE AND BLOCK. R 2 OF GENOTYPES (MIN. -MEAN -MAX.) IS THE MINIMUM, MEAN AND MAXIMUM R 2 VALUES IN EACH BAND FOR THE GENOTYPES, COMPUTED ON THE INTER-BLOCK VARIABILITY.

	Band	MAE	BIAS	RMSE	R 2	R 2 of genotypes (min. -mean-max.)
	Blue	0.0180 -0.0180 0.00106 0.41	0.0003-0.11 -0.79
	Green 0.0063 -0.0063 0.00223 0.43	0.0003-0.12 -0.88
	Red	0.0170 -0.0170 0.00104 0.51	0.0003-0.12 -0.75
	NIR	0.0194 -0.0044 0.02200 0.55	0.0023 -0.28-0.91

in : IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, n°11, 2017. p.4844-4852

ACKNOWLEDGMENT This study was funded by Forestry Science and Research Institute -IPEF (EUCFLUX project -funded by Arcelor Mittal, Cenibra, Copener, Duratex, Fibria, International Paper, Klabin, Suzano and Vallourec), the HyperTropik project (TOSCA program grant of the French Space Agency, CNES), CNPq, CAPES and the Centre de Coopération Internationale en Recherche Agronomique pour le Développement -CIRAD. Pleiades image was acquired in the frame of GEOSUD program, a project (ANR-10-EQPX-20) of the program "Investissementsd'avenir" managed by the French National Research Agency. We are grateful to Eloi Grau (IRSTEA) for support in running the program, José Guilherme Fronza for field work, and the staff at the Itatinga Experimental Station, in particular Rildo Moreira e Moreira (ESALQ, USP) and Eder Araujo da Silva (http://www.floragroapoio.com.br) for their technical support.


H Gómez
email: hgomez@apc.univ-paris7.fr

D Gibert

C Goy

K Jourde

Y Karyotakis

S Katsanevas

J Marteau

M Rosas-Carbajal

A Tonazzo

Forward scattering effects on muon imaging

Keywords:

A

: Muon imaging is one of the most promising non-invasive techniques for density structure scanning, specially for large objects reaching the kilometre scale. It has already interesting applications in different fields like geophysics or nuclear safety and has been proposed for some others like engineering or archaeology. One of the approaches of this technique is based on the well-known radiography principle, by reconstructing the incident direction of the detected muons after crossing the studied objects. In this case, muons detected after a previous forward scattering on the object surface represent an irreducible background noise, leading to a bias on the measurement and consequently on the reconstruction of the object mean density. Therefore, a prior characterization of this effect represents valuable information to conveniently correct the obtained results. Although the muon scattering process has been already theoretically described, a general study of this process has been carried out based on Monte Carlo simulations, resulting in a versatile tool to evaluate this effect for different object geometries and compositions. As an example, these simulations have been used to evaluate the impact of forward scattered muons on two different applications of muon imaging: archaeology and volcanology, revealing a significant impact on the latter case. The general way in which all the tools used have been developed can allow to make equivalent studies in the future for other muon imaging applications following the same procedure.

Introduction

The idea to use muons produced in the Earth's atmosphere by cosmic rays as a scanning method of anthropic or geological structures, the so-called muon imaging, was proposed soon after the discovery of these muons [START_REF] Neddermeyer | Note on the nature of cosmic-ray particles[END_REF][START_REF] Neddermeyer | Cosmic-ray particles of inter-mediate mass[END_REF][START_REF] Auger | Les rayons cosmiques[END_REF]. Muon imaging leverages the capability of cosmic muons to pass through hundreds of metres or even kilometres of ordinary matter with an attenuation mainly related to the length and density of this matter encountered by the muons along their trajectory before their detection [START_REF] Nagamine | Introductory Muon Science[END_REF]. As this attenuation is principally caused by muons absorption and scattering, muon imaging can be mainly performed using two different techniques. The first one is the so-called transmission and absorption muography [START_REF] Lesparre | Geophysical muon imaging: feasibility and limits[END_REF][START_REF] Marteau | Muons tomography applied to geosciences[END_REF]. This technique relies on the well-known radiography concept (widely used, for example, in medicine with X-rays), based on the muon energy loss and its consequently probability to cross a given amount of material. The second is known as deviation muography, which relies on the measurement of the muon track deviation to determine the object density [START_REF] Borozdin | Radiographic imaging with cosmic-ray muons[END_REF][START_REF] Procureur | Muon imaging: Principles, technologies and applications[END_REF]. For the first technique, studying all the directions for which muons go through the studied object, and knowing its external shape, it is possible to obtain a 2D mean density image. Thus, muon imaging provides a non-invasive and remote scanning technique utilisable even for large objects, where the detection set-up may be relatively far away from the -potentially dangeroustarget (e.g. domes of active volcanoes or damaged nuclear reactors).

One of the first studies performed based on muon imaging dates from 1955, being the scanning of the rock overburden over a tunnel in Australia [START_REF] George | Cosmic rays measure overburden of tunnel[END_REF]. Later, other applications from mining [START_REF] Malmqvist | Theoretical studies of in-situ rock density determination using cosmic-ray muon intensity measurements with application in mining geophysics[END_REF] to archaeology were proposed, both in the 70s. For the latter case, some measurements have been already performed, as the exploration of the Egyptian Chephren [START_REF] Alvarez | Search for hidden chambers in the pyramids using cosmic rays[END_REF] and Khufu [START_REF] Morishima | Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons[END_REF] pyramids. Nowadays, thanks to the improvements on the detector performance, and also to their autonomy and portability, muon imaging reveals itself to be a scanning technique competitive and complementary to others non-invasive methods as seismic and electrical resistivity tomography or gravimetry. This has led to its proposal and utilisation in a wide range of fields.

In addition to the above-mentioned applications (archaeology and mining), two others stand out. The first one is related with geophysics, more precisely with the monitoring of volcanoes. This A c c e p t e d m a n u s c r i p t has an important benefit both from a scientific and social point of view. The continuous monitoring of volcanoes helps to understand their internal dynamics, a key feature in the risk assessment. The other application, more related with particle physics, was motivated by the necessity to characterize the overburden of underground laboratories hosting various experiment detectors. It is worth mentioning other applications related to civil engineering and nuclear safety. For the first one, it will be possible for example to scan structures looking for defects. For the case of nuclear safety, set-ups looking for the transport of radioactive materials and wastes already work in cooperation with homeland security agencies. Moreover, the study of nuclear reactors looking for structural damages have been already used, as in the case of the recent Fukushima nuclear power plant accident [START_REF] Kume | Muon trackers for imaging a nuclear reactor[END_REF], and it is being considered as a remote scanning method.

As mentioned, the improvement on the detectors used for muon imaging has been one of the main reasons for the renewal of this technique. Better detectors provide a better angular resolution for the muon direction reconstruction and improve the precision of the density radiography. Nonetheless, the background muon flux rejection remains a key procedure for the structural imaging with muons. An important potential noise source, specially in the measurements based on the transmission and absorption muography, is the forward scattering of low energy muons on the object surface, reaching afterwards the detector. This effect mimics through-going particles since the reconstructed direction of the scattered particle points towards the target. The result is an increase of the total number of detected particles, as if the target's opacity was lower than its actual value, leading to a systematic underestimation of the density [START_REF] Nishiyama | Monte Carlo simulation for background study of geophysical inspection with cosmic-ray muons[END_REF][START_REF] Rosas-Carbajal | Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data[END_REF]. Being produced by muons, these events can not be rejected by particle identification techniques, representing an irreducible background. For this reason, an evaluation of the magnitude of this effect is mandatory to conveniently correct the reconstructed object density.

In this work, a general evaluation of forward scattering of muons has been performed by Monte Carlo simulations. The aim was to develop a versatile tool to be able to evaluate this process for different object geometries and compositions, due to the increasing number of proposed applications based on muon tomography. The main features and results are presented in section 2. Then the impact of this process on the muon imaging capability has been evaluated defining a signal to background parameter. Two physics cases have been studied in section 3. The first one concerns an archaeological target, the Apollonia tumulus near Thessaloniki in Greece, and the second one La Soufrière volcano in Guadeloupe Islands of the Lesser Antilles. Finally, a summary of the different results and the main conclusions extracted from them are compiled in section 4.

Evaluation of the forward scattering of muons

As mentioned in the introduction, low-energy cosmic muons can change their original direction after interacting with the target or any other object in the surroundings before their detection. As muon imaging is based on the reconstruction of the detected muons direction, these muons would forge the measurement. As a consequence, the determination of the target's internal structure and the corresponding reconstructed mean density will be affected.

Muons trajectory deviation is mainly driven by their interaction with matter via multiple Coulomb scattering. The resulting deflection angular distribution, theoretically described by the Molière theory [START_REF] Bethe | Molière's Theory of Multiple Scattering[END_REF], roughly follows a Gaussian, -2 -A c c e p t e d m a n u s c r i p t

dN dα = 1 √ 2πα M S e -α 2 2α 2 M S (2.1)
which is centred in zero (i.e. no deflection happens), having a standard deviation α M S :

α M S = 13.6MeV βcp Q x X 0 (1 + 0.038ln(x/X 0 )) (2.2)
where β is the relativistic factor, p the muon momentum in MeV/c, x the material thickness and Q the absolute electric charge of the muon. α M S also depends on the radiation length (X 0 ) which is empirically given by

X 0 ≈ 716.4g/cm 2 ρ A Z (Z + 1)log(287/ √ (Z )) (2.3)
with Z and A the atomic and mass numbers respectively and ρ the material density. This reveals the relationship of the multiple Coulomb scattering with the properties of the studied material. Different works (see for example [START_REF] Schenider | Coulomb scattering and spatial resolution in proton radiography[END_REF]) provide analytical solutions to the angular distribution of deflected muons after traversing an object with a determined geometry and composition. Besides, other relevant features, as the higher scattering probability for lower energy muons, are also demonstrated in these studies. However, the increasing number of different applications proposed for muon tomography, implies a large variety of objects dimensions, shapes and compositions, being less evident to obtain an analytical estimation of the forward scattering process suitable for all these cases. In this context, Monte Carlo simulations represent a useful tool for the study of muon scattering process, versatile enough to adapt them to the main features of each particular case. As first step on the development of these simulations, a general evaluation of the muons forward scattering has been carried using the Geant4 simulation tool-kit [START_REF] Agostinelli | GEANT4: A Simulation toolkit[END_REF]. It allows the simulation of the 3D muon transport through the defined geometry taking into account the energy loss and trajectory variations due to multiple Coulomb scattering as well as to ionization, bremsstrahlung, pair production and multiple inelastic scattering. Considering these possible processes, results can be compared with the estimations given by the analytical formulas above-mentioned. A scheme of the simulated set-up is shown in figure 1. For this case, generated muons are thrown to a fixed point on a standard rock surface (with a density of 2.5 g/cm 3 ). In the case of scattered muons, the direction changes, in zenith and/or azimuth angles, can be evaluated.

A first set of simulations were performed in order to evaluate the general features of the muon forward scattering. In the previously described set-up, muons up to 10 GeV, with a zenith incident angle (θ ini det ) between 70 • and 90 • and an azimuth incident angle ϕ ini det = 0 • were generated. It is worth mentioning that by the set-up definition θ ini det = 0 • implies muons perpendicular to the rock surface, while θ ini det = 90 • corresponds to tangential ones. Figure 2 summarizes the results of this general simulation, leading to some conclusions about the muon forward scattering studied in these simulations. First, it is observed that this process is negligible if the muon energy is higher than 5 GeV, independently of the incident direction. For the lower energy muons, most of the "efficient" scattering processes (i.e. when the scattered particle exits the medium) occur if θ ini det is higher than 85 • and do not exist if θ ini det is lower than 80 • . That means that only low energy muons with incident directions close to the surface tangent are likely to be scattered on the object surface and to induce a -3 - signal in the detector. For these muons the angular deviation can reach up to 25 • both for the zenith and azimuth angles. By the simulation set-up definition, only the azimuth scattering angle (∆ϕ det ) has been registered for the whole angular range. As presented in figure 2, the ∆ϕ det distribution for all the muon energies considered is in agreement with the Gaussian predicted by Molière theory, as well as the other extracted conclusions agree with the analytical predictions [START_REF] Patrignani | Review on Particle Physics[END_REF].

Taking into account this general information and having checked the agreement between the general simulations and the analytical predictions, a more detailed simulation, optimizing the initial muon sampling was performed. The objective was to establish a probability density function (PDF) to further estimate the background due to forward scattered muons that could be detected during a muon imaging measurement and should be considered in the image analysis. With this aim 10 8 muons homogeneously distributed up to 5 GeV and with θ ini det between 85 • and 90 • (all with ϕ ini det = 0 • ) were generated and simulated in the described Geant4 framework. The generated PDF provides a probability value P(θ ini det , θ f in det , E ini µ ) depending on the initial and final zenith angle and the initial muon energy. A summary plot of the generated PDF divided in 0.5 GeV energy windows is shown in figure 3.

At this point it is worth mentioning that for the studies presented in this work (summarized in section 3), the considered composition of the studied objects are the standard rock used to generate the PDF, but also a definition of soil with different composition and density than the rock (ρ = 2.2 g/cm 3 ). Moreover, there exist several types of rocks and soils with different compositions and densities typically, between 2.0 and 2.5 g/cm 3 . For this reason the influence of these two parameters in the PDF generation has been evaluated: a set of dedicated simulations have been performed changing the composition and the density of the target to compare their results. The obtained PDFs, including the standard soil case, agree to better than 97 %. Thus, the PDF presented in figure 3

Signal to background ratio estimations

The impact of the forward scattered muons in an imaging measurement for a particular object can be evaluated based on simulations as those presented in section 2. This impact can be expressed as a signal to background ratio (S/B) for a given direction θ z -ϕ z . These spherical coordinates correspond to those centred at the detector where θ z = 0 • is the vertical direction and ϕ z = 0 • points to the main axis of the studied object. The signal S(θ z , ϕ z ) is estimated as not scattered muon flux, so their reconstructed direction corresponds to their initial one. The background B(θ z , ϕ z ) represents the scattered muons for which the reconstructed direction points towards the target.

As mentioned, these evaluations allow the study of a particular object, with its corresponding composition. For this it is necessary to know its external shape, to assume the object mean density (since this is the observable that can be extracted from a muon tomography measurement), and to determine the muon detector position with respect to this object. This allows the estimation of the -5 - GeV windows between 0 and 5 GeV for the initial simulated muon energy (E ini µ ). These plots, correspond to 10 8 simulated muons with incident angles between 85 • and 90 • and energies between 0 and 5 GeV (both of them homogeneously distributed). They are used as PDF for further estimations on the forward scattered muon flux.

-6 -A c c e p t e d m a n u s c r i p t object length traversed by muons for each direction as well as its surfaces positions with respect to the detector and to the Earth's surface (required for the determination of θ z and ϕ z ).

For this work two cases have been considered, corresponding to two applications of the muon imaging: the archaeology and the volcanology. For the first one a Macedonian tumulus located near Apollonia (Greece) has been studied [START_REF] Gómez | Studies on muon tomography for archaeological internal structures scanning[END_REF]. For the second, La Soufrière volcano (Guadeloupe island in the Lesser Antilles), already explored by muon imaging, has been taken as reference [START_REF] Jourde | Experimental detection of upward going cosmic particles and consequences for correction of density radiography of volcanoes[END_REF][START_REF] Jourde | Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano[END_REF].

Archaeology: Apollonia tumulus

As quoted in section 1, the exploration of archaeological structures is one of the applications for which muon imaging has been proposed since it is non invasive and does not induce any harmful signals (contrary, for example, to vibrations used in seismic tomography). Already suggested in the 60s [START_REF] Alvarez | Search for hidden chambers in the pyramids using cosmic rays[END_REF], there exist at present different projects based on muon imaging devoted to the study of the internal structure of archaeological constructions (see for example [START_REF] Morishima | Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons[END_REF][START_REF]ScanPyramids project[END_REF]). The ARCHé project proposes to scan the Apollonia Macedonian tumulus [START_REF] Gómez | Muon imaging for archaeological applications: feasibility studies of the Apollonia Macedonian Tumulus[END_REF]. These tumuli are man-made burial structures where the tomb, placed on the ground, is covered by soil, creating a mound which can also contain internal corridors. The geometry and dimensions of these tumuli are variable but they can always be approximated to a truncated cone. In the case of Apollonia tumulus its height is 17 m while the radius of the base and the top are 46 m and 16 m respectively. With this geometry the angle of the slope of the lateral surface of the tumulus is 29.5 • .

In the present study a standard soil composition, with a density of 2.2 g/cm 3 , has been assumed. The detector has been placed 4 m beside the tumulus base (50 m from the tumulus base centre), so muons with zenith angles θ z > 63.4 • are those which will provide information about the structure of the tumulus. With all these features the signal and background, S(θ z , ϕ z ) and B(θ z , ϕ z ) respectively, can be estimated for a given direction θ z -ϕ z . As already described, these coordinates are centred at the detector and θ z = 0 • correspond to vertical muons, while ϕ z = 0 • points towards the centre of the tumulus base.

From the knowledge of the external shape, it is possible to determine the length of tumulus traversed by muons for a given direction L(θ z , ϕ z ) and thus, the corresponding opacity as the product of this length by the density ( = L ×ρ). The required minimal muon energy (E min ) to cross the target of opacity can be calculated as

E min = a b × e b× -1 (3.1)
where a(E) and b(E) represent the energy loss coefficients due to ionization and radiative losses respectively. In this case, coefficients corresponding to standard rock summarized in [START_REF] Patrignani | Review on Particle Physics[END_REF] have been used, obtaining E min values as a function of . As a cross-check, these E min values have been also estimated from the CSDA range values of standard rock [START_REF] Groom | Muon stopping-power and range tables[END_REF]. The agreement between both methods is better than 95 %. Hence, the expected signal S(θ z , ϕ z ) corresponds to the muon flux on the studied direction with energies higher than E min :

S(θ z , ϕ z ) = ∞ E=E mi n φ µ (E, θ z , ϕ z )dE (3.2)

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

To compute the background due to muons forward scattered in the same direction, B(θ z , ϕ z ), two assumptions have been done. First, a point-like detector is considered. This implies that for each scattering point on the tumulus surfaces there is a unique final direction reaching the detector. Second, scattering effects in the azimuth angle are neglected. Since the general muon scattering studies (section 2) show that these effects are symmetric and mostly below 5 • for the azimuth angle (see figure 2), a low influence on the overall estimation is expected, having fade-out effects among the different azimuth directions.

With these two assumptions, B(θ z , ϕ z ) corresponds to the product between the initial flux of muons which can be scattered by the corresponding probability to be scattered with a final zenith angle θ z . As already shown, only muons up to 5 GeV with an incident zenith angle higher than 85 • with respect to the surface normal need to be considered for the forward scattering studies. This delimits the energy and zenith angle ranges to estimate the initial muon flux. The scattering PDF,

P(θ ini det , θ f in det , E ini µ )
, corresponds to the one presented in figure 3. This PDF was generated using the coordinates θ det -ϕ det , centred in the scattering point and orthogonal to the surface. In order to be able to use this PDF with the θ z -ϕ z coordinates, it is necessary to define the relationship between θ det and θ z , which is given by θ det = α + θ z . α represents the elevation angle of the scattering surface (that means, with respect to the Earth's surface). θ det , θ z and α angles are presented in figure 4. For the case of ϕ z = 0 • , α corresponds to the slope of the lateral surface. For the cases where ϕ z 0 • , it is estimated from the tangent to the tumulus surface at the scattering point. Thus, the expected background B(θ z , ϕ z ) is calculated as:

B(θ z , ϕ z ) = 5 E=0 90-α θ=85-α P (θ, α + θ z , E)φ µ (E, θ, ϕ z )dEdθ (3.3)
For the different muon flux calculations required to obtain S(θ z , ϕ z ) and B(θ z , ϕ z ), the parametrization proposed in [START_REF] Shukla | Energy and angular distributions of atmospheric muons at the Earth[END_REF] has been used, corresponding to:

φ µ (θ, E) = I 0 (n -1)E n-1 0 (E 0 + E) -n 1 + E -1 D(θ) -(n-1) (3.4) D(θ) = R 2 d 2 cos 2 θ + 2 R d + 1 - R d cosθ (3.5)
where the experimental parameters, summarized in table 1 together with other constants used in the equations, have been obtained from the fit of different experimental measurements. This parametrization provides an analytical formula for the muon flux estimation valid for low energy muons and high incident zenith angles. With all these ingredients the S/B ratio for the Apollonia tumulus has been calculated scanning the ϕ z range in 10 • steps and the corresponding θ z values for each case in 1 • steps. The results from these calculations are summarized in figure 5 as a function of θ z and the opacity , which is a more significant variable than ϕ z since the muon flux is basically independent of the azimuth angle. The main conclusion is that for all the studied directions the S/B ratio is higher than 73.9, which means that at most the 1.3 % of the detected muons have been previously scattered on the object surface. It is observed that the directions with the lowest S/B values are those with high θ z values. For these directions lower values for the signal are expected since they correspond to -8 - the most horizontal ones (where the muon flux is lower) and, due to the tumulus geometry, these are cases for which longer tumulus length is traversed. Actually, for directions with θ z lower than 85 • , the S/B ratio is always higher than 254.8, reducing the contribution of the scattered muons to the total detected to less than 0.4 %. In this region, the obtained S/B values can be considered homogeneous. Differences between directions are basically associated to the uncertainties in the muon scattering PDF. As mentioned in section 2, even if the used PDF was generated with another target material than the assumed tumulus composition, it would have a limited effect on the results. This leads to consider that the forward scattered muons on the object surface do not significantly influence the results of the muon imaging for the case of tumuli and, by extension, of other objects with similar dimensions and composition.

Volcanology: La Soufrière

The use of muon imaging for the scanning of volcanoes is another application of this technique, which implies the study of objects with larger dimensions than for archaeology. With this purpose, some projects have already performed measurements in different locations. One of them is the -9 - DIAPHANE collaboration [START_REF]Diaphane project[END_REF], which surveys La Soufrière volcano paying special attention to possible variations of the inner liquid/vapour content that can be related to the hydrothermal system dynamics. For this work, La Soufrière volcano has been taken as reference to study the impact of the forward scattered muons on the muon imaging reconstruction in volcanology.

As for the case of tumuli, volcanoes geometry can also be approximated to a truncated cone. Based on the topographic plan of La Soufrière, their dimensions correspond to a height of 460 m and a base and top radius of 840 m and 160 m respectively. These dimensions lead to a lateral surface with a slope angle of 34.1 • . In this case a homogeneous composition of standard rock has been considered, with a density of ρ = 2.5 g/cm 3 , together with 3 different detector positions corresponding to real measurement points of the DIAPHANE project. They are labelled as h-270, h-170 and h-160 respectively because of the height where they are placed. These positions are summarized in table 2 taking as reference the centre of the volcano base. Main differences among the positions rely on the distance between the detector and the volcano, going from 5 to 25 m approximately, and on the height with respect to the volcano base, which has a direct influence on the length of the volcano traversed by muons before their detection and, consequently, on the signal S(θ z ,ϕ z ) computation.

Both tumulus and volcano have been approximated to the same geometrical shape with their -10 -A c c e p t e d m a n u s c r i p t corresponding dimensions. So for volcanoes, the procedure to determine the S/B ratio for different incident directions is equivalent to the described in section 3.1 for the tumulus case. The only difference is that for this case the assumed density is ρ = 2.5 g/cm 3 instead of ρ = 2.2 g/cm 3 (corresponding to the standard soil), affecting in the opacity estimation. Nevertheless, this density variation is expected to have a reduced impact on the results as it has been estimated in section 2.

As for the tumulus case, the S/B ratio have been evaluated scanning the ϕ z range in 10 • steps and the corresponding θ z values for each case in 1 • steps. Results have also been represented with respect to θ z and the opacity . They have been summarized for the three different detector positions in figure 6. For the three cases the S/B ratio takes values significantly lower than for the tumulus, although the corresponding distributions present similar features. For example, the S/B values for directions with θ z > 85 • are again lower than for the rest of the directions. Moreover, for all detector positions, directions with low opacity (corresponding to the volcano contour) present systematically higher values of S/B than those directions pointing to the bulk of the volcano.

Focused on each detector, for the h-270 case, incident directions with high θ z have S/B ratio values in general below 1, which implies that it is possible to detect more forward scattered muons than those emerging from the volcano, significantly influencing the object density reconstruction. On the opposite side, for the directions where the opacity is smaller than 50×10 3 g/cm 2 , S/B ratio takes values higher than 5, so no more than the 17 % of the detected muons have been previously scattered. If we consider all the other directions, with θ z < 85 • and > 50 ×10 3 g/cm 2 , the S/B distribution is more homogeneous, having a mean value of 2.7 with a standard deviation of 1.4. That means that in average about 27 % of the detected muons are low energy forward scattered muons. In this case, the scattered muons have a significant impact on the volcano density reconstruction. Assuming the percentage of scattered muons constant for all the scanned directions, and estimating the uncertainty of this percentage from the standard deviation of the S/B mean value, it would imply that the reconstructed density should be corrected by multiplying it by a factor 1.4 +0. 4 -0.1 . Results for the detector positions labelled as h-170 and h-160 are similar between them. This suggests that S/B ratio has more dependence with the detector height with respect to the volcano base (170 m for h-170 and 160 m for h-160) than with the distance between the detector and the volcano (6.72 m for h-170 and 23.99 m for h-160). Since both detectors are placed lower than the h-270 case, the mean volcano length traversed by non-scattered muons is longer in this case. This leads to smaller S/B ratios mainly because lower S(θ z , ϕ z ) values. As mentioned, the features of the distribution are equivalent: smaller S/B values for θ z > 85 • and higher for low opacities. If only muon directions with θ z < 85 • and > 50 ×10 3 g/cm 2 are considered, a mean S/B value of 1.1 is obtained for both cases with a standard deviation of 0.9 and 1.0 for the h-170 and h-160 detector positions respectively. These values reveal a high influence of the low energy forward scattered muons in the overall detection (almost half of the detected muons have been previously scattered). Keeping the assumption of a constant S/B value for all the considered directions, the density correction factors are 1.9 +4.1 -0.4 for the h-170 position and 1.9 +9.1 -0.4 for the h-160 case. Summarizing, for the case of volcanoes, where the length of material to be traversed by muons is longer than for archaeology, the forward scattering of low energy muons and their further detection has a clear influence on the results of the muon imaging. The three studied scenarios and the defined geometry of the volcano as a truncated cone, reveal that the S/B ratio mainly depends on the length of material traversed by non-scattered muons, those considered for S(θ z ,ϕ z ) -11 -A c c e p t e d m a n u s c r i p t Figure 6. Distribution of the ratio between the non-scattered and scattered detected muons (defined as S/B in the text) with respect to the reconstructed zenith incident angle (θ z ) and the opacity ( ). The distribution is showed for the 3 detectors installed in La Soufrière volcano. Numbers correspond to the bin value and have been placed next to the corresponding bin to ease their reading.

-12 -A c c e p t e d m a n u s c r i p t computation. Moreover, for a fixed detector position it can be considered S/B as homogeneous for all the incident directions corresponding to the volcano bulk volume, so a global correction factor for the reconstructed density can be applied. The main source of uncertainty of the S/B ratio estimation comes from those associated to the PDF and consequently, B(θ z , ϕ z ). For this reason, as deduced for the h-270 detection position, a higher S/B mean value translates to a more accurate determination of the correction factor for the reconstructed volcano density.

Summary and discussion

At present, muon imaging is being used and proposed for an increasing number of different applications. This implies that objects with quite different dimensions can be scanned, from some tens to several hundreds of meters as typical sizes. Furthermore, the composition and density of these objects can also vary from one to other. All the experimental approaches to do muon imaging, generally known as transmission and deviation tomography respectively, rely on the direction reconstruction of the detected muons. For this reason, specially for the transmissionbased technique, muons changing their direction because of a scattering on the object surface before their detection, represent an irreducible background noise, biasing the object mean density reconstruction. An estimation of the percentage of these forward scattered muons out of all detected, would allow the estimation of correction factors to reconstruct the proper density.

Muons trajectory deviation is mainly driven by the multiple Coulomb scattering. The resulting angular distribution due to this effect is theoretically described by the Molière theory. Besides, some analytical descriptions of the process have been already performed for particular objects and compositions. Nevertheless, the large variety of objects that are currently proposed to be studied by muon tomography, requires more versatile tools to evaluate the forward scattering muons, easily adaptable to each case. With this aim a set of Monte Carlo simulations have been performed using the Geant4 framework.

These simulations provide a general evaluation of the muon forward scattering probability depending on their energy and their incident angle, being in overall agreement with theoretical estimations. They revealed that muons with energies lower than 5 GeV and incident angles above 85 • with respect to the normal direction of the surface, are almost the only muons susceptible to be scattered and then detected. The simulations results have been used as PDF to evaluate the influence of scattered muons in different scenarios. To do that, the signal to background ratio (S/B) has been defined. S(θ z ,ϕ z ) corresponds to the flux of those muons reconstructed in a direction θ z -ϕ z without any previous scattering, while B(θ z ,ϕ z ) is the muon flux of muons reconstructed with the same direction but after a previous forward scattering on the object surface. The S/B evaluation has been presented for two particular cases, corresponding to two of the applications of muon imaging: archaeology and volcanology.

Taking the muon distribution at Earth's surface proposed in [START_REF] Shukla | Energy and angular distributions of atmospheric muons at the Earth[END_REF], for the archaeological applications, the Apollonia tumulus has been considered as reference, placing the detector beside the tumulus base. S/B estimations reveals that the percentage of scattered muons detected is never higher than 1.6 %, being lower than 0.5 % if the incident zenith angle is smaller than 85 • . This leads to conclude that the influence of scattered muons for these cases can be neglected.

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

This is not the case for volcanology applications. A model based on La Soufrière volcano, already scanned inside the DIAPHANE project, has been used together with three different detector positions, corresponding to real measurement points. It has been observed a significant influence of the forward scattered muons in the measurement, which can represent up to 50 % of the detected muons, and even more for incident zenith angles higher than 85 • . S/B values can be considered homogeneous for the directions corresponding to the bulk volume of the volcano. Main differences on S/B mainly depend on the height of the detector with respect to the volcano base. Due to the volcano geometry, defined as a truncated cone, this is directly related to the muon path length along the volcano. Other features, such as the distance between the detector and the volcano, seem to have smaller influence. With the estimations and numbers obtained in this work, correction factors for the density reconstruction have been computed, taking values from 1.4 to 1.9 depending on the detector position.

Forward scattered muons represent events that are in principle not taken into account, so their detection has a direct impact on the mean density reconstruction. Nonetheless, the observed homogeneity on the S/B ratio for all the considered directions, both in the tumulus and volcano case, leads to think that these muons would not significantly affect fading the resolution of the resulting image.

All these estimations and conclusions are based on simulations of scattered muons on standard rock, which has been demonstrated to produce equivalent results than the standard soil case. In any case, changing accordingly the material composition and properties, this simulation framework can be used to evaluate the influence of forward scattered muons for further muon imaging measurements of other objects and structures.

Figure 1 .

 1 Figure 1. Schema of the defined geometry to perform the general studies of forward scattering of muons.

  has been used for all the studies. -4 -A c c e p t e d m a n u s c r i p t

Figure 2 .

 2 Figure 2. Summary plots of the results of the general study of forward scattering of muons (see text for details about the study). Top left: Difference on the zenith angle (∆θ det = θ f in det -θ ini det ) with respect to the initial muon energy (E ini µ ). Top right: Correlation between the initial and final zenith angles (θ ini det vs. θ f in det ) for all the muon energies considered. Bottom left: Difference on the azimuth angle (∆ϕ det = ϕ f in det -ϕ ini det ) with respect to the initial muon energy (E ini µ ). Bottom right: ∆ϕ det distribution for all the muon energies considered. A Gaussian distribution as predicted by Molière theory (Equations 2.1 -2.3) is observed.

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

Figure 3 .

 3 Figure 3. Correlation between the initial and final zenith angles (θ ini det vs. θ f in det ) from the general study of forward scattering of muons (see text for details about the study). The correlation plots are divided in 0.5GeV windows between 0 and 5 GeV for the initial simulated muon energy (E ini µ ). These plots, correspond to 10 8 simulated muons with incident angles between 85 • and 90 • and energies between 0 and 5 GeV (both of them homogeneously distributed). They are used as PDF for further estimations on the forward scattered muon flux.

Figure 4 .Table 1 .

 41 Figure 4. Schema showing the relationship between θ det , θ z and α angles (see text for angles definition), for the use of the muon scattering PDF, P(θ ini det , θ f in det , E ini µ ), in the B(θ z , ϕ z ) calculation.

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

 5 Figure 5. Distribution of the ratio between the non-scattered and scattered detected muons (defined as S/B in the text) with respect to the reconstructed zenith incident angle (θ z ) and the opacity ( ) for the Apollonia tumulus case.

Table 2 .

 2 Summary of the detector positions with respect to the base centre of La Soufrière volcano model (see text for model details).

Acknowledgments

Authors would like to acknowledge the financial support from the UnivEarthS Labex program of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). Data from La Soufrière volcano are part of the ANR DIAPHANE project ANR-14-ce04-0001. Part of the project has been funded by the INSU/IN2P3 TelluS and "DEFI Instrumentation aux limites" programmes.